ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.34
Committed: Wed Aug 14 20:05:23 2024 UTC (8 months, 3 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.33: +26 -26 lines
Log Message:
feat(rxpict): Added new C++ picture rendering tool with multi-processing and spectral output

File Contents

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