code cleanup: unused defines, remove unused scanfill success value from BLI_scanfill_...
[blender.git] / intern / smoke / intern / spectrum.cpp
1 /*
2                 Colour Rendering of Spectra
3
4                        by John Walker
5                   http://www.fourmilab.ch/
6                   
7                  Last updated: March 9, 2003
8
9            This program is in the public domain.
10
11     For complete information about the techniques employed in
12     this program, see the World-Wide Web document:
13
14              http://www.fourmilab.ch/documents/specrend/
15              
16     The xyz_to_rgb() function, which was wrong in the original
17     version of this program, was corrected by:
18     
19             Andrew J. S. Hamilton 21 May 1999
20             Andrew.Hamilton@Colorado.EDU
21             http://casa.colorado.edu/~ajsh/
22
23     who also added the gamma correction facilities and
24     modified constrain_rgb() to work by desaturating the
25     colour by adding white.
26     
27     A program which uses these functions to plot CIE
28     "tongue" diagrams called "ppmcie" is included in
29     the Netpbm graphics toolkit:
30         http://netpbm.sourceforge.net/
31     (The program was called cietoppm in earlier
32     versions of Netpbm.)
33
34 */
35
36 #include <stdio.h>
37 #include <math.h>
38 #include "spectrum.h"
39
40 /* A colour system is defined by the CIE x and y coordinates of
41    its three primary illuminants and the x and y coordinates of
42    the white point. */
43
44 struct colourSystem {
45     const char *name;               /* Colour system name */
46     double xRed, yRed,              /* Red x, y */
47            xGreen, yGreen,          /* Green x, y */
48            xBlue, yBlue,            /* Blue x, y */
49            xWhite, yWhite,          /* White point x, y */
50            gamma;                   /* Gamma correction for system */
51 };
52
53 /* White point chromaticities. */
54
55 #if 0
56 #define IlluminantC     0.3101, 0.3162          /* For NTSC television */
57 #define IlluminantD65   0.3127, 0.3291          /* For EBU and SMPTE */
58 #endif
59 #define IlluminantE     0.33333333, 0.33333333  /* CIE equal-energy illuminant */
60
61 /*  Gamma of nonlinear correction.
62
63     See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at:
64     
65        http://www.poynton.com/ColorFAQ.html
66        http://www.poynton.com/GammaFAQ.html
67  
68 */
69
70 #define GAMMA_REC709    0               /* Rec. 709 */
71
72 static struct colourSystem
73                   /* Name                  xRed    yRed    xGreen  yGreen  xBlue  yBlue    White point        Gamma   */
74 #if 0 /* UNUSED */
75     NTSCsystem  =  { "NTSC",               0.67,   0.33,   0.21,   0.71,   0.14,   0.08,   IlluminantC,    GAMMA_REC709 },
76     EBUsystem   =  { "EBU (PAL/SECAM)",    0.64,   0.33,   0.29,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 },
77     SMPTEsystem =  { "SMPTE",              0.630,  0.340,  0.310,  0.595,  0.155,  0.070,  IlluminantD65,  GAMMA_REC709 },
78     HDTVsystem  =  { "HDTV",               0.670,  0.330,  0.210,  0.710,  0.150,  0.060,  IlluminantD65,  GAMMA_REC709 },
79 #endif
80
81     CIEsystem   =  { "CIE",                0.7355, 0.2645, 0.2658, 0.7243, 0.1669, 0.0085, IlluminantE,    GAMMA_REC709 };
82
83 #if 0 /* UNUSED */
84     Rec709system = { "CIE REC 709",        0.64,   0.33,   0.30,   0.60,   0.15,   0.06,   IlluminantD65,  GAMMA_REC709 };
85 #endif
86
87 /*                          UPVP_TO_XY
88
89     Given 1976 coordinates u', v', determine 1931 chromaticities x, y
90     
91 */
92
93 #if 0 /* UNUSED */
94 static void upvp_to_xy(double up, double vp, double *xc, double *yc)
95 {
96     *xc = (9 * up) / ((6 * up) - (16 * vp) + 12);
97     *yc = (4 * vp) / ((6 * up) - (16 * vp) + 12);
98 }
99 #endif
100
101 /*                          XY_TO_UPVP
102
103     Given 1931 chromaticities x, y, determine 1976 coordinates u', v'
104     
105 */
106
107 #if 0 /* UNUSED */
108 static void xy_to_upvp(double xc, double yc, double *up, double *vp)
109 {
110     *up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3);
111     *vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3);
112 }
113 #endif
114
115 /*                             XYZ_TO_RGB
116
117     Given an additive tricolour system CS, defined by the CIE x
118     and y chromaticities of its three primaries (z is derived
119     trivially as 1-(x+y)), and a desired chromaticity (XC, YC,
120     ZC) in CIE space, determine the contribution of each
121     primary in a linear combination which sums to the desired
122     chromaticity.  If the  requested chromaticity falls outside
123     the Maxwell  triangle (colour gamut) formed by the three
124     primaries, one of the r, g, or b weights will be negative. 
125
126     Caller can use constrain_rgb() to desaturate an
127     outside-gamut colour to the closest representation within
128     the available gamut and/or norm_rgb to normalise the RGB
129     components so the largest nonzero component has value 1.
130     
131 */
132
133 static void xyz_to_rgb(struct colourSystem *cs,
134                        double xc, double yc, double zc,
135                        double *r, double *g, double *b)
136 {
137     double xr, yr, zr, xg, yg, zg, xb, yb, zb;
138     double xw, yw, zw;
139     double rx, ry, rz, gx, gy, gz, bx, by, bz;
140     double rw, gw, bw;
141
142     xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
143     xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
144     xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);
145
146     xw = cs->xWhite;  yw = cs->yWhite;  zw = 1 - (xw + yw);
147
148     /* xyz -> rgb matrix, before scaling to white. */
149     
150     rx = (yg * zb) - (yb * zg);  ry = (xb * zg) - (xg * zb);  rz = (xg * yb) - (xb * yg);
151     gx = (yb * zr) - (yr * zb);  gy = (xr * zb) - (xb * zr);  gz = (xb * yr) - (xr * yb);
152     bx = (yr * zg) - (yg * zr);  by = (xg * zr) - (xr * zg);  bz = (xr * yg) - (xg * yr);
153
154     /* White scaling factors.
155        Dividing by yw scales the white luminance to unity, as conventional. */
156        
157     rw = ((rx * xw) + (ry * yw) + (rz * zw)) / yw;
158     gw = ((gx * xw) + (gy * yw) + (gz * zw)) / yw;
159     bw = ((bx * xw) + (by * yw) + (bz * zw)) / yw;
160
161     /* xyz -> rgb matrix, correctly scaled to white. */
162     
163     rx = rx / rw;  ry = ry / rw;  rz = rz / rw;
164     gx = gx / gw;  gy = gy / gw;  gz = gz / gw;
165     bx = bx / bw;  by = by / bw;  bz = bz / bw;
166
167     /* rgb of the desired point */
168
169     *r = (rx * xc) + (ry * yc) + (rz * zc);
170     *g = (gx * xc) + (gy * yc) + (gz * zc);
171     *b = (bx * xc) + (by * yc) + (bz * zc);
172 }
173
174 /*                            INSIDE_GAMUT
175
176      Test whether a requested colour is within the gamut
177      achievable with the primaries of the current colour
178      system.  This amounts simply to testing whether all the
179      primary weights are non-negative. */
180
181 #if 0 /* UNUSED */
182 static int inside_gamut(double r, double g, double b)
183 {
184     return (r >= 0) && (g >= 0) && (b >= 0);
185 }
186 #endif
187
188 /*                          CONSTRAIN_RGB
189
190     If the requested RGB shade contains a negative weight for
191     one of the primaries, it lies outside the colour gamut 
192     accessible from the given triple of primaries.  Desaturate
193     it by adding white, equal quantities of R, G, and B, enough
194     to make RGB all positive.  The function returns 1 if the
195     components were modified, zero otherwise.
196     
197 */
198
199 static int constrain_rgb(double *r, double *g, double *b)
200 {
201     double w;
202
203     /* Amount of white needed is w = - min(0, *r, *g, *b) */
204     
205     w = (0 < *r) ? 0 : *r;
206     w = (w < *g) ? w : *g;
207     w = (w < *b) ? w : *b;
208     w = -w;
209
210     /* Add just enough white to make r, g, b all positive. */
211     
212     if (w > 0) {
213         *r += w;  *g += w; *b += w;
214         return 1;                     /* Colour modified to fit RGB gamut */
215     }
216
217     return 0;                         /* Colour within RGB gamut */
218 }
219
220 /*                          GAMMA_CORRECT_RGB
221
222     Transform linear RGB values to nonlinear RGB values. Rec.
223     709 is ITU-R Recommendation BT. 709 (1990) ``Basic
224     Parameter Values for the HDTV Standard for the Studio and
225     for International Programme Exchange'', formerly CCIR Rec.
226     709. For details see
227     
228        http://www.poynton.com/ColorFAQ.html
229        http://www.poynton.com/GammaFAQ.html
230 */
231
232 #if 0 /* UNUSED */
233 static void gamma_correct(const struct colourSystem *cs, double *c)
234 {
235     double gamma;
236
237     gamma = cs->gamma;
238
239     if (gamma == GAMMA_REC709) {
240         /* Rec. 709 gamma correction. */
241         double cc = 0.018;
242         
243         if (*c < cc) {
244             *c *= ((1.099 * pow(cc, 0.45)) - 0.099) / cc;
245         } else {
246             *c = (1.099 * pow(*c, 0.45)) - 0.099;
247         }
248     } else {
249         /* Nonlinear colour = (Linear colour)^(1/gamma) */
250         *c = pow(*c, 1.0 / gamma);
251     }
252 }
253
254 static void gamma_correct_rgb(const struct colourSystem *cs, double *r, double *g, double *b)
255 {
256     gamma_correct(cs, r);
257     gamma_correct(cs, g);
258     gamma_correct(cs, b);
259 }
260 #endif
261
262 /*                          NORM_RGB
263
264     Normalise RGB components so the most intense (unless all
265     are zero) has a value of 1.
266     
267 */
268
269 static void norm_rgb(double *r, double *g, double *b)
270 {
271 #define Max(a, b)   (((a) > (b)) ? (a) : (b))
272     double greatest = Max(*r, Max(*g, *b));
273
274     if (greatest > 0) {
275         *r /= greatest;
276         *g /= greatest;
277         *b /= greatest;
278     }
279 #undef Max
280 }
281
282 /*                          SPECTRUM_TO_XYZ
283
284     Calculate the CIE X, Y, and Z coordinates corresponding to
285     a light source with spectral distribution given by  the
286     function SPEC_INTENS, which is called with a series of
287     wavelengths between 380 and 780 nm (the argument is 
288     expressed in meters), which returns emittance at  that
289     wavelength in arbitrary units.  The chromaticity
290     coordinates of the spectrum are returned in the x, y, and z
291     arguments which respect the identity:
292
293             x + y + z = 1.
294 */
295
296 static void spectrum_to_xyz(double (*spec_intens)(double wavelength),
297                             double *x, double *y, double *z)
298 {
299     int i;
300     double lambda, X = 0, Y = 0, Z = 0, XYZ;
301
302     /* CIE colour matching functions xBar, yBar, and zBar for
303        wavelengths from 380 through 780 nanometers, every 5
304        nanometers.  For a wavelength lambda in this range:
305
306             cie_colour_match[(lambda - 380) / 5][0] = xBar
307             cie_colour_match[(lambda - 380) / 5][1] = yBar
308             cie_colour_match[(lambda - 380) / 5][2] = zBar
309
310         To save memory, this table can be declared as floats
311         rather than doubles; (IEEE) float has enough 
312         significant bits to represent the values. It's declared
313         as a double here to avoid warnings about "conversion
314         between floating-point types" from certain persnickety
315         compilers. */
316
317     static double cie_colour_match[81][3] = {
318         {0.0014,0.0000,0.0065}, {0.0022,0.0001,0.0105}, {0.0042,0.0001,0.0201},
319         {0.0076,0.0002,0.0362}, {0.0143,0.0004,0.0679}, {0.0232,0.0006,0.1102},
320         {0.0435,0.0012,0.2074}, {0.0776,0.0022,0.3713}, {0.1344,0.0040,0.6456},
321         {0.2148,0.0073,1.0391}, {0.2839,0.0116,1.3856}, {0.3285,0.0168,1.6230},
322         {0.3483,0.0230,1.7471}, {0.3481,0.0298,1.7826}, {0.3362,0.0380,1.7721},
323         {0.3187,0.0480,1.7441}, {0.2908,0.0600,1.6692}, {0.2511,0.0739,1.5281},
324         {0.1954,0.0910,1.2876}, {0.1421,0.1126,1.0419}, {0.0956,0.1390,0.8130},
325         {0.0580,0.1693,0.6162}, {0.0320,0.2080,0.4652}, {0.0147,0.2586,0.3533},
326         {0.0049,0.3230,0.2720}, {0.0024,0.4073,0.2123}, {0.0093,0.5030,0.1582},
327         {0.0291,0.6082,0.1117}, {0.0633,0.7100,0.0782}, {0.1096,0.7932,0.0573},
328         {0.1655,0.8620,0.0422}, {0.2257,0.9149,0.0298}, {0.2904,0.9540,0.0203},
329         {0.3597,0.9803,0.0134}, {0.4334,0.9950,0.0087}, {0.5121,1.0000,0.0057},
330         {0.5945,0.9950,0.0039}, {0.6784,0.9786,0.0027}, {0.7621,0.9520,0.0021},
331         {0.8425,0.9154,0.0018}, {0.9163,0.8700,0.0017}, {0.9786,0.8163,0.0014},
332         {1.0263,0.7570,0.0011}, {1.0567,0.6949,0.0010}, {1.0622,0.6310,0.0008},
333         {1.0456,0.5668,0.0006}, {1.0026,0.5030,0.0003}, {0.9384,0.4412,0.0002},
334         {0.8544,0.3810,0.0002}, {0.7514,0.3210,0.0001}, {0.6424,0.2650,0.0000},
335         {0.5419,0.2170,0.0000}, {0.4479,0.1750,0.0000}, {0.3608,0.1382,0.0000},
336         {0.2835,0.1070,0.0000}, {0.2187,0.0816,0.0000}, {0.1649,0.0610,0.0000},
337         {0.1212,0.0446,0.0000}, {0.0874,0.0320,0.0000}, {0.0636,0.0232,0.0000},
338         {0.0468,0.0170,0.0000}, {0.0329,0.0119,0.0000}, {0.0227,0.0082,0.0000},
339         {0.0158,0.0057,0.0000}, {0.0114,0.0041,0.0000}, {0.0081,0.0029,0.0000},
340         {0.0058,0.0021,0.0000}, {0.0041,0.0015,0.0000}, {0.0029,0.0010,0.0000},
341         {0.0020,0.0007,0.0000}, {0.0014,0.0005,0.0000}, {0.0010,0.0004,0.0000},
342         {0.0007,0.0002,0.0000}, {0.0005,0.0002,0.0000}, {0.0003,0.0001,0.0000},
343         {0.0002,0.0001,0.0000}, {0.0002,0.0001,0.0000}, {0.0001,0.0000,0.0000},
344         {0.0001,0.0000,0.0000}, {0.0001,0.0000,0.0000}, {0.0000,0.0000,0.0000}
345     };
346
347     for (i = 0, lambda = 380; lambda < 780.1; i++, lambda += 5) {
348         double Me;
349
350         Me = (*spec_intens)(lambda);
351         X += Me * cie_colour_match[i][0];
352         Y += Me * cie_colour_match[i][1];
353         Z += Me * cie_colour_match[i][2];
354     }
355     XYZ = (X + Y + Z);
356     *x = X / XYZ;
357     *y = Y / XYZ;
358     *z = Z / XYZ;
359 }
360
361 /*                            BB_SPECTRUM
362
363     Calculate, by Planck's radiation law, the emittance of a black body
364     of temperature bbTemp at the given wavelength (in metres).  */
365
366 double bbTemp = 5000;                 /* Hidden temperature argument
367                                          to BB_SPECTRUM. */
368 static double bb_spectrum(double wavelength)
369 {
370     double wlm = wavelength * 1e-9;   /* Wavelength in meters */
371
372     return (3.74183e-16 * pow(wlm, -5.0)) /
373            (exp(1.4388e-2 / (wlm * bbTemp)) - 1.0);
374 }
375
376 static void xyz_to_lms(double x, double y, double z, double* l, double* m, double* s)
377 {
378         *l =  0.3897*x + 0.6890*y - 0.0787*z;
379         *m = -0.2298*x + 1.1834*y + 0.0464*z;
380         *s = z;
381 }
382
383 static void lms_to_xyz(double l, double m, double s, double* x, double *y, double* z)
384 {
385         *x =  1.9102*l - 1.1121*m + 0.2019*s;
386         *y =  0.3709*l + 0.6290*m + 0.0000*s;
387         *z = s;
388 }
389
390 void spectrum(double t1, double t2, int N, unsigned char *d)
391 {
392         int i,j,dj;
393         double X,Y,Z,R,G,B,L,M,S, Lw, Mw, Sw;
394         struct colourSystem *cs = &CIEsystem;
395
396         j = 0; dj = 1;
397         if (t1<t2) {
398                 double t = t1;
399                 t1 = t2;
400                 t2 = t;
401                 j = N-1; dj=-1;
402         }
403
404         for (i=0; i<N; i++) {
405                 bbTemp = t1 + (t2-t1)/N*i;
406
407                         // integrate blackbody radiation spectrum to XYZ
408                 spectrum_to_xyz(bb_spectrum, &X, &Y, &Z);
409
410                         // normalize highest temperature to white (in LMS system)
411                 xyz_to_lms(X,Y,Z,&L,&M,&S);
412                 if (i==0) {
413                         Lw=1/L; Mw=1/M; Sw=1/S;
414                 }
415                 L *= Lw; M *= Mw; S *= Sw;
416                 lms_to_xyz(L,M,S,&X,&Y,&Z);
417
418                         // convert to RGB
419                 xyz_to_rgb(cs, X, Y, Z, &R, &G, &B);
420                 constrain_rgb(&R, &G, &B);
421                 norm_rgb(&R, &G, &B);
422                 d[(j<<2)] = (unsigned char) ((double)R*255);
423                 d[(j<<2)+1] = (unsigned char) ((double)G*255);
424                 d[(j<<2)+2] = (unsigned char) ((double)B*255);
425                 d[(j<<2)+3] = (B>0.1)? B*255 : 0;
426                 j += dj;
427         }
428 }