; Ported to IDL from "Colour Rendering of Spectra" by John Walker ; http://www.fourmilab.ch/documents/specrend/specrend.c pro setup_color, FORCE=force dontrecurse=1 @./TOOLS/common_color.inc if KEYWORD_SET(force) or not KEYWORD_SET(colorSystem) then begin ; A colour system is defined by the CIE x and y coordinates of ; its three primary illuminants and the x and y coordinates of ; the white point. colorSystem = {colorSystem, name:'', xy:DBLARR(2,4), gamma:0.D} ; White point chromaticities. IlluminantC = [0.3101D, 0.3162D] ; For NTSC television IlluminantD65 = [0.3127D, 0.3291D] ; For EBU and SMPTE IlluminantE = [0.33333333D, 0.33333333D] ; CIE equal-energy illuminant ; Gamma of nonlinear correction. ; See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at: ; http://www.poynton.com/ColorFAQ.html ; http://www.poynton.com/GammaFAQ.html GAMMA_REC709 = 0 ; Rec. 709 NTSCsystem = {colorSystem, name:"NTSC", xy:[[ 0.67D, 0.33D], [ 0.21D, 0.71D], [ 0.14D, 0.08D], [IlluminantC ]], gamma:GAMMA_REC709} EBUsystem = {colorSystem, name:"EBU (PAL/SECAM)",xy:[[ 0.64D, 0.33D], [ 0.29D, 0.60D], [ 0.15D, 0.06D], [IlluminantD65]], gamma:GAMMA_REC709} SMPTEsystem = {colorSystem, name:"SMPTE", xy:[[ 0.630D, 0.340D], [ 0.310D, 0.595D], [ 0.155D, 0.070D], [IlluminantD65]], gamma:GAMMA_REC709} HDTVsystem = {colorSystem, name:"HDTV", xy:[[ 0.670D, 0.330D], [ 0.210D, 0.710D], [ 0.150D, 0.060D], [IlluminantD65]], gamma:GAMMA_REC709} CIEsystem = {colorSystem, name:"CIE", xy:[[0.7355D, 0.2645D], [0.2658D, 0.7243D], [0.1669D, 0.0085D], [IlluminantE ]], gamma:GAMMA_REC709} Rec709system = {colorSystem, name:"CIE REC 709", xy:[[ 0.64D, 0.33D], [ 0.30D, 0.60D], [ 0.15D, 0.06D], [IlluminantD65]], gamma:GAMMA_REC709} ; CIE colour matching functions xBar, yBar, and zBar for ; wavelengths from 380 through 780 nanometers, every 5 ; nanometers. For a wavelength lambda in this range: ; cie_colour_match[(lambda - 380) / 5][0] = xBar ; cie_colour_match[(lambda - 380) / 5][1] = yBar ; cie_colour_match[(lambda - 380) / 5][2] = zBar ; To save memory, this table can be declared as floats ; rather than doubles; (IEEE) float has enough ; significant bits to represent the values. It's declared ; as a double here to avoid warnings about "conversion ; between floating-point types" from certain persnickety ; compilers. cie_colour_match = [ $ [0.0014,0.0000,0.0065], [0.0022,0.0001,0.0105], [0.0042,0.0001,0.0201],$ [0.0076,0.0002,0.0362], [0.0143,0.0004,0.0679], [0.0232,0.0006,0.1102],$ [0.0435,0.0012,0.2074], [0.0776,0.0022,0.3713], [0.1344,0.0040,0.6456],$ [0.2148,0.0073,1.0391], [0.2839,0.0116,1.3856], [0.3285,0.0168,1.6230],$ [0.3483,0.0230,1.7471], [0.3481,0.0298,1.7826], [0.3362,0.0380,1.7721],$ [0.3187,0.0480,1.7441], [0.2908,0.0600,1.6692], [0.2511,0.0739,1.5281],$ [0.1954,0.0910,1.2876], [0.1421,0.1126,1.0419], [0.0956,0.1390,0.8130],$ [0.0580,0.1693,0.6162], [0.0320,0.2080,0.4652], [0.0147,0.2586,0.3533],$ [0.0049,0.3230,0.2720], [0.0024,0.4073,0.2123], [0.0093,0.5030,0.1582],$ [0.0291,0.6082,0.1117], [0.0633,0.7100,0.0782], [0.1096,0.7932,0.0573],$ [0.1655,0.8620,0.0422], [0.2257,0.9149,0.0298], [0.2904,0.9540,0.0203],$ [0.3597,0.9803,0.0134], [0.4334,0.9950,0.0087], [0.5121,1.0000,0.0057],$ [0.5945,0.9950,0.0039], [0.6784,0.9786,0.0027], [0.7621,0.9520,0.0021],$ [0.8425,0.9154,0.0018], [0.9163,0.8700,0.0017], [0.9786,0.8163,0.0014],$ [1.0263,0.7570,0.0011], [1.0567,0.6949,0.0010], [1.0622,0.6310,0.0008],$ [1.0456,0.5668,0.0006], [1.0026,0.5030,0.0003], [0.9384,0.4412,0.0002],$ [0.8544,0.3810,0.0002], [0.7514,0.3210,0.0001], [0.6424,0.2650,0.0000],$ [0.5419,0.2170,0.0000], [0.4479,0.1750,0.0000], [0.3608,0.1382,0.0000],$ [0.2835,0.1070,0.0000], [0.2187,0.0816,0.0000], [0.1649,0.0610,0.0000],$ [0.1212,0.0446,0.0000], [0.0874,0.0320,0.0000], [0.0636,0.0232,0.0000],$ [0.0468,0.0170,0.0000], [0.0329,0.0119,0.0000], [0.0227,0.0082,0.0000],$ [0.0158,0.0057,0.0000], [0.0114,0.0041,0.0000], [0.0081,0.0029,0.0000],$ [0.0058,0.0021,0.0000], [0.0041,0.0015,0.0000], [0.0029,0.0010,0.0000],$ [0.0020,0.0007,0.0000], [0.0014,0.0005,0.0000], [0.0010,0.0004,0.0000],$ [0.0007,0.0002,0.0000], [0.0005,0.0002,0.0000], [0.0003,0.0001,0.0000],$ [0.0002,0.0001,0.0000], [0.0002,0.0001,0.0000], [0.0001,0.0000,0.0000],$ [0.0001,0.0000,0.0000], [0.0001,0.0000,0.0000], [0.0000,0.0000,0.0000] $ ] endif end ; UPVP_TO_XY ; Given 1976 coordinates u', v', determine 1931 chromaticities x, y pro upvp_to_xy, up, vp, XC=xc, YC=yc xc = (9 * up) / ((6 * up) - (16 * vp) + 12) yc = (4 * vp) / ((6 * up) - (16 * vp) + 12) end ; XY_TO_UPVP ; Given 1931 chromaticities x, y, determine 1976 coordinates u', v' pro xy_to_upvp, xc, yc, UP=up, VP=vp up = (4 * xc) / ((-2 * xc) + (12 * yc) + 3) vp = (9 * yc) / ((-2 * xc) + (12 * yc) + 3) end ; XYZ_TO_RGB ; Given an additive tricolour system CS, defined by the CIE x ; and y chromaticities of its three primaries (z is derived ; trivially as 1-(x+y)), and a desired chromaticity (XC, YC, ; ZC) in CIE space, determine the contribution of each ; primary in a linear combination which sums to the desired ; chromaticity. If the requested chromaticity falls outside ; the Maxwell triangle (colour gamut) formed by the three ; primaries, one of the r, g, or b weights will be negative. ; ; Caller can use constrain_rgb() to desaturate an ; outside-gamut colour to the closest representation within ; the available gamut and/or norm_rgb to normalise the RGB ; components so the largest nonzero component has value 1. function xyz_to_rgb, cs, c xyzrgb = [cs.xy, 1 - (cs.xy[0,*]+cs.xy[1,*])] ; find the matrix of rgb coordinates of x, y, and z ; xyz -> rgb matrix, before scaling to white. rgbxyz = MATRIX_POWER(xyzrgb[*,0:2], -1) w = rgbxyz ## transpose(xyzrgb[*,3]) / xyzrgb[1,3] ; White scaling factors. Dividing by yw scales the white luminance to unity, as conventional. rgbxyz /= [w,w,w] ; xyz -> rgb matrix, correctly scaled to white. return, rgbxyz ## transpose(c) ; rgb of the desired point end ; INSIDE_GAMUT ; Test whether a requested colour is within the gamut ; achievable with the primaries of the current colour ; system. This amounts simply to testing whether all the ; primary weights are non-negative. function inside_gamut, rgb return, min(rgb) lt 0 end ; CONSTRAIN_RGB ; If the requested RGB shade contains a negative weight for ; one of the primaries, it lies outside the colour gamut ; accessible from the given triple of primaries. Desaturate ; it by adding white, equal quantities of R, G, and B, enough ; to make RGB all positive. The function returns 1 if the ; components were modified, zero otherwise. function constrain_rgb, rgb return, rgb + (-min(rgb) > 0) ; Amount of white needed is w = - min(0, *r, *g, *b) end ; GAMMA_CORRECT_RGB ; Transform linear RGB values to nonlinear RGB values. Rec. ; 709 is ITU-R Recommendation BT. 709 (1990) ``Basic ; Parameter Values for the HDTV Standard for the Studio and ; for International Programme Exchange'', formerly CCIR Rec. ; 709. For details see ; http://www.poynton.com/ColorFAQ.html ; http://www.poynton.com/GammaFAQ.html pro gamma_correct, cs, c @./TOOLS/common_color.inc if cs.gamma eq GAMMA_REC709 then begin ; Rec. 709 gamma correction. cc = 0.018D linear = where(c lt cc, COMPLEMENT=power) c[linear] *= ((1.099D * cc^0.45D) - 0.099D) / cc c[power] = (1.099D * c[power]^0.45D) - 0.099D endif else begin ; Nonlinear colour = (Linear colour)^(1/gamma) c = c^(1.0D / gamma) endelse end ; NORM_RGB ; Normalise RGB components so the most intense (unless all ; are zero) has a value of 1. function norm_rgb, rgb m = max(rgb) return, rgb / ((m eq 0.D) ? 1.D : m) end ; SPECTRUM_TO_XYZ ; Calculate the CIE X, Y, and Z coordinates corresponding to ; a light source with spectral distribution given by the ; function SPEC_INTENS, which is called with a series of ; wavelengths between 380 and 780 nm (the argument is ; expressed in meters), which returns emittance at that ; wavelength in arbitrary units. The chromaticity ; coordinates of the spectrum are returned in the x, y, and z ; arguments which respect the identity: ; x + y + z = 1. function spectrum_to_xyz, spec_intens @./TOOLS/common_color.inc XYZ = DBLARR(3) spectrum = REBIN(spec_intens[380:784], 81) XYZ = spectrum ## cie_colour_match return, XYZ / TOTAL(XYZ) end ; BB_SPECTRUM ; ; Calculate, by Planck's radiation law, the emittance of a black body ; of temperature bbTemp at the given wavelength (in metres). ; function blackbody_spectrum, T lambda = FINDGEN(1000) return, (3.74183e-16 * lambda^(-5.0)) / $ (exp(1.4388e-2 / (lambda * T)) - 1.0) end ; ;/* Built-in test program which displays the x, y, and Z and RGB ; values for black body spectra from 1000 to 10000 degrees kelvin. ; When run, this program should produce the following output: ; ; Temperature x y z R G B ; ----------- ------ ------ ------ ----- ----- ----- ; 1000 K 0.6528 0.3444 0.0028 1.000 0.007 0.000 (Approximation) ; 1500 K 0.5857 0.3931 0.0212 1.000 0.126 0.000 (Approximation) ; 2000 K 0.5267 0.4133 0.0600 1.000 0.234 0.010 ; 2500 K 0.4770 0.4137 0.1093 1.000 0.349 0.067 ; 3000 K 0.4369 0.4041 0.1590 1.000 0.454 0.151 ; 3500 K 0.4053 0.3907 0.2040 1.000 0.549 0.254 ; 4000 K 0.3805 0.3768 0.2428 1.000 0.635 0.370 ; 4500 K 0.3608 0.3636 0.2756 1.000 0.710 0.493 ; 5000 K 0.3451 0.3516 0.3032 1.000 0.778 0.620 ; 5500 K 0.3325 0.3411 0.3265 1.000 0.837 0.746 ; 6000 K 0.3221 0.3318 0.3461 1.000 0.890 0.869 ; 6500 K 0.3135 0.3237 0.3628 1.000 0.937 0.988 ; 7000 K 0.3064 0.3166 0.3770 0.907 0.888 1.000 ; 7500 K 0.3004 0.3103 0.3893 0.827 0.839 1.000 ; 8000 K 0.2952 0.3048 0.4000 0.762 0.800 1.000 ; 8500 K 0.2908 0.3000 0.4093 0.711 0.766 1.000 ; 9000 K 0.2869 0.2956 0.4174 0.668 0.738 1.000 ; 9500 K 0.2836 0.2918 0.4246 0.632 0.714 1.000 ; 10000 K 0.2807 0.2884 0.4310 0.602 0.693 1.000 ;*/ ;int main() ;{ ; double t, x, y, z, r, g, b; ; struct colourSystem *cs = &SMPTEsystem; ; ; printf("Temperature x y z R G B\n"); ; printf("----------- ------ ------ ------ ----- ----- -----\n"); ; ; for (t = 1000; t <= 10000; t+= 500) { ; bbTemp = t; ; spectrum_to_xyz(bb_spectrum, &x, &y, &z); ; xyz_to_rgb(cs, x, y, z, &r, &g, &b); ; printf(" %5.0f K %.4f %.4f %.4f ", t, x, y, z); ; if (constrain_rgb(&r, &g, &b)) { ; norm_rgb(&r, &g, &b); ; printf("%.3f %.3f %.3f (Approximation)\n", r, g, b); ; } else { ; norm_rgb(&r, &g, &b); ; printf("%.3f %.3f %.3f\n", r, g, b); ; } ; } ; return 0; ;} ; pro set_spectral_palette, lambda @TOOLS/common_color.inc @common_blocks.inc common colors, r_orig, g_orig, b_orig, r_curr, g_curr, b_curr spectrum = DBLARR(1000) & spectrum[lambda > 380 < 779] = 1.D color = norm_rgb(constrain_rgb(xyz_to_rgb(CIEsystem, spectrum_to_xyz(spectrum)))) ; fade from 0:black to 192:color to 255:white t = [FINDGEN(192)/192., (63. - FINDGEN(64))/63.] a = [FLTARR(192), 1. + FLTARR(64)] r_curr = FIX(255*(color[0] * t + a * (1.-t))) g_curr = FIX(255*(color[1] * t + a * (1.-t))) b_curr = FIX(255*(color[2] * t + a * (1.-t))) TVLCT, r_curr,g_curr,b_curr end