ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.30
Committed: Mon Nov 27 21:00:14 2023 UTC (5 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.29: +40 -7 lines
Log Message:
refactor: Added number of spectral samples and wavelength splits to arguments

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: spec_rgb.c,v 2.29 2023/11/21 01:30:20 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 scolor2photopic( /* compute scotopic integral for spectral color */
350 SCOLOR scol,
351 int ncs,
352 const float wlpt[4]
353 )
354 {
355 if (ncs == 3)
356 return bright(scol);
357
358 return(spec_dot(scol, ncs, wlpt, cie_y_cumul, CIE_Y_WLMIN, CIE_Y_WLMAX));
359 }
360
361
362 double
363 scolor2scotopic( /* compute Y channel for spectral color */
364 SCOLOR scol,
365 int ncs,
366 const float wlpt[4]
367 )
368 {
369 return(spec_dot(scol, ncs, wlpt, scotopic_cumul, SCOTOPIC_WLMIN, SCOTOPIC_WLMAX));
370 }
371
372
373 double
374 scolor2melanopic( /* compute melanopic integral for spectral color */
375 SCOLOR scol,
376 int ncs,
377 const float wlpt[4]
378 )
379 {
380 return(spec_dot(scol, ncs, wlpt, melanopic_cumul, MELANOPIC_WLMIN, MELANOPIC_WLMAX));
381 }
382
383
384 double
385 scolor_photopic( /* compute scotopic integral for spectral color */
386 SCOLOR scol
387 )
388 {
389 return(scolor2photopic(scol, NCSAMP, WLPART));
390 }
391
392
393 double
394 scolor_scotopic( /* compute Y channel for spectral color */
395 SCOLOR scol
396 )
397 {
398 return(scolor2scotopic(scol, NCSAMP, WLPART));
399 }
400
401
402 double
403 scolor_melanopic( /* compute melanopic integral for spectral color */
404 SCOLOR scol
405 )
406 {
407 return(scolor2melanopic(scol, NCSAMP, WLPART));
408 }
409
410
411 void
412 cie_rgb( /* convert CIE color to standard RGB */
413 COLOR rgb,
414 COLOR xyz
415 )
416 {
417 colortrans(rgb, xyz2rgbmat, xyz);
418 clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
419 }
420
421
422 int
423 clipgamut( /* clip to gamut cube */
424 COLOR col,
425 double brt,
426 int gamut,
427 COLOR lower,
428 COLOR upper
429 )
430 {
431 int rflags = 0;
432 double brtmin, brtmax, v, vv;
433 COLOR cgry;
434 int i;
435 /* check for no check */
436 if (gamut == 0) return(0);
437 /* check brightness limits */
438 brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
439 if (gamut & CGAMUT_LOWER && brt < brtmin) {
440 copycolor(col, lower);
441 return(CGAMUT_LOWER);
442 }
443 brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
444 if (gamut & CGAMUT_UPPER && brt > brtmax) {
445 copycolor(col, upper);
446 return(CGAMUT_UPPER);
447 }
448 /* compute equivalent grey */
449 v = (brt - brtmin)/(brtmax - brtmin);
450 for (i = 0; i < 3; i++)
451 cgry[i] = v*upper[i] + (1.-v)*lower[i];
452 vv = 1.; /* check each limit */
453 for (i = 0; i < 3; i++)
454 if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
455 v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
456 if (v < vv) vv = v;
457 rflags |= CGAMUT_LOWER;
458 } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
459 v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
460 if (v < vv) vv = v;
461 rflags |= CGAMUT_UPPER;
462 }
463 if (rflags) /* desaturate to cube face */
464 for (i = 0; i < 3; i++)
465 col[i] = vv*col[i] + (1.-vv)*cgry[i];
466 return(rflags);
467 }
468
469
470 void
471 colortrans( /* convert c1 by mat and put into c2 */
472 COLOR c2,
473 COLORMAT mat,
474 COLOR c1
475 )
476 {
477 COLOR cout;
478
479 cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
480 cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
481 cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];
482
483 copycolor(c2, cout);
484 }
485
486
487 void
488 multcolormat( /* multiply m1 by m2 and put into m3 */
489 COLORMAT m3, /* m3 can be either m1 or m2 w/o harm */
490 COLORMAT m2,
491 COLORMAT m1
492 )
493 {
494 COLORMAT mt;
495 int i, j;
496
497 for (i = 0; i < 3; i++)
498 for (j = 0; j < 3; j++)
499 mt[i][j] = m1[i][0]*m2[0][j] +
500 m1[i][1]*m2[1][j] +
501 m1[i][2]*m2[2][j] ;
502 cpcolormat(m3, mt);
503 }
504
505
506 int
507 colorprimsOK( /* are color primaries reasonable? */
508 RGBPRIMS pr
509 )
510 {
511 int i, j;
512 /* check white point */
513 if ((pr[3][CIEX] <= CEPS) | (pr[3][CIEX] >= 1.-CEPS) |
514 (pr[3][CIEY] <= CEPS) | (pr[3][CIEY] >= 1.-CEPS))
515 return(0);
516 for (i = 3; i--; ) /* check for XYZ color primaries */
517 if (!CEQ(pr[i][CIEX],(i==0)) | !CEQ(pr[i][CIEY],(i==1)))
518 break;
519 if (i < 0)
520 return(-1); /* flag as XYZ color space */
521 /* check color primaries */
522 for (i = 0; i < 3; i++) {
523 if ((pr[i][CIEX] <= -2.) | (pr[i][CIEY] <= -2.))
524 return(0);
525 if ((pr[i][CIEX] >= 3.) | (pr[i][CIEY] >= 3.))
526 return(0);
527 if (pr[i][CIEX] + pr[i][CIEY] <= -2.)
528 return(0);
529 if (pr[i][CIEX] + pr[i][CIEY] >= 3.)
530 return(0);
531 }
532 for (i = 0; i < 4; i++) /* make sure space is 3-dimensional */
533 for (j = i+1; j < 4; j++)
534 if (CEQ(pr[i][CIEX],pr[j][CIEX]) &
535 CEQ(pr[i][CIEY],pr[j][CIEY]))
536 return(0);
537 return(1);
538 }
539
540
541
542 int
543 compxyz2rgbmat( /* compute conversion from CIE to RGB space */
544 COLORMAT mat,
545 RGBPRIMS pr
546 )
547 {
548 double C_rD, C_gD, C_bD;
549
550 if (pr == stdprims) { /* can use xyz2rgbmat */
551 cpcolormat(mat, xyz2rgbmat);
552 return(1);
553 }
554 if (pr == xyzprims) { /* identity */
555 memset(mat, 0, sizeof(COLORMAT));
556 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
557 return(1);
558 }
559 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
560 return(0);
561 C_rD = (1./pr[WHT][CIEY]) *
562 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
563 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
564 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
565 if (CEQ(C_rD,0.))
566 return(0);
567 C_gD = (1./pr[WHT][CIEY]) *
568 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
569 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
570 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
571 if (CEQ(C_gD,0.))
572 return(0);
573 C_bD = (1./pr[WHT][CIEY]) *
574 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
575 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
576 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
577 if (CEQ(C_bD,0.))
578 return(0);
579 mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
580 pr[BLU][CIEX]*pr[GRN][CIEY] +
581 pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
582 mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
583 pr[BLU][CIEX]*pr[GRN][CIEY] +
584 pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
585 mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
586 pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
587 mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
588 pr[BLU][CIEY]*pr[RED][CIEX] +
589 pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
590 mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
591 pr[RED][CIEX]*pr[BLU][CIEY] +
592 pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
593 mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
594 pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
595 mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
596 pr[RED][CIEY]*pr[GRN][CIEX] +
597 pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
598 mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
599 pr[GRN][CIEX]*pr[RED][CIEY] +
600 pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
601 mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
602 pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
603 return(1);
604 }
605
606
607 int
608 comprgb2xyzmat( /* compute conversion from RGB to CIE space */
609 COLORMAT mat,
610 RGBPRIMS pr
611 )
612 {
613 double C_rD, C_gD, C_bD, D;
614
615 if (pr == stdprims) { /* can use rgb2xyzmat */
616 cpcolormat(mat, rgb2xyzmat);
617 return(1);
618 }
619 if (pr == xyzprims) { /* identity */
620 memset(mat, 0, sizeof(COLORMAT));
621 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
622 return(1);
623 }
624 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
625 return(0);
626 C_rD = (1./pr[WHT][CIEY]) *
627 ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
628 pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
629 pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
630 C_gD = (1./pr[WHT][CIEY]) *
631 ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
632 pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
633 pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
634 C_bD = (1./pr[WHT][CIEY]) *
635 ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
636 pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
637 pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
638 D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
639 pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
640 pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
641 if (CEQ(D,0.))
642 return(0);
643 mat[0][0] = pr[RED][CIEX]*C_rD/D;
644 mat[0][1] = pr[GRN][CIEX]*C_gD/D;
645 mat[0][2] = pr[BLU][CIEX]*C_bD/D;
646 mat[1][0] = pr[RED][CIEY]*C_rD/D;
647 mat[1][1] = pr[GRN][CIEY]*C_gD/D;
648 mat[1][2] = pr[BLU][CIEY]*C_bD/D;
649 mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
650 mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
651 mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
652 return(1);
653 }
654
655
656 int
657 comprgb2rgbmat( /* compute conversion from RGB1 to RGB2 */
658 COLORMAT mat,
659 RGBPRIMS pr1,
660 RGBPRIMS pr2
661 )
662 {
663 COLORMAT pr1toxyz, xyztopr2;
664
665 if (pr1 == pr2) {
666 memset(mat, 0, sizeof(COLORMAT));
667 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
668 return(1);
669 }
670 if (!comprgb2xyzmat(pr1toxyz, pr1))
671 return(0);
672 if (!compxyz2rgbmat(xyztopr2, pr2))
673 return(0);
674 /* combine transforms */
675 multcolormat(mat, pr1toxyz, xyztopr2);
676 return(1);
677 }
678
679
680 int
681 compxyzWBmat( /* CIE von Kries transform from wht1 to wht2 */
682 COLORMAT mat,
683 float wht1[2],
684 float wht2[2]
685 )
686 {
687 COLOR cw1, cw2;
688 if (XYEQ(wht1,wht2)) {
689 memset(mat, 0, sizeof(COLORMAT));
690 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
691 return(1);
692 }
693 if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
694 return(0);
695 cw1[RED] = wht1[CIEX]/wht1[CIEY];
696 cw1[GRN] = 1.;
697 cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
698 colortrans(cw1, vkmat, cw1);
699 if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
700 return(0);
701 cw2[RED] = wht2[CIEX]/wht2[CIEY];
702 cw2[GRN] = 1.;
703 cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
704 colortrans(cw2, vkmat, cw2);
705 if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
706 return(0);
707 mat[0][0] = cw2[RED]/cw1[RED];
708 mat[1][1] = cw2[GRN]/cw1[GRN];
709 mat[2][2] = cw2[BLU]/cw1[BLU];
710 mat[0][1] = mat[0][2] = mat[1][0] =
711 mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
712 multcolormat(mat, vkmat, mat);
713 multcolormat(mat, mat, ivkmat);
714 return(1);
715 }
716
717
718 int
719 compxyz2rgbWBmat( /* von Kries conversion from CIE to RGB space */
720 COLORMAT mat,
721 RGBPRIMS pr
722 )
723 {
724 COLORMAT wbmat;
725
726 if (!compxyz2rgbmat(mat, pr))
727 return(0);
728 if (XYEQ(pr[WHT],xyneu))
729 return(1);
730 if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
731 return(0);
732 multcolormat(mat, wbmat, mat);
733 return(1);
734 }
735
736 int
737 comprgb2xyzWBmat( /* von Kries conversion from RGB to CIE space */
738 COLORMAT mat,
739 RGBPRIMS pr
740 )
741 {
742 COLORMAT wbmat;
743
744 if (!comprgb2xyzmat(mat, pr))
745 return(0);
746 if (XYEQ(pr[WHT],xyneu))
747 return(1);
748 if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
749 return(0);
750 multcolormat(mat, mat, wbmat);
751 return(1);
752 }
753
754 int
755 comprgb2rgbWBmat( /* von Kries conversion from RGB1 to RGB2 */
756 COLORMAT mat,
757 RGBPRIMS pr1,
758 RGBPRIMS pr2
759 )
760 {
761 COLORMAT pr1toxyz, xyztopr2, wbmat;
762
763 if (pr1 == pr2) {
764 memset(mat, 0, sizeof(COLORMAT));
765 mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
766 return(1);
767 }
768 if (!comprgb2xyzmat(pr1toxyz, pr1))
769 return(0);
770 if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
771 return(0);
772 if (!compxyz2rgbmat(xyztopr2, pr2))
773 return(0);
774 /* combine transforms */
775 multcolormat(mat, pr1toxyz, wbmat);
776 multcolormat(mat, mat, xyztopr2);
777 return(1);
778 }