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 (4 months, 2 weeks 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

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.32 static const char RCSid[] = "$Id: spec_rgb.c,v 2.31 2023/12/08 18:48:09 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 2.32 #include <math.h>
16 greg 1.1 #include "color.h"
17 greg 2.11
18     #define CEPS 1e-4 /* color epsilon */
19    
20 greg 2.26 #define CEQ(v1,v2) (((v1) <= (v2)+CEPS) & ((v2) <= (v1)+CEPS))
21 greg 1.1
22 greg 2.26 #define XYEQ(c1,c2) (CEQ((c1)[CIEX],(c2)[CIEX]) & CEQ((c1)[CIEY],(c2)[CIEY]))
23 greg 2.5
24 gwlarson 2.10
25 greg 2.11 RGBPRIMS stdprims = STDPRIMS; /* standard primary chromaticities */
26 greg 2.27 RGBPRIMS xyzprims = {1,0, 0,1, 0,0, 1.f/3.f,1.f/3.f};
27 greg 2.5
28 greg 2.8 COLOR cblack = BLKCOLOR; /* global black color */
29     COLOR cwhite = WHTCOLOR; /* global white color */
30    
31 greg 2.11 float xyneu[2] = {1./3., 1./3.}; /* neutral xy chromaticities */
32    
33 greg 1.1 /*
34 greg 2.27 * 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 greg 1.1 */
38    
39 greg 2.27 #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 greg 1.1 };
201    
202 greg 2.11 COLORMAT xyz2rgbmat = { /* XYZ to RGB (no white balance) */
203 greg 2.3 {(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 greg 1.2 };
213 greg 1.1
214 greg 2.11 COLORMAT rgb2xyzmat = { /* RGB to XYZ (no white balance) */
215 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},
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 greg 1.2
222 greg 2.11 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 greg 1.2
234 greg 2.4
235 greg 2.11 void
236 greg 2.27 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 greg 2.15 spec_rgb( /* compute RGB color from spectral range */
266     COLOR col,
267     int s,
268     int e
269     )
270 greg 1.1 {
271     COLOR ciecolor;
272    
273     spec_cie(ciecolor, s, e);
274 greg 2.31 colortrans(col, xyz2rgbmat, ciecolor);
275 greg 1.1 }
276    
277    
278 greg 2.27 static double
279     spec_dot( /* spectrum dot-product with cumulative observer */
280     SCOLOR scol,
281     int ncs,
282 greg 2.29 const float wlpt[4],
283 greg 2.27 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 greg 2.28 wl1 = (int)(ncs==3 ? wlpt[2-n] : wlpt[3] + (n+1)*wlstp);
302 greg 2.27 if (wl1 >= wlmax) {
303 greg 2.28 sum += (65535 - cumul[wl0-wlmin]) * scol[ncs-1-n];
304 greg 2.27 break;
305     }
306 greg 2.28 sum += (cumul[wl1-wlmin] - cumul[wl0-wlmin]) * scol[ncs-1-n++];
307 greg 2.27 }
308     return(sum * (1./65535.));
309     }
310    
311    
312 greg 2.11 void
313 greg 2.27 scolor2cie( /* accurate conversion from spectrum to XYZ */
314     COLOR col,
315     SCOLOR scol,
316     int ncs,
317 greg 2.29 const float wlpt[4]
318 greg 2.27 )
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 greg 2.29 const float wlpt[4]
336 greg 2.15 )
337 greg 1.1 {
338 greg 2.27 COLOR ciecolor;
339    
340     if (ncs == 3) { /* not really a spectrum, so we punt... */
341     copycolor(col, scol);
342 greg 1.2 return;
343     }
344 greg 2.27 scolor2cie(ciecolor, scol, ncs, wlpt);
345     cie_rgb(col, ciecolor);
346     }
347    
348 greg 1.1
349 greg 2.27 double
350 greg 2.30 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 greg 2.27 scolor_photopic( /* compute scotopic integral for spectral color */
387     SCOLOR scol
388     )
389     {
390 greg 2.30 return(scolor2photopic(scol, NCSAMP, WLPART));
391 greg 2.27 }
392    
393    
394     double
395     scolor_scotopic( /* compute Y channel for spectral color */
396     SCOLOR scol
397     )
398     {
399 greg 2.30 return(scolor2scotopic(scol, NCSAMP, WLPART));
400 greg 2.27 }
401 greg 1.1
402    
403 greg 2.27 double
404     scolor_melanopic( /* compute melanopic integral for spectral color */
405     SCOLOR scol
406     )
407     {
408 greg 2.30 return(scolor2melanopic(scol, NCSAMP, WLPART));
409 greg 1.1 }
410    
411    
412 greg 2.11 void
413 greg 2.32 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 greg 2.15 cie_rgb( /* convert CIE color to standard RGB */
441     COLOR rgb,
442     COLOR xyz
443     )
444 greg 2.8 {
445     colortrans(rgb, xyz2rgbmat, xyz);
446     clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
447     }
448    
449    
450     int
451 greg 2.15 clipgamut( /* clip to gamut cube */
452     COLOR col,
453     double brt,
454     int gamut,
455     COLOR lower,
456     COLOR upper
457     )
458 greg 2.8 {
459     int rflags = 0;
460     double brtmin, brtmax, v, vv;
461     COLOR cgry;
462 greg 2.25 int i;
463 greg 2.8 /* 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 greg 2.14 v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
484 greg 2.8 if (v < vv) vv = v;
485     rflags |= CGAMUT_LOWER;
486     } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
487 greg 2.14 v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
488 greg 2.8 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 greg 2.11 void
499 greg 2.15 colortrans( /* convert c1 by mat and put into c2 */
500 greg 2.25 COLOR c2,
501     COLORMAT mat,
502     COLOR c1
503 greg 2.15 )
504 greg 1.1 {
505 greg 2.9 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 greg 2.4 }
513    
514    
515 greg 2.11 void
516 greg 2.15 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 greg 2.4 {
522 greg 2.5 COLORMAT mt;
523 greg 2.25 int i, j;
524 greg 2.4
525     for (i = 0; i < 3; i++)
526 greg 2.5 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 greg 2.15 int
535     colorprimsOK( /* are color primaries reasonable? */
536     RGBPRIMS pr
537     )
538     {
539 greg 2.20 int i, j;
540 greg 2.26 /* 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 greg 2.22 for (i = 0; i < 3; i++) {
551 greg 2.24 if ((pr[i][CIEX] <= -2.) | (pr[i][CIEY] <= -2.))
552 greg 2.15 return(0);
553 greg 2.24 if ((pr[i][CIEX] >= 3.) | (pr[i][CIEY] >= 3.))
554 greg 2.17 return(0);
555 greg 2.24 if (pr[i][CIEX] + pr[i][CIEY] <= -2.)
556 greg 2.21 return(0);
557 greg 2.24 if (pr[i][CIEX] + pr[i][CIEY] >= 3.)
558 greg 2.15 return(0);
559     }
560 greg 2.26 for (i = 0; i < 4; i++) /* make sure space is 3-dimensional */
561 greg 2.20 for (j = i+1; j < 4; j++)
562 greg 2.26 if (CEQ(pr[i][CIEX],pr[j][CIEX]) &
563 greg 2.20 CEQ(pr[i][CIEY],pr[j][CIEY]))
564     return(0);
565 greg 2.15 return(1);
566     }
567    
568    
569    
570     int
571     compxyz2rgbmat( /* compute conversion from CIE to RGB space */
572     COLORMAT mat,
573 greg 2.25 RGBPRIMS pr
574 greg 2.15 )
575 greg 2.5 {
576     double C_rD, C_gD, C_bD;
577    
578 greg 2.6 if (pr == stdprims) { /* can use xyz2rgbmat */
579     cpcolormat(mat, xyz2rgbmat);
580 greg 2.16 return(1);
581 greg 2.6 }
582 greg 2.27 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 greg 2.15 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
588     return(0);
589 greg 2.5 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 greg 2.15 if (CEQ(C_rD,0.))
594     return(0);
595 greg 2.5 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 greg 2.15 if (CEQ(C_gD,0.))
600     return(0);
601 greg 2.5 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 greg 2.15 if (CEQ(C_bD,0.))
606     return(0);
607 greg 2.5 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 greg 2.15 return(1);
632 greg 2.5 }
633    
634    
635 greg 2.15 int
636     comprgb2xyzmat( /* compute conversion from RGB to CIE space */
637     COLORMAT mat,
638 greg 2.25 RGBPRIMS pr
639 greg 2.15 )
640 greg 2.5 {
641     double C_rD, C_gD, C_bD, D;
642    
643 greg 2.6 if (pr == stdprims) { /* can use rgb2xyzmat */
644     cpcolormat(mat, rgb2xyzmat);
645 greg 2.16 return(1);
646 greg 2.6 }
647 greg 2.27 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 greg 2.15 if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
653     return(0);
654 greg 2.5 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 greg 2.15 if (CEQ(D,0.))
670     return(0);
671 greg 2.5 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 greg 2.15 return(1);
681 greg 2.5 }
682    
683    
684 greg 2.15 int
685     comprgb2rgbmat( /* compute conversion from RGB1 to RGB2 */
686     COLORMAT mat,
687     RGBPRIMS pr1,
688     RGBPRIMS pr2
689     )
690 greg 2.5 {
691     COLORMAT pr1toxyz, xyztopr2;
692    
693 greg 2.6 if (pr1 == pr2) {
694 greg 2.27 memset(mat, 0, sizeof(COLORMAT));
695     mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
696 greg 2.16 return(1);
697 greg 2.6 }
698 greg 2.15 if (!comprgb2xyzmat(pr1toxyz, pr1))
699     return(0);
700     if (!compxyz2rgbmat(xyztopr2, pr2))
701     return(0);
702 greg 2.5 /* combine transforms */
703     multcolormat(mat, pr1toxyz, xyztopr2);
704 greg 2.15 return(1);
705 greg 2.11 }
706    
707    
708 greg 2.15 int
709     compxyzWBmat( /* CIE von Kries transform from wht1 to wht2 */
710     COLORMAT mat,
711     float wht1[2],
712     float wht2[2]
713     )
714 greg 2.11 {
715     COLOR cw1, cw2;
716     if (XYEQ(wht1,wht2)) {
717 greg 2.27 memset(mat, 0, sizeof(COLORMAT));
718     mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
719 greg 2.15 return(1);
720 greg 2.11 }
721 greg 2.15 if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
722     return(0);
723 greg 2.11 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 greg 2.15 if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
728     return(0);
729 greg 2.11 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 greg 2.15 if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
734     return(0);
735 greg 2.11 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 greg 2.15 return(1);
743 greg 2.11 }
744    
745    
746 greg 2.15 int
747     compxyz2rgbWBmat( /* von Kries conversion from CIE to RGB space */
748     COLORMAT mat,
749     RGBPRIMS pr
750     )
751 greg 2.11 {
752     COLORMAT wbmat;
753    
754 greg 2.15 if (!compxyz2rgbmat(mat, pr))
755     return(0);
756 greg 2.11 if (XYEQ(pr[WHT],xyneu))
757 greg 2.15 return(1);
758     if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
759     return(0);
760 greg 2.11 multcolormat(mat, wbmat, mat);
761 greg 2.15 return(1);
762 greg 2.11 }
763    
764 greg 2.15 int
765     comprgb2xyzWBmat( /* von Kries conversion from RGB to CIE space */
766     COLORMAT mat,
767     RGBPRIMS pr
768     )
769 greg 2.11 {
770     COLORMAT wbmat;
771    
772 greg 2.15 if (!comprgb2xyzmat(mat, pr))
773     return(0);
774 greg 2.11 if (XYEQ(pr[WHT],xyneu))
775 greg 2.15 return(1);
776     if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
777     return(0);
778 greg 2.11 multcolormat(mat, mat, wbmat);
779 greg 2.15 return(1);
780 greg 2.11 }
781    
782 greg 2.15 int
783     comprgb2rgbWBmat( /* von Kries conversion from RGB1 to RGB2 */
784     COLORMAT mat,
785     RGBPRIMS pr1,
786     RGBPRIMS pr2
787     )
788 greg 2.11 {
789     COLORMAT pr1toxyz, xyztopr2, wbmat;
790    
791     if (pr1 == pr2) {
792 greg 2.27 memset(mat, 0, sizeof(COLORMAT));
793     mat[0][0] = mat[1][1] = mat[2][2] = 1.f;
794 greg 2.15 return(1);
795 greg 2.11 }
796 greg 2.15 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 greg 2.11 /* combine transforms */
803     multcolormat(mat, pr1toxyz, wbmat);
804     multcolormat(mat, mat, xyztopr2);
805 greg 2.15 return(1);
806 greg 1.1 }