ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/spec_rgb.c
Revision: 2.27
Committed: Wed Nov 15 18:02:52 2023 UTC (5 months, 3 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.26: +301 -60 lines
Log Message: 
feat(rpict,rtrace,rcontrib,rtpict): Hyperspectral rendering (except photon map)

File Contents

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