ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.32
Committed: Wed Jan 17 17:36:20 2024 UTC (3 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.31: +29 -1 lines
Log Message:
feat: Added convertscolorcol() to generate better RGB result

File Contents

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