ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.29
Committed: Tue Nov 21 01:30:20 2023 UTC (5 months, 2 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.28: +4 -4 lines
Log Message:
feat(rmtxop): Added handling of Radiance spectral pictures as matrices

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: spec_rgb.c,v 2.28 2023/11/17 20:02:07 greg Exp $";
3 #endif
4 /*
5 * Convert colors and spectral ranges.
6 * Added von Kries white-balance calculations 10/01 (GW).
7 *
8 * Externals declared in color.h
9 */
10
11 #include "copyright.h"
12
13 #include <stdio.h>
14 #include <string.h>
15 #include "color.h"
16
17 #define CEPS 1e-4 /* color epsilon */
18
19 #define CEQ(v1,v2) (((v1) <= (v2)+CEPS) & ((v2) <= (v1)+CEPS))
20
21 #define XYEQ(c1,c2) (CEQ((c1)[CIEX],(c2)[CIEX]) & CEQ((c1)[CIEY],(c2)[CIEY]))
22
23
24 RGBPRIMS stdprims = STDPRIMS; /* standard primary chromaticities */
25 RGBPRIMS xyzprims = {1,0, 0,1, 0,0, 1.f/3.f,1.f/3.f};
26
27 COLOR cblack = BLKCOLOR; /* global black color */
28 COLOR cwhite = WHTCOLOR; /* global white color */
29
30 float xyneu[2] = {1./3., 1./3.}; /* neutral xy chromaticities */
31
32 /*
33 * The tables below encode the CIE 1931 2° standard observer
34 * functions for X, Y, and Z. These tables are cumulative, as are
35 * ones that follow.
36 */
37
38 #define CIE_X_WLMIN 362
39 #define CIE_X_WLMAX 774
40
41 static const unsigned short cie_x_cumul[CIE_X_WLMAX-CIE_X_WLMIN+1] = {
42 0,1,1,1,1,1,1,2,2,2,3,3,3,4,5,5,6,7,8,9,10,11,12,14,16,18,20,23,26,29,
43 33,37,41,47,53,60,68,77,86,97,108,121,135,151,170,190,214,241,271,304,
44 342,385,432,486,545,612,686,768,860,961,1073,1195,1326,1467,1618,
45 1776,1943,2117,2298,2485,2677,2875,3076,3281,3489,3700,3912,4125,4340,4555,
46 4769,4983,5197,5409,5620,5830,6038,6244,6448,6651,6851,7049,7245,7437,
47 7627,7813,7995,8173,8347,8517,8682,8842,8996,9143,9284,9418,9545,9665,
48 9778,9884,9984,10077,10164,10245,10321,10390,10454,10513,10566,10615,
49 10659,10698,10734,10766,10794,10819,10842,10861,10879,10893,10906,10917,
50 10926,10933,10939,10944,10948,10951,10953,10955,10957,10958,10960,10961,
51 10964,10967,10971,10977,10984,10994,11006,11020,11038,11060,11085,11114,
52 11148,11187,11231,11280,11335,11396,11464,11537,11618,11705,11799,11901,
53 12010,12126,12249,12380,12518,12664,12818,12980,13150,13328,13515,13709,
54 13913,14124,14345,14574,14813,15060,15317,15582,15858,16142,16437,16741,
55 17055,17379,17713,18057,18412,18776,19151,19536,19931,20337,20753,21180,
56 21616,22063,22520,22988,23465,23953,24450,24957,25474,26000,26535,27080,
57 27633,28195,28765,29343,29929,30522,31122,31729,32342,32961,33585,34215,
58 34849,35487,36129,36775,37423,38073,38724,39375,40027,40679,41329,41978,
59 42625,43270,43911,44548,45181,45808,46429,47044,47652,48253,48845,49430,
60 50005,50571,51128,51674,52209,52733,53245,53745,54232,54706,55167,55614,
61 56048,56468,56876,57269,57651,58019,58376,58720,59052,59373,59681,59979,
62 60265,60539,60803,61056,61298,61529,61750,61962,62163,62355,62538,62712,
63 62877,63034,63183,63325,63459,63586,63706,63820,63927,64028,64123,64213,
64 64297,64376,64451,64520,64586,64647,64704,64758,64808,64855,64899,64941,
65 64980,65016,65051,65083,65114,65143,65169,65194,65218,65240,65260,65278,
66 65296,65312,65327,65341,65354,65366,65377,65387,65397,65406,65415,65423,
67 65430,65437,65444,65450,65455,65461,65466,65470,65475,65479,65483,65486,
68 65489,65493,65495,65498,65501,65503,65505,65507,65509,65511,65513,65514,
69 65516,65517,65518,65519,65520,65521,65522,65523,65524,65525,65526,65526,
70 65527,65527,65528,65528,65529,65529,65530,65530,65530,65531,65531,65531,
71 65532,65532,65532,65532,65532,65533,65533,65533,65533,65533,65533,65533,
72 65533,65534,65534,65534,65534,65534,65534,65534,65534,65534,65534,65534,
73 65534,65534,65534,65534,65535
74 };
75
76 #define CIE_Y_WLMIN 386
77 #define CIE_Y_WLMAX 760
78
79 static const unsigned short cie_y_cumul[CIE_Y_WLMAX-CIE_Y_WLMIN+1] = {
80 0,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,3,4,4,5,5,6,7,8,8,10,11,12,14,
81 15,17,19,22,25,28,31,35,40,45,50,56,63,70,78,86,95,105,115,126,
82 138,150,164,178,193,208,225,242,260,280,300,321,343,367,391,
83 417,443,472,501,532,564,598,633,670,708,748,790,833,879,926,
84 975,1027,1080,1136,1194,1255,1318,1384,1453,1525,1601,1679,1761,
85 1846,1935,2027,2123,2223,2327,2435,2548,2665,2787,2915,3048,3187,
86 3332,3484,3643,3808,3981,4162,4352,4550,4757,4975,5203,5442,5691,
87 5952,6225,6509,6805,7114,7435,7769,8116,8476,8849,9235,9634,10045,
88 10469,10904,11351,11808,12275,12752,13239,13734,14239,14752,15273,
89 15801,16337,16880,17429,17984,18546,19112,19684,20260,20841,21426,
90 22015,22608,23203,23802,24403,25007,25612,26220,26828,27439,28050,
91 28662,29275,29888,30501,31114,31727,32340,32951,33561,34170,34777,
92 35382,35985,36585,37182,37777,38368,38955,39539,40119,40695,41266,
93 41832,42393,42950,43501,44046,44586,45119,45646,46167,46682,47189,
94 47690,48183,48670,49149,49621,50085,50542,50991,51432,51866,52293,
95 52711,53122,53524,53919,54306,54685,55057,55420,55775,56123,56463,
96 56795,57119,57435,57743,58044,58337,58623,58901,59172,59435,59691,
97 59939,60180,60414,60640,60859,61070,61274,61471,61661,61844,62019,
98 62188,62351,62507,62657,62802,62940,63073,63201,63323,63441,63553,
99 63660,63763,63861,63954,64043,64128,64208,64285,64358,64427,64493,
100 64555,64614,64670,64723,64773,64820,64865,64907,64947,64984,65019,
101 65052,65083,65113,65140,65166,65190,65212,65233,65253,65271,65288,
102 65304,65319,65334,65347,65360,65371,65383,65393,65403,65412,65420,
103 65428,65435,65442,65449,65454,65460,65465,65470,65474,65478,65482,
104 65485,65488,65492,65494,65497,65500,65502,65504,65506,65508,65510,
105 65512,65513,65515,65516,65517,65519,65520,65521,65522,65523,65523,
106 65524,65525,65526,65526,65527,65527,65528,65528,65529,65529,65530,
107 65530,65530,65531,65531,65531,65532,65532,65532,65532,65532,65533,
108 65533,65533,65533,65533,65533,65533,65534,65534,65534,65534,65534,
109 65534,65534,65534,65534,65534,65534,65534,65534,65534,65534,65534,65535
110 };
111
112 #define CIE_Z_WLMIN 359
113 #define CIE_Z_WLMAX 624
114
115 static const unsigned short cie_z_cumul[CIE_Z_WLMAX-CIE_Z_WLMIN+1] = {
116 0,1,1,2,2,3,4,5,6,7,8,9,11,12,14,16,19,22,25,28,32,37,41,47,52,59,66,75,84,
117 95,107,121,136,154,173,196,221,250,283,321,362,408,458,512,573,640,717,804,
118 903,1015,1142,1285,1446,1627,1830,2058,2313,2598,2917,3272,3668,4108,4597,
119 5135,5723,6360,7044,7773,8544,9356,10205,11090,12006,12952,13923,14918,
120 15934,16967,18016,19076,20148,21227,22312,23401,24492,25585,26678,27771,
121 28861,29950,31037,32121,33202,34281,35355,36424,37487,38542,39588,40624,
122 41647,42657,43652,44631,45590,46527,47438,48321,49173,49993,50783,51542,
123 52270,52968,53636,54275,54885,55465,56018,56543,57042,57514,57961,58384,
124 58784,59162,59519,59856,60175,60477,60762,61032,61287,61529,61757,61974,
125 62179,62374,62559,62734,62901,63059,63211,63355,63492,63622,63745,63862,
126 63971,64075,64172,64263,64347,64427,64500,64569,64632,64692,64747,64798,
127 64846,64891,64933,64973,65010,65045,65078,65109,65139,65166,65192,65216,
128 65239,65260,65280,65298,65315,65331,65345,65359,65371,65383,65393,65403,
129 65412,65420,65428,65435,65441,65447,65452,65457,65462,65466,65470,65473,
130 65476,65479,65482,65485,65487,65489,65491,65493,65495,65497,65498,65500,
131 65501,65503,65504,65505,65506,65508,65509,65510,65511,65512,65513,65514,
132 65515,65516,65517,65518,65519,65520,65520,65521,65522,65523,65523,65524,
133 65525,65525,65526,65527,65527,65528,65528,65529,65529,65530,65530,65531,
134 65531,65531,65532,65532,65532,65532,65533,65533,65533,65533,65533,65533,
135 65534,65534,65534,65534,65534,65534,65534,65534,65534,65535
136 };
137
138 #define SCOTOPIC_WLMIN 382
139 #define SCOTOPIC_WLMAX 684
140
141 static const unsigned short scotopic_cumul[SCOTOPIC_WLMAX-SCOTOPIC_WLMIN+1] = {
142 0, 1, 1, 2, 3, 4, 5, 6, 7, 9, 10, 12, 15, 17, 20, 24, 28, 33, 38, 44,
143 52, 60, 69, 80, 93, 107, 123, 142, 163, 186, 213, 242, 275, 312, 353,
144 398, 448, 502, 562, 627, 698, 775, 859, 949, 1046, 1150, 1261, 1380,
145 1507, 1642, 1785, 1936, 2096, 2265, 2442, 2628, 2823, 3027, 3239, 3461,
146 3691, 3930, 4178, 4435, 4700, 4974, 5257, 5548, 5847, 6154, 6469, 6793,
147 7123, 7462, 7809, 8162, 8524, 8892, 9268, 9651, 10041, 10438, 10843, 11254,
148 11673, 12099, 12532, 12973, 13422, 13878, 14342, 14814, 15293, 15780, 16276,
149 16779, 17290, 17809, 18336, 18872, 19415, 19967, 20526, 21093, 21667, 22249,
150 22839, 23436, 24039, 24649, 25266, 25890, 26519, 27154, 27795, 28441, 29092,
151 29747, 30405, 31068, 31734, 32402, 33074, 33747, 34420, 35095, 35771, 36446,
152 37119, 37793, 38464, 39132, 39798, 40460, 41118, 41772, 42421, 43064, 43701,
153 44332, 44957, 45575, 46185, 46787, 47381, 47967, 48543, 49110, 49668, 50215,
154 50753, 51280, 51797, 52302, 52797, 53281, 53754, 54215, 54665, 55104, 55531,
155 55947, 56352, 56744, 57125, 57495, 57853, 58200, 58536, 58860, 59174, 59477,
156 59769, 60051, 60322, 60583, 60834, 61075, 61306, 61528, 61741, 61944, 62139,
157 62326, 62504, 62674, 62836, 62991, 63138, 63278, 63412, 63538, 63659, 63773,
158 63881, 63983, 64080, 64172, 64259, 64340, 64418, 64490, 64559, 64623, 64684,
159 64742, 64795, 64846, 64893, 64937, 64978, 65017, 65053, 65087, 65119, 65149,
160 65176, 65202, 65226, 65248, 65269, 65289, 65307, 65323, 65339, 65353, 65367,
161 65379, 65391, 65402, 65412, 65421, 65430, 65438, 65445, 65452, 65458, 65464,
162 65469, 65474, 65479, 65483, 65487, 65491, 65494, 65497, 65500, 65503, 65505,
163 65507, 65509, 65511, 65513, 65514, 65516, 65517, 65519, 65520, 65521, 65522,
164 65523, 65524, 65525, 65525, 65526, 65527, 65527, 65528, 65528, 65529, 65529,
165 65529, 65530, 65530, 65530, 65531, 65531, 65531, 65531, 65532, 65532, 65532,
166 65532, 65532, 65533, 65533, 65533, 65533, 65533, 65533, 65533, 65533, 65533,
167 65533, 65534, 65534, 65534, 65534, 65534, 65534, 65534, 65534, 65535
168 };
169
170 #define MELANOPIC_WLMIN 380
171 #define MELANOPIC_WLMAX 651
172
173 static const unsigned short melanopic_cumul[MELANOPIC_WLMAX-MELANOPIC_WLMIN+1] = {
174 0, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 17, 20, 23, 27, 32, 37, 42, 49, 56,
175 65, 75, 86, 99, 114, 131, 150, 173, 200, 230, 265, 303, 347, 395, 448,
176 508, 574, 650, 734, 828, 931, 1041, 1158, 1283, 1415, 1554, 1703, 1862,
177 2031, 2211, 2400, 2600, 2809, 3028, 3257, 3497, 3748, 4012, 4288, 4576,
178 4876, 5187, 5509, 5842, 6185, 6540, 6905, 7283, 7673, 8075, 8489, 8915,
179 9351, 9799, 10259, 10729, 11211, 11705, 12211, 12729, 13258, 13799,
180 14351, 14915, 15491, 16078, 16676, 17286, 17908, 18541, 19184, 19837,
181 20498, 21168, 21846, 22532, 23226, 23928, 24637, 25353, 26075, 26802,
182 27532, 28267, 29005, 29746, 30488, 31233, 31979, 32726, 33473, 34220,
183 34966, 35711, 36455, 37196, 37935, 38671, 39403, 40129, 40851, 41568,
184 42279, 42983, 43680, 44369, 45050, 45724, 46388, 47043, 47688, 48322,
185 48945, 49557, 50156, 50743, 51317, 51879, 52428, 52964, 53487, 53997,
186 54493, 54976, 55445, 55901, 56343, 56771, 57186, 57587, 57976, 58351,
187 58712, 59061, 59397, 59720, 60031, 60330, 60616, 60891, 61154, 61405,
188 61646, 61875, 62094, 62303, 62501, 62690, 62870, 63040, 63201, 63354,
189 63498, 63634, 63763, 63884, 63998, 64105, 64206, 64300, 64389, 64472,
190 64549, 64622, 64689, 64753, 64811, 64866, 64917, 64964, 65008, 65049,
191 65086, 65121, 65154, 65184, 65211, 65237, 65260, 65282, 65302, 65321,
192 65338, 65354, 65368, 65381, 65394, 65405, 65415, 65425, 65434, 65442,
193 65449, 65456, 65463, 65468, 65474, 65479, 65483, 65487, 65491, 65494,
194 65498, 65501, 65503, 65506, 65508, 65510, 65512, 65514, 65515, 65517,
195 65518, 65520, 65521, 65522, 65523, 65524, 65525, 65525, 65526, 65527,
196 65527, 65528, 65528, 65529, 65529, 65530, 65530, 65530, 65531, 65531,
197 65531, 65531, 65532, 65532, 65532, 65532, 65532, 65533, 65533, 65533,
198 65533, 65533, 65533, 65533, 65533, 65533, 65534, 65534, 65534, 65535
199 };
200
201 COLORMAT xyz2rgbmat = { /* XYZ to RGB (no white balance) */
202 {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
203 (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
204 (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
205 {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
206 (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
207 (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
208 {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
209 (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
210 (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
211 };
212
213 COLORMAT rgb2xyzmat = { /* RGB to XYZ (no white balance) */
214 {CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
215 {CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
216 {(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
217 (1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
218 (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
219 };
220
221 COLORMAT vkmat = { /* Sharp primary matrix */
222 { 1.2694, -0.0988, -0.1706},
223 {-0.8364, 1.8006, 0.0357},
224 { 0.0297, -0.0315, 1.0018}
225 };
226
227 COLORMAT ivkmat = { /* inverse Sharp primary matrix */
228 { 0.8156, 0.0472, 0.1372},
229 { 0.3791, 0.5769, 0.0440},
230 {-0.0123, 0.0167, 0.9955}
231 };
232
233
234 void
235 spec_cie( /* compute XYZ color from a spectral range */
236 COLOR col, /* returned color */
237 int s, /* starting and ending wavelengths */
238 int e
239 )
240 {
241 if (s > e) { int i = s; s = e; e = i; }
242
243 if ((s >= CIE_X_WLMAX) | (e <= CIE_X_WLMIN))
244 col[CIEX] = 0;
245 else
246 col[CIEX] = (cie_x_cumul[e>=CIE_X_WLMAX ? CIE_X_WLMAX-CIE_X_WLMIN : e-CIE_X_WLMIN] -
247 cie_x_cumul[(s>CIE_X_WLMIN)*(s-CIE_X_WLMIN)])*(1./65535.);
248
249 if ((s >= CIE_Y_WLMAX) | (e <= CIE_Y_WLMIN))
250 col[CIEY] = 0;
251 else
252 col[CIEY] = (cie_y_cumul[e>=CIE_Y_WLMAX ? CIE_Y_WLMAX-CIE_Y_WLMIN : e-CIE_Y_WLMIN] -
253 cie_y_cumul[(s>CIE_Y_WLMIN)*(s-CIE_Y_WLMIN)])*(1./65535.);
254
255 if ((s >= CIE_Z_WLMAX) | (e <= CIE_Z_WLMIN))
256 col[CIEZ] = 0;
257 else
258 col[CIEZ] = (cie_z_cumul[e>=CIE_Z_WLMAX ? CIE_Z_WLMAX-CIE_Z_WLMIN : e-CIE_Z_WLMIN] -
259 cie_z_cumul[(s>CIE_Z_WLMIN)*(s-CIE_Z_WLMIN)])*(1./65535.);
260 }
261
262
263 void
264 spec_rgb( /* compute RGB color from spectral range */
265 COLOR col,
266 int s,
267 int e
268 )
269 {
270 COLOR ciecolor;
271
272 spec_cie(ciecolor, s, e);
273 cie_rgb(col, ciecolor);
274 }
275
276
277 static double
278 spec_dot( /* spectrum dot-product with cumulative observer */
279 SCOLOR scol,
280 int ncs,
281 const float wlpt[4],
282 const unsigned short cumul[],
283 int wlmin,
284 int wlmax
285 )
286 {
287 const double wlstp = (wlpt[0] - wlpt[3])/(double)ncs;
288 int wl1 = (int)wlpt[3];
289 int n = 0;
290 double sum = 0;
291
292 if (wl1 < wlmin) {
293 n += (int)((wlmin - wl1)/wlstp);
294 wl1 = wlmin;
295 } else if (wl1 >= wlmax)
296 return(.0);
297
298 while (n < ncs) {
299 int wl0 = wl1;
300 wl1 = (int)(ncs==3 ? wlpt[2-n] : wlpt[3] + (n+1)*wlstp);
301 if (wl1 >= wlmax) {
302 sum += (65535 - cumul[wl0-wlmin]) * scol[ncs-1-n];
303 break;
304 }
305 sum += (cumul[wl1-wlmin] - cumul[wl0-wlmin]) * scol[ncs-1-n++];
306 }
307 return(sum * (1./65535.));
308 }
309
310
311 void
312 scolor2cie( /* accurate conversion from spectrum to XYZ */
313 COLOR col,
314 SCOLOR scol,
315 int ncs,
316 const float wlpt[4]
317 )
318 {
319 if (ncs == 3) { /* not a spectrum! */
320 rgb_cie(col, scol);
321 return;
322 }
323 col[CIEX] = spec_dot(scol, ncs, wlpt, cie_x_cumul, CIE_X_WLMIN, CIE_X_WLMAX);
324 col[CIEY] = spec_dot(scol, ncs, wlpt, cie_y_cumul, CIE_Y_WLMIN, CIE_Y_WLMAX);
325 col[CIEZ] = spec_dot(scol, ncs, wlpt, cie_z_cumul, CIE_Z_WLMIN, CIE_Z_WLMAX);
326 }
327
328
329 void
330 scolor2rgb( /* accurate conversion from spectrum to RGB */
331 COLOR col,
332 SCOLOR scol,
333 int ncs,
334 const float wlpt[4]
335 )
336 {
337 COLOR ciecolor;
338
339 if (ncs == 3) { /* not really a spectrum, so we punt... */
340 copycolor(col, scol);
341 return;
342 }
343 scolor2cie(ciecolor, scol, ncs, wlpt);
344 cie_rgb(col, ciecolor);
345 }
346
347
348 double
349 scolor_photopic( /* compute scotopic integral for spectral color */
350 SCOLOR scol
351 )
352 {
353 if (NCSAMP == 3)
354 return bright(scol);
355
356 return(spec_dot(scol, NCSAMP, WLPART, cie_y_cumul, CIE_Y_WLMIN, CIE_Y_WLMAX));
357 }
358
359
360 double
361 scolor_scotopic( /* compute Y channel for spectral color */
362 SCOLOR scol
363 )
364 {
365 return(spec_dot(scol, NCSAMP, WLPART, scotopic_cumul, SCOTOPIC_WLMIN, SCOTOPIC_WLMAX));
366 }
367
368
369 double
370 scolor_melanopic( /* compute melanopic integral for spectral color */
371 SCOLOR scol
372 )
373 {
374 return(spec_dot(scol, NCSAMP, WLPART, melanopic_cumul, MELANOPIC_WLMIN, MELANOPIC_WLMAX));
375 }
376
377
378 void
379 cie_rgb( /* convert CIE color to standard RGB */
380 COLOR rgb,
381 COLOR xyz
382 )
383 {
384 colortrans(rgb, xyz2rgbmat, xyz);
385 clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
386 }
387
388
389 int
390 clipgamut( /* clip to gamut cube */
391 COLOR col,
392 double brt,
393 int gamut,
394 COLOR lower,
395 COLOR upper
396 )
397 {
398 int rflags = 0;
399 double brtmin, brtmax, v, vv;
400 COLOR cgry;
401 int i;
402 /* check for no check */
403 if (gamut == 0) return(0);
404 /* check brightness limits */
405 brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
406 if (gamut & CGAMUT_LOWER && brt < brtmin) {
407 copycolor(col, lower);
408 return(CGAMUT_LOWER);
409 }
410 brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
411 if (gamut & CGAMUT_UPPER && brt > brtmax) {
412 copycolor(col, upper);
413 return(CGAMUT_UPPER);
414 }
415 /* compute equivalent grey */
416 v = (brt - brtmin)/(brtmax - brtmin);
417 for (i = 0; i < 3; i++)
418 cgry[i] = v*upper[i] + (1.-v)*lower[i];
419 vv = 1.; /* check each limit */
420 for (i = 0; i < 3; i++)
421 if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
422 v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
423 if (v < vv) vv = v;
424 rflags |= CGAMUT_LOWER;
425 } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
426 v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
427 if (v < vv) vv = v;
428 rflags |= CGAMUT_UPPER;
429 }
430 if (rflags) /* desaturate to cube face */
431 for (i = 0; i < 3; i++)
432 col[i] = vv*col[i] + (1.-vv)*cgry[i];
433 return(rflags);
434 }
435
436
437 void
438 colortrans( /* convert c1 by mat and put into c2 */
439 COLOR c2,
440 COLORMAT mat,
441 COLOR c1
442 )
443 {
444 COLOR cout;
445
446 cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
447 cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
448 cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
449
450 copycolor(c2, cout);
451 }
452
453
454 void
455 multcolormat( /* multiply m1 by m2 and put into m3 */
456 COLORMAT m3, /* m3 can be either m1 or m2 w/o harm */
457 COLORMAT m2,
458 COLORMAT m1
459 )
460 {
461 COLORMAT mt;
462 int i, j;
463
464 for (i = 0; i < 3; i++)
465 for (j = 0; j < 3; j++)
466 mt[i][j] = m1[i][0]*m2[0][j] +
467 m1[i][1]*m2[1][j] +
468 m1[i][2]*m2[2][j] ;
469 cpcolormat(m3, mt);
470 }
471
472
473 int
474 colorprimsOK( /* are color primaries reasonable? */
475 RGBPRIMS pr
476 )
477 {
478 int i, j;
479 /* check white point */
480 if ((pr[3][CIEX] <= CEPS) | (pr[3][CIEX] >= 1.-CEPS) |
481 (pr[3][CIEY] <= CEPS) | (pr[3][CIEY] >= 1.-CEPS))
482 return(0);
483 for (i = 3; i--; ) /* check for XYZ color primaries */
484 if (!CEQ(pr[i][CIEX],(i==0)) | !CEQ(pr[i][CIEY],(i==1)))
485 break;
486 if (i < 0)
487 return(-1); /* flag as XYZ color space */
488 /* check color primaries */
489 for (i = 0; i < 3; i++) {
490 if ((pr[i][CIEX] <= -2.) | (pr[i][CIEY] <= -2.))
491 return(0);
492 if ((pr[i][CIEX] >= 3.) | (pr[i][CIEY] >= 3.))
493 return(0);
494 if (pr[i][CIEX] + pr[i][CIEY] <= -2.)
495 return(0);
496 if (pr[i][CIEX] + pr[i][CIEY] >= 3.)
497 return(0);
498 }
499 for (i = 0; i < 4; i++) /* make sure space is 3-dimensional */
500 for (j = i+1; j < 4; j++)
501 if (CEQ(pr[i][CIEX],pr[j][CIEX]) &
502 CEQ(pr[i][CIEY],pr[j][CIEY]))
503 return(0);
504 return(1);
505 }
506
507
508
509 int
510 compxyz2rgbmat( /* compute conversion from CIE to RGB space */
511 COLORMAT mat,
512 RGBPRIMS pr
513 )
514 {
515 double C_rD, C_gD, C_bD;
516
517 if (pr == stdprims) { /* can use xyz2rgbmat */
518 cpcolormat(mat, xyz2rgbmat);
519 return(1);
520 }
521 if (pr == xyzprims) { /* identity */
522 memset(mat, 0, sizeof(COLORMAT));
523 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
524 return(1);
525 }
526 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
527 return(0);
528 C_rD = (1./pr[WHT][CIEY]) *
529 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
530 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
531 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
532 if (CEQ(C_rD,0.))
533 return(0);
534 C_gD = (1./pr[WHT][CIEY]) *
535 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
536 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
537 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
538 if (CEQ(C_gD,0.))
539 return(0);
540 C_bD = (1./pr[WHT][CIEY]) *
541 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
542 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
543 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
544 if (CEQ(C_bD,0.))
545 return(0);
546 mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
547 pr[BLU][CIEX]*pr[GRN][CIEY] +
548 pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
549 mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
550 pr[BLU][CIEX]*pr[GRN][CIEY] +
551 pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
552 mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
553 pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
554 mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
555 pr[BLU][CIEY]*pr[RED][CIEX] +
556 pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
557 mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
558 pr[RED][CIEX]*pr[BLU][CIEY] +
559 pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
560 mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
561 pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
562 mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
563 pr[RED][CIEY]*pr[GRN][CIEX] +
564 pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
565 mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
566 pr[GRN][CIEX]*pr[RED][CIEY] +
567 pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
568 mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
569 pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
570 return(1);
571 }
572
573
574 int
575 comprgb2xyzmat( /* compute conversion from RGB to CIE space */
576 COLORMAT mat,
577 RGBPRIMS pr
578 )
579 {
580 double C_rD, C_gD, C_bD, D;
581
582 if (pr == stdprims) { /* can use rgb2xyzmat */
583 cpcolormat(mat, rgb2xyzmat);
584 return(1);
585 }
586 if (pr == xyzprims) { /* identity */
587 memset(mat, 0, sizeof(COLORMAT));
588 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
589 return(1);
590 }
591 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
592 return(0);
593 C_rD = (1./pr[WHT][CIEY]) *
594 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
595 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
596 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
597 C_gD = (1./pr[WHT][CIEY]) *
598 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
599 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
600 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
601 C_bD = (1./pr[WHT][CIEY]) *
602 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
603 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
604 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
605 D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
606 pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
607 pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
608 if (CEQ(D,0.))
609 return(0);
610 mat[0][0] = pr[RED][CIEX]*C_rD/D;
611 mat[0][1] = pr[GRN][CIEX]*C_gD/D;
612 mat[0][2] = pr[BLU][CIEX]*C_bD/D;
613 mat[1][0] = pr[RED][CIEY]*C_rD/D;
614 mat[1][1] = pr[GRN][CIEY]*C_gD/D;
615 mat[1][2] = pr[BLU][CIEY]*C_bD/D;
616 mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
617 mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
618 mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
619 return(1);
620 }
621
622
623 int
624 comprgb2rgbmat( /* compute conversion from RGB1 to RGB2 */
625 COLORMAT mat,
626 RGBPRIMS pr1,
627 RGBPRIMS pr2
628 )
629 {
630 COLORMAT pr1toxyz, xyztopr2;
631
632 if (pr1 == pr2) {
633 memset(mat, 0, sizeof(COLORMAT));
634 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
635 return(1);
636 }
637 if (!comprgb2xyzmat(pr1toxyz, pr1))
638 return(0);
639 if (!compxyz2rgbmat(xyztopr2, pr2))
640 return(0);
641 /* combine transforms */
642 multcolormat(mat, pr1toxyz, xyztopr2);
643 return(1);
644 }
645
646
647 int
648 compxyzWBmat( /* CIE von Kries transform from wht1 to wht2 */
649 COLORMAT mat,
650 float wht1[2],
651 float wht2[2]
652 )
653 {
654 COLOR cw1, cw2;
655 if (XYEQ(wht1,wht2)) {
656 memset(mat, 0, sizeof(COLORMAT));
657 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
658 return(1);
659 }
660 if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
661 return(0);
662 cw1[RED] = wht1[CIEX]/wht1[CIEY];
663 cw1[GRN] = 1.;
664 cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
665 colortrans(cw1, vkmat, cw1);
666 if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
667 return(0);
668 cw2[RED] = wht2[CIEX]/wht2[CIEY];
669 cw2[GRN] = 1.;
670 cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
671 colortrans(cw2, vkmat, cw2);
672 if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
673 return(0);
674 mat[0][0] = cw2[RED]/cw1[RED];
675 mat[1][1] = cw2[GRN]/cw1[GRN];
676 mat[2][2] = cw2[BLU]/cw1[BLU];
677 mat[0][1] = mat[0][2] = mat[1][0] =
678 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
679 multcolormat(mat, vkmat, mat);
680 multcolormat(mat, mat, ivkmat);
681 return(1);
682 }
683
684
685 int
686 compxyz2rgbWBmat( /* von Kries conversion from CIE to RGB space */
687 COLORMAT mat,
688 RGBPRIMS pr
689 )
690 {
691 COLORMAT wbmat;
692
693 if (!compxyz2rgbmat(mat, pr))
694 return(0);
695 if (XYEQ(pr[WHT],xyneu))
696 return(1);
697 if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
698 return(0);
699 multcolormat(mat, wbmat, mat);
700 return(1);
701 }
702
703 int
704 comprgb2xyzWBmat( /* von Kries conversion from RGB to CIE space */
705 COLORMAT mat,
706 RGBPRIMS pr
707 )
708 {
709 COLORMAT wbmat;
710
711 if (!comprgb2xyzmat(mat, pr))
712 return(0);
713 if (XYEQ(pr[WHT],xyneu))
714 return(1);
715 if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
716 return(0);
717 multcolormat(mat, mat, wbmat);
718 return(1);
719 }
720
721 int
722 comprgb2rgbWBmat( /* von Kries conversion from RGB1 to RGB2 */
723 COLORMAT mat,
724 RGBPRIMS pr1,
725 RGBPRIMS pr2
726 )
727 {
728 COLORMAT pr1toxyz, xyztopr2, wbmat;
729
730 if (pr1 == pr2) {
731 memset(mat, 0, sizeof(COLORMAT));
732 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
733 return(1);
734 }
735 if (!comprgb2xyzmat(pr1toxyz, pr1))
736 return(0);
737 if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
738 return(0);
739 if (!compxyz2rgbmat(xyztopr2, pr2))
740 return(0);
741 /* combine transforms */
742 multcolormat(mat, pr1toxyz, wbmat);
743 multcolormat(mat, mat, xyztopr2);
744 return(1);
745 }