ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/common/spec_rgb.c
Revision: 2.30
Committed: Mon Nov 27 21:00:14 2023 UTC (23 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

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.30 static const char RCSid[] = "$Id: spec_rgb.c,v 2.29 2023/11/21 01:30:20 greg Exp $";
3 greg 1.1 #endif
4     /*
5     * Convert colors and spectral ranges.
6 greg 2.11 * Added von Kries white-balance calculations 10/01 (GW).
7     *
8     * Externals declared in color.h
9     */
10    
11 greg 2.12 #include "copyright.h"
12 greg 1.1
13 greg 2.13 #include <stdio.h>
14     #include <string.h>
15 greg 1.1 #include "color.h"
16 greg 2.11
17     #define CEPS 1e-4 /* color epsilon */
18    
19 greg 2.26 #define CEQ(v1,v2) (((v1) <= (v2)+CEPS) & ((v2) <= (v1)+CEPS))
20 greg 1.1
21 greg 2.26 #define XYEQ(c1,c2) (CEQ((c1)[CIEX],(c2)[CIEX]) & CEQ((c1)[CIEY],(c2)[CIEY]))
22 greg 2.5
23 gwlarson 2.10
24 greg 2.11 RGBPRIMS stdprims = STDPRIMS; /* standard primary chromaticities */
25 greg 2.27 RGBPRIMS xyzprims = {1,0, 0,1, 0,0, 1.f/3.f,1.f/3.f};
26 greg 2.5
27 greg 2.8 COLOR cblack = BLKCOLOR; /* global black color */
28     COLOR cwhite = WHTCOLOR; /* global white color */
29    
30 greg 2.11 float xyneu[2] = {1./3., 1./3.}; /* neutral xy chromaticities */
31    
32 greg 1.1 /*
33 greg 2.27 * 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 greg 1.1 */
37    
38 greg 2.27 #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 greg 1.1 };
200    
201 greg 2.11 COLORMAT xyz2rgbmat = { /* XYZ to RGB (no white balance) */
202 greg 2.3 {(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 greg 1.2 };
212 greg 1.1
213 greg 2.11 COLORMAT rgb2xyzmat = { /* RGB to XYZ (no white balance) */
214 greg 2.4 {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 greg 1.2
221 greg 2.11 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 greg 1.2
233 greg 2.4
234 greg 2.11 void
235 greg 2.27 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 greg 2.15 spec_rgb( /* compute RGB color from spectral range */
265     COLOR col,
266     int s,
267     int e
268     )
269 greg 1.1 {
270     COLOR ciecolor;
271    
272     spec_cie(ciecolor, s, e);
273     cie_rgb(col, ciecolor);
274     }
275    
276    
277 greg 2.27 static double
278     spec_dot( /* spectrum dot-product with cumulative observer */
279     SCOLOR scol,
280     int ncs,
281 greg 2.29 const float wlpt[4],
282 greg 2.27 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 greg 2.28 wl1 = (int)(ncs==3 ? wlpt[2-n] : wlpt[3] + (n+1)*wlstp);
301 greg 2.27 if (wl1 >= wlmax) {
302 greg 2.28 sum += (65535 - cumul[wl0-wlmin]) * scol[ncs-1-n];
303 greg 2.27 break;
304     }
305 greg 2.28 sum += (cumul[wl1-wlmin] - cumul[wl0-wlmin]) * scol[ncs-1-n++];
306 greg 2.27 }
307     return(sum * (1./65535.));
308     }
309    
310    
311 greg 2.11 void
312 greg 2.27 scolor2cie( /* accurate conversion from spectrum to XYZ */
313     COLOR col,
314     SCOLOR scol,
315     int ncs,
316 greg 2.29 const float wlpt[4]
317 greg 2.27 )
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 greg 2.29 const float wlpt[4]
335 greg 2.15 )
336 greg 1.1 {
337 greg 2.27 COLOR ciecolor;
338    
339     if (ncs == 3) { /* not really a spectrum, so we punt... */
340     copycolor(col, scol);
341 greg 1.2 return;
342     }
343 greg 2.27 scolor2cie(ciecolor, scol, ncs, wlpt);
344     cie_rgb(col, ciecolor);
345     }
346    
347 greg 1.1
348 greg 2.27 double
349 greg 2.30 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 greg 2.27 scolor_photopic( /* compute scotopic integral for spectral color */
386     SCOLOR scol
387     )
388     {
389 greg 2.30 return(scolor2photopic(scol, NCSAMP, WLPART));
390 greg 2.27 }
391    
392    
393     double
394     scolor_scotopic( /* compute Y channel for spectral color */
395     SCOLOR scol
396     )
397     {
398 greg 2.30 return(scolor2scotopic(scol, NCSAMP, WLPART));
399 greg 2.27 }
400 greg 1.1
401    
402 greg 2.27 double
403     scolor_melanopic( /* compute melanopic integral for spectral color */
404     SCOLOR scol
405     )
406     {
407 greg 2.30 return(scolor2melanopic(scol, NCSAMP, WLPART));
408 greg 1.1 }
409    
410    
411 greg 2.11 void
412 greg 2.15 cie_rgb( /* convert CIE color to standard RGB */
413     COLOR rgb,
414     COLOR xyz
415     )
416 greg 2.8 {
417     colortrans(rgb, xyz2rgbmat, xyz);
418     clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
419     }
420    
421    
422     int
423 greg 2.15 clipgamut( /* clip to gamut cube */
424     COLOR col,
425     double brt,
426     int gamut,
427     COLOR lower,
428     COLOR upper
429     )
430 greg 2.8 {
431     int rflags = 0;
432     double brtmin, brtmax, v, vv;
433     COLOR cgry;
434 greg 2.25 int i;
435 greg 2.8 /* 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 greg 2.14 v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
456 greg 2.8 if (v < vv) vv = v;
457     rflags |= CGAMUT_LOWER;
458     } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
459 greg 2.14 v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
460 greg 2.8 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 greg 2.11 void
471 greg 2.15 colortrans( /* convert c1 by mat and put into c2 */
472 greg 2.25 COLOR c2,
473     COLORMAT mat,
474     COLOR c1
475 greg 2.15 )
476 greg 1.1 {
477 greg 2.9 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 greg 2.4 }
485    
486    
487 greg 2.11 void
488 greg 2.15 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 greg 2.4 {
494 greg 2.5 COLORMAT mt;
495 greg 2.25 int i, j;
496 greg 2.4
497     for (i = 0; i < 3; i++)
498 greg 2.5 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 greg 2.15 int
507     colorprimsOK( /* are color primaries reasonable? */
508     RGBPRIMS pr
509     )
510     {
511 greg 2.20 int i, j;
512 greg 2.26 /* 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 greg 2.22 for (i = 0; i < 3; i++) {
523 greg 2.24 if ((pr[i][CIEX] <= -2.) | (pr[i][CIEY] <= -2.))
524 greg 2.15 return(0);
525 greg 2.24 if ((pr[i][CIEX] >= 3.) | (pr[i][CIEY] >= 3.))
526 greg 2.17 return(0);
527 greg 2.24 if (pr[i][CIEX] + pr[i][CIEY] <= -2.)
528 greg 2.21 return(0);
529 greg 2.24 if (pr[i][CIEX] + pr[i][CIEY] >= 3.)
530 greg 2.15 return(0);
531     }
532 greg 2.26 for (i = 0; i < 4; i++) /* make sure space is 3-dimensional */
533 greg 2.20 for (j = i+1; j < 4; j++)
534 greg 2.26 if (CEQ(pr[i][CIEX],pr[j][CIEX]) &
535 greg 2.20 CEQ(pr[i][CIEY],pr[j][CIEY]))
536     return(0);
537 greg 2.15 return(1);
538     }
539    
540    
541    
542     int
543     compxyz2rgbmat( /* compute conversion from CIE to RGB space */
544     COLORMAT mat,
545 greg 2.25 RGBPRIMS pr
546 greg 2.15 )
547 greg 2.5 {
548     double C_rD, C_gD, C_bD;
549    
550 greg 2.6 if (pr == stdprims) { /* can use xyz2rgbmat */
551     cpcolormat(mat, xyz2rgbmat);
552 greg 2.16 return(1);
553 greg 2.6 }
554 greg 2.27 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 greg 2.15 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
560     return(0);
561 greg 2.5 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 greg 2.15 if (CEQ(C_rD,0.))
566     return(0);
567 greg 2.5 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 greg 2.15 if (CEQ(C_gD,0.))
572     return(0);
573 greg 2.5 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 greg 2.15 if (CEQ(C_bD,0.))
578     return(0);
579 greg 2.5 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 greg 2.15 return(1);
604 greg 2.5 }
605    
606    
607 greg 2.15 int
608     comprgb2xyzmat( /* compute conversion from RGB to CIE space */
609     COLORMAT mat,
610 greg 2.25 RGBPRIMS pr
611 greg 2.15 )
612 greg 2.5 {
613     double C_rD, C_gD, C_bD, D;
614    
615 greg 2.6 if (pr == stdprims) { /* can use rgb2xyzmat */
616     cpcolormat(mat, rgb2xyzmat);
617 greg 2.16 return(1);
618 greg 2.6 }
619 greg 2.27 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 greg 2.15 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
625     return(0);
626 greg 2.5 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 greg 2.15 if (CEQ(D,0.))
642     return(0);
643 greg 2.5 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 greg 2.15 return(1);
653 greg 2.5 }
654    
655    
656 greg 2.15 int
657     comprgb2rgbmat( /* compute conversion from RGB1 to RGB2 */
658     COLORMAT mat,
659     RGBPRIMS pr1,
660     RGBPRIMS pr2
661     )
662 greg 2.5 {
663     COLORMAT pr1toxyz, xyztopr2;
664    
665 greg 2.6 if (pr1 == pr2) {
666 greg 2.27 memset(mat, 0, sizeof(COLORMAT));
667     mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
668 greg 2.16 return(1);
669 greg 2.6 }
670 greg 2.15 if (!comprgb2xyzmat(pr1toxyz, pr1))
671     return(0);
672     if (!compxyz2rgbmat(xyztopr2, pr2))
673     return(0);
674 greg 2.5 /* combine transforms */
675     multcolormat(mat, pr1toxyz, xyztopr2);
676 greg 2.15 return(1);
677 greg 2.11 }
678    
679    
680 greg 2.15 int
681     compxyzWBmat( /* CIE von Kries transform from wht1 to wht2 */
682     COLORMAT mat,
683     float wht1[2],
684     float wht2[2]
685     )
686 greg 2.11 {
687     COLOR cw1, cw2;
688     if (XYEQ(wht1,wht2)) {
689 greg 2.27 memset(mat, 0, sizeof(COLORMAT));
690     mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
691 greg 2.15 return(1);
692 greg 2.11 }
693 greg 2.15 if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
694     return(0);
695 greg 2.11 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 greg 2.15 if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
700     return(0);
701 greg 2.11 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 greg 2.15 if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
706     return(0);
707 greg 2.11 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 greg 2.15 return(1);
715 greg 2.11 }
716    
717    
718 greg 2.15 int
719     compxyz2rgbWBmat( /* von Kries conversion from CIE to RGB space */
720     COLORMAT mat,
721     RGBPRIMS pr
722     )
723 greg 2.11 {
724     COLORMAT wbmat;
725    
726 greg 2.15 if (!compxyz2rgbmat(mat, pr))
727     return(0);
728 greg 2.11 if (XYEQ(pr[WHT],xyneu))
729 greg 2.15 return(1);
730     if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
731     return(0);
732 greg 2.11 multcolormat(mat, wbmat, mat);
733 greg 2.15 return(1);
734 greg 2.11 }
735    
736 greg 2.15 int
737     comprgb2xyzWBmat( /* von Kries conversion from RGB to CIE space */
738     COLORMAT mat,
739     RGBPRIMS pr
740     )
741 greg 2.11 {
742     COLORMAT wbmat;
743    
744 greg 2.15 if (!comprgb2xyzmat(mat, pr))
745     return(0);
746 greg 2.11 if (XYEQ(pr[WHT],xyneu))
747 greg 2.15 return(1);
748     if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
749     return(0);
750 greg 2.11 multcolormat(mat, mat, wbmat);
751 greg 2.15 return(1);
752 greg 2.11 }
753    
754 greg 2.15 int
755     comprgb2rgbWBmat( /* von Kries conversion from RGB1 to RGB2 */
756     COLORMAT mat,
757     RGBPRIMS pr1,
758     RGBPRIMS pr2
759     )
760 greg 2.11 {
761     COLORMAT pr1toxyz, xyztopr2, wbmat;
762    
763     if (pr1 == pr2) {
764 greg 2.27 memset(mat, 0, sizeof(COLORMAT));
765     mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
766 greg 2.15 return(1);
767 greg 2.11 }
768 greg 2.15 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 greg 2.11 /* combine transforms */
775     multcolormat(mat, pr1toxyz, wbmat);
776     multcolormat(mat, mat, xyztopr2);
777 greg 2.15 return(1);
778 greg 1.1 }