Previous: Acknowledgments Up: A Computational Model of Color Perception and Color Naming Next: Getting the software
This appendix lists some equations for deriving various color spaces from the CIE XYZ standard, and from each other, in addition to some CIE chromaticity-related equations and functions. The Mathematica code is listed verbatim below, I believe it is not too hard to read.
These are color space transforms of various kinds.
(* Transforms of fundamentals etc. Includes both linear and non-linear transforms. The former are represented as 3x3 matrices, the latter algorithmically. *)BeginPackage["transforms`", {"common`", "CIE`", "chromaticity`"}]
(* The symbols appearing below, before the start of the private part, will be exported from the package defined above. The usage information will be displayed by the help command. *)
XyzToRgbTV::usage = "Linear transform from XYZ to RGB TV coordinates" RgbTVToXyz::usage = "Linear transform from RGB TV to XYZ coordinates" XyzToRgbCie::usage = "Linear transform from XYZ to RGB CIE coordinates" RgbCieToXyz::usage = "Linear transform from RGB CIE to XYZ coordinates" RgbCieToLum::usage = "Linear transform from RGB CIE to luminance" ChromRgbCie::usage = "RGB CIE chromaticity coordinates" XyzToRgbFCC::usage = "Linear transform from XYZ to RGB FCC coordinates" RgbFCCToXyz::usage = "Linear transform from RGB FCC to XYZ coordinates" RgbFCCToLum::usage = "Linear transform from RGB FCC to luminance" ChromRgbNTSC::usage = "RGB NTSC chromaticity coordinates" XyzToRgbNTSC::usage = "Linear transform from XYZ to RGB NTSC coordinates" RgbNTSCToXyz::usage = "Linear transform from RGB NTSC to XYZ coordinates" ChromRgbCRT::usage = "RGB CRT chromaticity coordinates" XyzToRgbCRT::usage = "Linear transform from XYZ to RGB CRT coordinates" RgbCRTToXyz::usage = "Linear transform from RGB CRT to XYZ coordinates" RgbToI123::usage = "Linear transform from generic RGB to Ohta's I123 coordinates" I123ToRgb::usage = "Linear transform from Ohta's I123 to generic RGB coordinates" XyzToI123::usage = "Linear transform from XYZ to I123 coordinates, using RGB CRT as the intermediate norm" I123ToXyz::usage = "Linear transform from I123 to XYZ coordinates, using RGB CRT as the intermediate norm" YiqToRgb::usage = "Linear transform from YIQ to generic RGB coordinates" RgbToYiq::usage = "Linear transform from generic RGB to YIQ coordinates" XyzToYiq::usage = "Linear transform from XYZ to YIQ coordinates, using RGB NTSC as the intermediate norm" YiqToXyz::usage = "Linear transform from YIQ to XYZ coordinates, using RGB NTSC as the intermediate norm" ChromRgbSony::usage = "Typical Sony monitor RGB chromaticities" ChromWSony::usage = "Typical Sony monitor white chromaticities" TriStimWSony::usage = "Typical Sony monitor while tristimulus values" XyzToRgbSony::usage = "Linear transform from XYZ to RGB Sony coordinates" RgbSonyToXyz::usage = "Linear transform from RGB Sony to XYZ coordinates" RgbCRTToRgbSony::usage = "Linear transform from RGB CRT to RGB Sony coordinates" XyzToLmsVW::usage = "Linear transform from XYZ to LMS Vos & Walraven coordinates" LmsVWToXyz::usage = "Linear transform from LMS Vos & Walraven to XYZ coordinates" XyzToLmsE::usage = "Linear transform from XYZ to LMS Estevez coordinates" LmsEToXyz::usage = "Linear transform from LMS Estevez to XYZ coordinates" XyzToLmsSP::usage = "Linear transform from XYZ to LMS Smith & Pokorny coordinates" LmsSPToXyz::usage = "Linear transform from LMS Smith & Pokorny to XYZ coordinates" XyzToH81::usage = "Linear transform from XYZ to Hurvich & Jameson opponent coordinates" H81ToXyz::usage = "Linear transform from Hurvich & Jameson opponent to XYZ coordinates" RgbCRTToH81::usage = "Linear transform from RGB CRT to Hurvich & Jameson opponent coordinates" H81ToRgbCRT::usage = "Linear transform from Hurvich & Jameson opponent to RGB CRT coordinates" LmsVWToH81::usage = "Linear transform from LMS Vos & Walraven to Hurvich & Jameson opponent coordinates" H81ToLmsVW::usage = "Linear transform from Hurvich & Jameson to LMS Vos & Walraven coordinates" LmsVWToWWAW::usage = "Linear transform from LMS Vos & Walraven to Werner & Wooten subject AW opponent coordinates" WWAWToLmsVW::usage = "Linear transform from Werner & Wooten subject AW opponent to LMS Vos & Walraven coordinates" RgbCRTToWWAW::usage = "Linear transform from RGB CRT to Werner & Wooten subject AW opponent coordinates" WWAWToRgbCRT::usage = "Linear transform from Werner & Wooten subject AW opponent to RGB CRT coordinates" XyzToWWAW::usage = "Linear transform from XYZ to Werner & Wooten subject AW opponent coordinates" WWAWToXyz::usage = "Linear transform from Werner & Wooten subject AW opponent to XYZ coordinates" RgbToHsi::usage = "RgbToHsi[{r_,g_,b_}] returns {h,s,i} in {[0,2 Pi], [0,1], [0,1]} corresponding to {r_,g_,b_} in [0,1]." XyzToUvl::usage = "XyzToUvl[{{x_,y_,z_}, {xI_,yI_,zI_}}] returns {u,v,l} in {R,R,[0,100]}, corresponding to {x_,y_,z_} in R+ and {xI_,yI_,zI_} in R+ (the latter representing illuminant color). This is the CIE 1976 L*u*v* space (CIELUV)." diffUvl::usage = "diffUvl[{u1_,v1_,l1_}, {u2_,v2_,l2_}] returns the CIELUV color difference between the colors given by {u1_,v1_,l1_} and {u2_,v2_,l2_}, both in {R,R,[0,100]}. This amounts to Eucidean distance in CIELUV space." UvlToHcl::usage = "UvlToHcl[{u_,v_,l_}] returns {hue, chroma, lightness} in {[0,2Pi], R+, [0,100]} corresponding to {u_,v_,l_} in {R,R,[0,100]}, representing a color in CIELUV coordinates. Chroma changes with changing lightness and constant chromaticity." UvlToHsl::usage = "UvlToHsl[{u_,v_,l_}] returns {hue, saturation, lightness} in {[0,2Pi], R+, [0,100]} corresponding to {u_,v_,l_} in {R,R,[0,100]}, representing a color in CIELUV coordinates. Saturation does not change with changing lightness and constant chromaticity." XyzToAbl::usage = "XyzToAbl[{{x_,y_,z_}, {xI_,yI_,zI_}}] returns {a,b,l} in {R,R,[0,100]}, corresponding to {x_,y_,z_} in R+ and {xI_,yI_,zI_} in R+ (the latter representing illuminant color). This is the CIE 1976 L*a*b* space (CIELAB)." diffAbl::usage = "diffAbl[{a1_,b1_,l1_}, {a2_,b2_,l2_}] returns the CIELAB color difference between the colors given by {a1_,b1_,l1_} and {a2_,b2_,l2_}, both in {R,R,[0,100]}. This amounts to Eucidean distance in CIELAB space." AblToHcl::usage = "AblToHcl[{a_,b_,l_}] returns {hue, chroma, lightness} in {[0,2Pi], R+, [0,100]} corresponding to {a_,b_,l_} in {R,R,[0,100]}, representing a color in CIELAB coordinates. Chroma changes with changing lightness and constant chromaticity." AblToHsl::usage = "AblToHsl[{a_,b_,l_}] returns {hue, saturation, lightness} in {[0,2Pi], R+, [0,100]} corresponding to {a_,b_,l_} in {R,R,[0,100]}, representing a color in CIELAB coordinates. Saturation does not change with changing lightness and constant chromaticity." XyzToSVF::usage = "XyzToSVF[{{X_,Y_,Z_}, {XW_,YW_,ZW_}}] returns {F1,F2,VY} in {R,R,R+} corresponding to {X_,Y_,Z_} in R+ and {XW_,YW_,ZW_} in R+. The latter represents the XYZ values of white. {F1,F2,VY} are the opponent coordinates F1 and F2 and the lightness magnitude VY of the SVF uniform color space." Rgb2Hls::usage = "Rgb2Hls[{r_,g_,b_}] returns {h,l,s} in {[0,2Pi], [0,1], [0,1]} corresponding to {r_,g_,b_} in [0,1]." Hls2Rgb::usage = "Hls2Rgb[{h_,l_,s_}] returns {r,g,b} in [0,1] corresponding to {h_,l_,s_} in {[0,2Pi], [0,1], [0,1]}." Rgb2Hsv::usage = "Rgb2Hsv[{r_,g_,b_}] returns {h,s,v} in {[0,2Pi], [0,1], [0,1]} corresponding to {r_,g_,b_} in [0,1]." Hsv2Rgb::usage = "Hsv2Rgb[{h_,s_,v_}] returns {r,g,b} in [0,1] corresponding to {h_,s_,v_} in {[0,2Pi], [0,1], [0,1]}."
Begin["`Private`"]
(* The symbols appearing below are private to this package, and will not be exported. *)
(* LINEAR TRANSFORMS *)
(* Ency of AI *)
(* Note: EAI 1992 edition CV article says linear transforms should be scaled s.t. each component has a max of 1; I don't do this. May be better for statistical feature counting, but not for my purpose. This will distort the relative positions of colors in the space. *)
XyzToRgbTV := {r0TV, g0TV, b0TV} * {{0.587, -0.164, -0.089}, {-0.301, 0.611, -0.0087}, {0.0178, -0.0362, 0.274}}
RgbTVToXyz := Inverse[XyzToRgbTV]
(* Compute scaling factors to preserve E response relative to Cie XYZ functions. This does not preserve absolute magnitudes of transformed primaries, so they are renormalized afterwards. The transforms are defined with delayed evaluation, so that the scaling factors can be redefined at any time and the transforms will be modified accordingly. Note that TV fundamentals are equalized on standard illuminant C, not on equal-energy white as the XYZ functions are. *)
{r0TV, g0TV, b0TV} = {1,1,1} {r0TV, g0TV, b0TV} = 1 / (XyzToRgbTV . TriStimC) {r0TV, g0TV, b0TV} = RgbNormFactorFromTxr[XyzToRgbTV, Cie31List] {r0TV, g0TV, b0TV}
(* relation Cie 1931 XYZ to Cie real primaries 700, 546.1, 435.8 nm, see McIlwain & Dean 1956 Eq 4-4 p. 48 *)
RgbCieToXyz1 = {{2.7690, 1.7518, 1.1300}, {1.0000, 4.5907, 0.0601}, {0.0000, 0.0565, 5.5943}}
XyzToRgbCie1 = Inverse[RgbCieToXyz1]
XyzToRgbCie := {r0Cie, g0Cie, b0Cie} * XyzToRgbCie1
RgbCieToXyz := Inverse[XyzToRgbCie]
{r0Cie, g0Cie, b0Cie} = {1,1,1} {r0Cie, g0Cie, b0Cie} = 1 / (XyzToRgbCie . TriStimE) (* note eq on E *) {r0Cie, g0Cie, b0Cie} = RgbNormFactorFromTxr[XyzToRgbCie, Cie31List] {r0Cie, g0Cie, b0Cie}
RgbCieToLum = {1, 4.5907, 0.0601} (* luminance contributions *)
ChromRgbCie = {{0.735, 0.265}, (* from rogers85 *) {0.274, 0.717}, {0.167, 0.009}}
(* I have used ChromRgbCie to compute matching functions from the chromaticities, with TxrFromChrom[ChromRgbCie, ChromE], which gives the same results as the transformation above. That means that, most likely, TxrFromChrom is working ok, I can trust it to compute matching functions based on chromaticities. *)
(* McIlwain & Dean conversions between Cie 1931 XYZ and FCC RGB used as the primaries for color tv. These are equalized on standard illuminant C as reference white, with luminance of this white taken as unity. Eqns 4-11 p. 62.
After scaling, RgbFCC is vitually identical to RgbTV. *)
RgbFCCToXyz1 = {{0.608, 0.174, 0.200}, {0.299, 0.587, 0.114}, {0.000, 0.0662, 1.112}}
(* compute inverse transform rather than using eq 4-12, which is apparently wrong: 1.191X should read 1.907X, perhaps 1.911X *)
XyzToRgbFCC1 = Inverse[RgbFCCToXyz1]
XyzToRgbFCC := {r0FCC, g0FCC, b0FCC} * XyzToRgbFCC1
RgbFCCToXyz := Inverse[XyzToRgbFCC]
{r0FCC, g0FCC, b0FCC} = {1,1,1} {r0FCC, g0FCC, b0FCC} = 1 / (XyzToRgbFCC . TriStimC) {r0FCC, g0FCC, b0FCC} = RgbNormFactorFromTxr[XyzToRgbFCC, Cie31List] {r0FCC, g0FCC, b0FCC}
(* luminance signal as used in color tv, eq 8-2 p. 121, without gamma correction *)
RgbFCCToLum = {0.30, 0.59, 0.11}
(* Sun VideoPix doc: 0.299, 0.587, 0.114 for NTSC *)
(* NTSC from Xyz chromaticity coordinates, from Hill 1990, note 3 p. 574 After scaling, RgbNTSC is vitually identical to RgbTV. *)
ChromRgbNTSC = {{0.670, 0.330}, {0.210, 0.710}, {0.140, 0.080}}
XyzToRgbNTSC = TxrFromChrom[ChromRgbNTSC, ChromC]
RgbNTSCToXyz = Inverse[XyzToRgbNTSC]
(* typical CRT color monitor from chromaticity coordinates, from Hill 1990, table 16.2 p. 574, and rest of the chapter. Note this is equalized on CIE D65, not on C as the NTSC functions are. These are slightly different from the values given in rogers85. *)
ChromRgbCRT = {{0.628, 0.330}, (* rogers85: 0.628 0.346 *) {0.285, 0.590}, (* 0.268 0.588 *) {0.1507, 0.060}} (* 0.150 0.070 *)
(* these are fairly close to chromaticities computed on positive parts only of RgbFCC matching curves, using this kind of computation: {Function[l, If[R[l] < 0, 0, R[l]]], Function[l, If[G[l] < 0, 0, G[l]]], Function[l, If[B[l] < 0, 0, B[l]]]} *)
(* used a transform computed from chromaticity coordinates rather than the one given in Hill. The differences are considerable. *)
XyzToRgbCRT = TxrFromChrom[ChromRgbCRT, ChromD65] (* Hill has the following: {r0CRT, g0CRT, b0CRT} * {{2.739, -1.119, 0.138}, {-1.145, 2.2029, -0.333}, {-0.424, 0.033, 1.105}} *)
RgbCRTToXyz = Inverse[XyzToRgbCRT]
(* equalize on source D65 for white *)
{r0CRT, g0CRT, b0CRT} = {1,1,1} {r0CRT, g0CRT, b0CRT} = 1 / (XyzToRgbCRT . TriStimD65) {r0CRT, g0CRT, b0CRT} = RgbNormFactorFromTxr[XyzToRgbCRT, Cie31List] {r0CRT, g0CRT, b0CRT}
(* EAI Ohta's I1,I2,I3 space. This is a linear transform of RGB. Note the purpose of this trafo, described in the article. *)
RgbToI123 = {{1/3, 1/3, 1/3}, {1, 0, -1}, {-1/2, 1, -1/2}}
I123ToRgb = Inverse[RgbToI123]
XyzToI123 = RgbToI123 . XyzToRgbCRT
XyzToI123 = NormFactorFromTxr[XyzToI123, 550, 610, 450, Cie31List] XyzToI123
I123ToXyz = Inverse[XyzToI123]
(* RGB<-->YIQ transforms from mma Packages/Graphics/Colors.m *)
YiqToRgb = {{1, .95, .625}, {1, -.28, -.64}, {1, -1.11, 1.73}} (* rogers85 has a typo omitting the 1 from the last number: {{1, .956, .623}, {1, -.272, -.648}, {1, -1.105, .705}} *)
RgbToYiq = Inverse[YiqToRgb] (* rogers85: {{.299, .587, .114}, {.596, -.274, -.322}, {.211, -.522, .311}} *)
XyzToYiq = RgbToYiq . XyzToRgbNTSC
XyzToYiq = NormFactorFromTxr[XyzToYiq, 550, 610, 540, Cie31List] XyzToYiq
YiqToXyz = Inverse[XyzToYiq]
(* Sony CRT color monitor from Xyz chromaticity coordinates, from twj@wri.com. These are the values used in Mathematics as "approximations to typical color monitors", i.e. the Sony tubes as used in Mac's etc. Typical gamma for these tubes is 1.8. Agrees with values gotten from Sun. *)
ChromRgbSony = {{0.625, 0.340}, {0.280, 0.595}, {0.155, 0.070}}
ChromWSony = {0.283, 0.298} TriStimWSony = {0.283, 0.298, 0.419}
XyzToRgbSony = TxrFromChrom[ChromRgbSony, ChromWSony]
RgbSonyToXyz = Inverse[XyzToRgbSony]
(* use the following transform to display RgbCRT colors on a Sony monitor *)
RgbCRTToRgbSony = XyzToRgbSony . RgbCRTToXyz
(* Vos & Walraven 1978 fundamentals. W/o scaling factors, the resulting functions are on luminance basis. From Wyszecki & Stiles 2nd ed, p. 612 ff. The transformation is relative to the Vos-Judd 1978 modified CIE standard functions (Cie31VxBar, Cie31VyBar, Cie31VzBar). *)
XyzToLmsVW := {r0VW, g0VW, b0VW} * {{0.1551646, 0.5430763, -0.0370161}, {-0.1551646, 0.4569237, 0.0296946}, {0.0, 0.0, 0.0073215}}
LmsVWToXyz := Inverse[XyzToLmsVW]
(* equalized on E *)
{r0VW, g0VW, b0VW} = {1,1,1}; {r0VW, g0VW, b0VW} = 1 / (XyzToLmsVW . TriStimE) {r0VW, g0VW, b0VW} = RgbNormFactorFromTxr[XyzToLmsVW, Cie31VList] {r0VW, g0VW, b0VW}
(* Estevez 1979 fundamentals, according to Valberg et al 1986. The transformation is relative to the CIE 1931 standard functions (Cie31VxBar, Cie31VyBar, Cie31VzBar). *)
XyzToLmsE := {r0E, g0E, b0E} * {{0.3841, 0.7391, -0.0650}, {-0.3471, 1.1463, 0.0870}, {0.0, 0.0, 0.5610}}
LmsEToXyz := Inverse[XyzToLmsE]
(* equalized on E *)
{r0E, g0E, b0E} = {1,1,1}; {r0E, g0E, b0E} = 1 / (XyzToLmsE . TriStimE) {r0E, g0E, b0E} = RgbNormFactorFromTxr[XyzToLmsE, Cie31VList] {r0E, g0E, b0E}
(* Smith & Pokorny 1975 fundamentals. W/o scaling factors, the resulting functions are on luminance basis. From Wyszecki & Stiles 2nd ed, p. 612 ff. The transformation is relative to the Judd 1951 modified CIE standard functions (Cie31JxBar, Cie31JyBar, Cie31JzBar). *)
XyzToLmsSP := {r0SP, g0SP, b0SP} * {{0.15514, 0.54312, -0.03286}, {-0.15514, 0.45684, 0.03286}, {0, 0, 0.00801}}
LmsSPToXyz := Inverse[XyzToLmsSP]
(* equalized on E *)
{r0SP, g0SP, b0SP} = {1,1,1}; {r0SP, g0SP, b0SP} = 1 / (XyzToLmsSP . TriStimE) {r0SP, g0SP, b0SP} = RgbNormFactorFromTxr[XyzToLmsSP, Cie31JList] {r0SP, g0SP, b0SP}
(* Hurvich 81 opponent responses *)
XyzToH81 := {gr0H81, by0H81, wb0H81} {{1, -1, 0}, {0, 0.4, -0.4}, {0, 1, 0}}
H81ToXyz := Inverse[XyzToH81]
RgbCRTToH81 := XyzToH81 . RgbCRTToXyz
H81ToRgbCRT := XyzToRgbCRT . H81ToXyz
LmsVWToH81 := XyzToH81 . LmsVWToXyz
H81ToLmsVW := XyzToLmsVW . H81ToXyz
{gr0H81, by0H81, wb0H81} = {1,1,1} {gr0H81, by0H81, wb0H81} = OppNormVectorFromTxr[XyzToH81] {gr0H81, by0H81, wb0H81}
(* W&W 79, w.r.t. normalized V&W fundamentals, observer AW *)
LmsVWToWWAWOrig = {{2.38, -2.87, 0.54}, {0.93, -0.36, -1.03}, {0.85, 0.15, 0.01}}
(* this is the normalizing matrix for VW fundamentals. The argument to DiagonalMatrix is computed as 1 / {rBarVWMax, gBarVWMax, bBarVWMax} *)
normVW = DiagonalMatrix[{1.03759, 0.844756, 0.618546}]
(* derive WWAW fns from unnormalized VW fundamentals and from rgb *)
LmsVWToWWAW := LmsVWToWWAWOrig . normVW
WWAWToLmsVW := Inverse[LmsVWToWWAW]
RgbCRTToWWAW := LmsVWToWWAW . XyzToLmsVW . RgbCRTToXyz
WWAWToRgbCRT := Inverse[RgbCRTToWWAW]
XyzToWWAW := {gr0WWAW, by0WWAW, wb0WWAW} (LmsVWToWWAW . XyzToLmsVW)
WWAWToXyz := Inverse[XyzToWWAW]
{gr0WWAW, by0WWAW, wb0WWAW} = {1,1,1} {gr0WWAW, by0WWAW, wb0WWAW} = OppNormVectorFromTxr[XyzToWWAW] {gr0WWAW, by0WWAW, wb0WWAW}
(* NON-LINEAR TRANSFORMS *)
(* EAI: HSI space used in color vision work. Note this is not a linear transform! Note also the singularities in this fn at low intensities. Modified the transform for for i below .03, otherwise get strange results with black, e.g. At low intensities, H and S are meaningless, almost random, using the unmodified transform. Assume r,g,b, in [0,1]. This version is better for black, but still unreliable at low intensities, e.g. for brown. *)
RgbToHsi[{r_,g_,b_}] := Module[ {h, s, i, rn, gn, bn}, ( i = (r+g+b)/3; If[i<.02, {0,0,i}, ( {rn,gn,bn} = {r,g,b}/i; s = 1 - Min[{rn,gn,bn}]; x = ArcCos[(2 rn-gn-bn) / (Sqrt[6] Sqrt[(rn-1/3)^2 + (gn-1/3)^2 + (bn-1/3)^2])]; h = If[bn<=gn, x, 2 Pi-x]; {h,s,i} )] )]
(* EAI: L*u*v* perceptually uniform color space, superseding the earliest uniform space UVW. Quantities denoted by postfix I in the definition of L*u*v* refer to the incident illumination color (!). Note this is not a linear transform of XYZ space. Note the singularities for low intensity here too. This incorporates the modification for low values of y/yI as described in W&S p. 165. This is CieLUV 1976. Note these function actually use the order U, v, l to be compatible with H81 etc, which have brightness as the third dimension (vertical). *)
XyzToUvl[{{x_,y_,z_}, {xI_,yI_,zI_}}] := Module[ {ls, us, vs, u, v, uI, vI, d, dI}, ( d = x + 15 y + 3 z; dI = xI + 15 yI + 3 zI; u = If[d>0, 4 x / d, 0]; uI = If[dI>0, 4 xI / dI, 0]; v = If[d>0, 9 y / d, 0]; vI = If[dI>0, 9 yI / dI, 0]; ls = If[(y/yI)>0.008856, 116 (y/yI)^(1/3) - 16, 903.3 (y/yI)]; us = 13 ls (u-uI); vs = 13 ls (v-vI); {us,vs,ls} )]
(* color difference in CieLUV 1976 is just Euclidean distance *)
diffUvl[{u1_,v1_,l1_}, {u2_,v2_,l2_}] := Sqrt[(u1-u2)^2 + (v1-v2)^2 + (l1-l2)^2]
(* Derive perceived lightness, chroma, and hue from uvl coordinates. Chroma is similar to saturation, except the latter (as defined in W&S) does not change with changing lightness and constant chromaticity, but the former does. *)
UvlToHcl[{u_,v_,l_}] := If[u==0 && v==0, {0,0,l}, {ArcTan[u,v], Sqrt[u^2 + v^2], l}] (* chroma *)
UvlToHsl[{u_,v_,l_}] := If[u==0 && v==0, {0,0,l}, If[l==0, {0,0,l}, {ArcTan[u,v], Sqrt[u^2 + v^2] / l, l}]] (* sat *)
(* Cie 1976 L*a*b* space and color-difference formula, from W&S p. 166. Note again the order a, b, L of the components. *)
XyzToAbl[{{x_,y_,z_}, {xI_,yI_,zI_}}] := Module[ {ls,as,bs}, ( ls = If[(y/yI)>0.008856, 116 (y/yI)^(1/3) - 16, 903.3 (y/yI)]; as = 500 (If[(x/xI)>0.008856, (x/xI)^(1/3), 7.787 (x/xI) + 16/116] - If[(y/yI)>0.008856, (y/yI)^(1/3), 7.787 (y/yI) + 16/116]); bs = 200 (If[(y/yI)>0.008856, (y/yI)^(1/3), 7.787 (y/yI) + 16/116] - If[(z/zI)>0.008856, (z/zI)^(1/3), 7.787 (z/zI) + 16/116]); {as,bs,ls} )]
(* color difference in CieLAB 1976 is just Euclidean distance *)
diffAbl[{a1_,b1_,l1_}, {a2_,b2_,l2_}] := Sqrt[(a1-a2)^2 + (b1-b2)^2 + (l1-l2)^2]
(* Derive perceived lightness, chroma, and hue from Abl coordinates. *)
AblToHcl[{a_,b_,l_}] := If[a==0 && b==0, {0,0,l}, {ArcTan[a,b], Sqrt[a^2 + b^2], l}] (* chrom *)
(* same but with saturation in stead of chroma -- not in W&S? analogous to Uvl to Hcl, Hsl transforms above. *)
AblToHsl[{a_,b_,l_}] := If[a==0 && b==0, {0,0,l}, If[l==0, {0,0,l}, {ArcTan[a,b], Sqrt[a^2 + b^2] / l, l}]] (* sat *)
(* SVF uniform color space, from Valberg et al 1986, app. A *)
(* helper funcions *) v1[Y_] := If[Y<=0.0043, 0, ((100 Y - 0.43)^0.51)/((100 Y - 0.43)^0.51 + 31.75)]
k[VY_] := 0.140 + 0.175 VY
v2[Y_] := Module[ {VY}, VY = 40 v1[Y]; If[Y<=0.001 k[VY], 0, (((100 Y / k[VY]) - 0.1)^0.86) / (((100 Y / k[VY]) - 0.1)^0.86 + 103.2) ] ]
XyzToSVF[{{X_,Y_,Z_}, {XWhite_,YWhite_,ZWhite_}}] := Module[ {S1,S2,S3,VY,p1,p2,F1,F2}, (* center relative cone absorptions of color sample to the white stimulus (or the light source), aka von Kries transformation. Trafo A1 given in Valberg et al 86 p. 1732, differs by scaling factors only from XyzToLmsE. Since the sample values are scaled w.r.t. the white values, this doesn't matter. *) {S1,S2,S3} = (XyzToLmsE . {X,Y,Z}) / (XyzToLmsE . {XWhite,YWhite,ZWhite}); (* lighness magnitude VY *) VY = 40 v1[Y]; (* first opponent stage coordinates p1, p2 *) p1 = v1[S1] - v1[Y]; p2 = If[S3<=Y, v1[Y] - v1[S3], v2[Y] - v2[S3]]; (* second opponent stage coordinates F1 and F2 *) F1 = 700 p1 - 54 p2; F2 = 96.5 p2; {F1,F2,VY} ]
(* rgb to hls triplets, from Hill 90. In this hls space, all pure colors (red, green, blue, yellow, cyan, magenta) lie on the 0.5 lightness plane and are all equally saturated (1.0), with black and white at the bottom and top of the double cone. This is not very realistic in psychophysical terms. For greys, h=0 although strictly speaking it is undefined. Otherwise, h is given in radians, with 0 = red.
Foley & Van Dam 1982(84) have essentially the same routine, except they differ on the sign of one of the terms! (see code). *)
Rgb2Hls[{r_,g_,b_}] := Module[ {mx, mn, rc, gc, bc, h, l, s}, ( (* convert {r,g,b}, each in [0,1], to {h,l,s}, in {[0,2Pi], [0,1], [0,1]} *) mx = Max[r,g,b]; mn = Min[r,g,b]; l = (mx + mn) / 2.0; (* lightness *) (* compute saturation *) If[mx == mn, s = h = 0, ( (* grey *) (* chromatic color *) If[l<=0.5, s = (mx-mn) / (mx+mn), s = (mx-mn) / (2-mx-mn)]; (* Hill has (2-mx+mn) for the last term, which is wrong as evidenced by computing some transforms to HLS and back to RGB. *) rc = (mx-r) / (mx-mn); (* hue *) gc = (mx-g) / (mx-mn); bc = (mx-b) / (mx-mn); If[r==mx, h = bc-gc, If[g==mx, h = 2+rc-bc, If[b==mx, h = 4+gc-rc]]]; h = 60 h; If[h<0, h = h+360] )]; {h (Pi/180), l, s} )]
(* this is from Foley & vDam 82, p 619 *)
Hls2Rgb[{h_,l_,s_}] := Module[ {value, h1, r, g, b, m1, m2}, ( (* convert {h,l,s} in {[0,2Pi], [0,1], [0,1]} to {r,g,b}, each in [0,1] *) value[n1_,n2_,hue_] := Module[ {hue1}, ( If[hue>360, hue1=hue-360, If[hue<0, hue1=hue+360, hue1=hue]]; Which[ hue1<60, n1+(n2-n1) hue1/60, hue1<180, n2, hue1<240, n1+(n2-n1) (240-hue1)/60, True, n1])]; h1 = h (180/Pi); (* convert to degrees *) If[l<=0.5, m2=l(1+s), m2=l+s-l s]; m1=2 l-m2; If[s==0, {r,g,b} = {l,l,l}, {r,g,b} = {value[m1,m2,h1+120], value[m1,m2,h1], value[m1,m2,h1-120]}]; {r,g,b} )]
(* HSV hexcone after Smith 78, from F&vD 82, 613ff. As in HLS above, saturation is relative to gamut represented by model, not to Cie chart, i.e. not the same as purity. HSV = Hue, Saturation, Value (lightness, after Munsell value I suppose). For the following conversion functions, {r,g,b} are in [0,1], {h, s, v} in {[0, 2Pi], [0,1], [0,1]}. For undefined values of h, 0 is returned. *)
Rgb2Hsv[{r_,g_,b_}] := Module[ {mx, mn, rc, gc, bc, h, s, v}, ( (* convert {r,g,b}, each in [0,1], to {h,s,v}, in {[0,2Pi], [0,1], [0,1]} *) mx = Max[r,g,b]; mn = Min[r,g,b]; v = mx; (* value *) If[mx != 0, s = (mx-mn) / mx, s = 0]; (* saturation *) If[s==0, h = 0, ( (* grey *) (* chromatic color *) rc = (mx-r) / (mx-mn); (* distance from red *) gc = (mx-g) / (mx-mn); bc = (mx-b) / (mx-mn); If[r==mx, h = bc-gc, (* between Y, Magenta *) If[g==mx, h = 2+rc-bc, (* between Cyan, Y *) If[b==mx, h = 4+gc-rc]]]; (* between Magenta, Cyan *) h = 60 h; (* convert to degrees *) If[h<0, h = h+360] (* make nonnegative *) )]; {h (Pi/180), s, v} (* convert to radians *) )]
Hsv2Rgb[{h_,s_,v_}] := Module[ {r, g, b, h1, i, f, p, q, t}, ( (* convert {h,s,v}, in {[0,2Pi], [0,1], [0,1]} to {r,g,b}, each in [0,1] *) h1 = h (180/Pi); (* convert to degrees *) If[s==0, {r,g,b} = {v,v,v}, ( If[h1==360, h1=0]; (* must use local variable for h *) h1 = h1/60; (* convert to [0,6] *) i = Floor[h1]; (* integer part of h *) f = h1-i; (* fractional part of h *) p = v (1-s); q = v (1-(s f)); t = v (1-(s (1-f))); Switch[i, 0, {r,g,b} = {v,t,p}, 1, {r,g,b} = {q,v,p}, 2, {r,g,b} = {p,v,t}, 3, {r,g,b} = {p,q,v}, 4, {r,g,b} = {t,p,v}, 5, {r,g,b} = {v,p,q}] )]; {r,g,b} )]
End[] (* private *) EndPackage[] (* package *)
These are CIE chromaticity-related equations and functions.
(* Cie chromaticity functions and related stuff *)BeginPackage["chromaticity`", {"common`", "CIE`", "mygraphics`"}]
(* The symbols appearing below, before the start of the private part, will be exported from the package defined above. The usage information will be displayed by the help command. *)
SetDefaultCieFns::usage = "SetDefaultCieFns[CieFnsList] sets the default CIE basis functions (x, y, z) to the list CieFnsList."
ChromDiagram::usage = "Parametric plot of the CIE chromaticity diagram"
ChromDiagramColor::usage = "Color parametric plot of the CIE chromaticity diagram"
PlotChromList::usage = "PlotChromList[{{x,y}, ...}] plots the CIE chromaticity diagram and points with chromaticities {x,y}."
PlotChromListColor::usage = "PlotChromList[{{x,y}, ...}] plots the CIE chromaticity diagram and points with chromaticities {x,y}, in color."
TxrFromChrom::usage = "TxrFromChrom[{{Rx, Ry}, {Gx, Gy}, {Bx, By}}, {Wx,Wy}] computes an XYZ to RGB transform from RGB chromaticity coordinates {{Rx, Ry}, ...}, equalized on the white point given by chromaticity coordinates {Wx, Wy}. The result is normalized for max[rgb]=1."
ChromFromTxr::usage = "ChromFromTxr[txr] returns chromaticity coordinates {x,y} of primaries, given a linear transform txr from Cie XYZ color matching functions to the color matching functions of those primaries."
ChromFromSPD::usage = "ChromFromSPD[{f1, f2, ...}] returns chromaticity coordinates {x, y} of SPD's f1, f2, ... (functions of wavelength)."
LBlue::usage = "Typical wavelength resulting in perception of blue." LGreen::usage = "Typical wavelength resulting in perception of green." LYellow::usage = "Typical wavelength resulting in perception of yellow." LRed::usage = "Typical wavelength resulting in perception of red."
PlotTxr::usage = "PlotTxr[matrix, label:\"\", style:rgbColors, CieFns:Cie31List] plots the primaries defined by matrix, a linear transform from XYZ (with basis functions CieFns) to the primaries, and labels the plot with label. Uses style to plot the individual functions."
PlotTxrP::usage = "PlotTxrP[matrix, label:\"\", style:rgbColors, CieFns:Cie31List] plots the primaries defined by matrix, a linear transform from XYZ (with basis functions CieFns) to the primaries, and labels the plot with label. Plots only the positive lobes of the primaries (useful for RGB primaries only). Uses style to plot the individual functions."
PlotTxrPM::usage = "PlotTxrPM[matrix, label:\"\", style:rgbColors, CieFns:Cie31List] plots the primaries defined by matrix, a linear transform from XYZ (with basis functions CieFns) to the primaries, and labels the plot with label. Plots only the major positive lobes of the primaries (useful for RGB primaries only). Uses style to plot the individual functions."
NIntTxr::usage = "NIntTxr[matrix] computes the integrals over the visible wavelength range of the RGB primaries defined by matrix, a linear transform from XYZ coordinates."
NIntTxrP::usage = "NIntTxr[matrix] computes the integrals over the visible wavelength range of the RGB primaries defined by matrix, a linear transform from XYZ coordinates. Uses only the positive lobes of the primaries."
NIntTxrPM::usage = "NIntTxr[matrix] computes the integrals over the visible wavelength range of the RGB primaries defined by matrix, a linear transform from XYZ coordinates. Uses only the major positive lobes of the primaries."
Txr2ChromGamut::usage = "Txr2ChromGamut[txr,range,points,label] plots the gamut of the primaries defined by txr (a linear transform from XYZ to the primaries) in a square region of the chromaticity diagram with sides (0,range), with a vertical and horizontal resolution of points. The plot is labeled with label. Actually displayed colors are approximations for CRTs. Good values for range and points are 0.85 and 64."
GamutNTSC::usage = "NTSC gamut as computed by function Txr2ChromGamut" GamutCRT::usage = "CRT gamut as computed by function Txr2ChromGamut" GamutLMS::usage = "LMS gamut as computed by function Txr2ChromGamut"
Tor2Gamut::usage = "Tor2Gamut[tor, brightness, hrange, vrange, points, label] computes the gamut of the opponent primaries defined by tor (a linear transform from opponent functions to rgb) at the specified brightness level. The plot is displayed in 3D, with x and y (opponent) coordinates going from -hrange to hrange and z (brightness) coordinates going from 0 to vrange. The x and y resolution is given by points, and the plot is labeled label. Varying the brightness level while keeping vrange constant can be used to generate successive equal brightness planes through the color space for animation. The acually displayed colors are approximations for typical CRTs."
Tor2GamutN::usage = "Like Tor2Gamut, but with RGB intensities normalized for max(r,g,b)=1. This looses all intensity information, but can be used to judge hues better than with Tor2Gamut."
AnimH81CRT::usage = "Animated plot of constant brightness planes through the Hurvich 81 color opponent space, transformed to CRT RGB coordinates."
AnimH81Lms::usage = "Animated plot of constant brightness planes through the Hurvich 81 color opponent space, transformed to LMS coordinates."
AnimWWAWCRT::usage = "Animated plot of constant brightness planes through the Werner&Wooten AW color opponent space, transformed to CRT RGB coordinates."
Begin["`Private`"]
(* The symbols appearing below are private to this package, and will not be exported. *)
(* default basis set to use for CIE XYZ space *)
XyzList := Cie31List
SetDefaultCieFns[CieFns_:Cie31List] := XyzList = CieFns
(* chromaticity coordinates of monochromatic light *)
z[l_] := XyzList[[3]][l] / (XyzList[[1]][l] + XyzList[[2]][l] + XyzList[[3]][l]) y[l_] := XyzList[[2]][l] / (XyzList[[1]][l] + XyzList[[2]][l] + XyzList[[3]][l]) x[l_] := XyzList[[1]][l] / (XyzList[[1]][l] + XyzList[[2]][l] + XyzList[[3]][l])
chrom[l_] := {x[l], y[l]} chrom3[l_] := {x[l], y[l], 1 - x[l] - y[l]}
cPoint[{x_, y_}] := Graphics[Circle[{x,y}, 0.008]] cPointL[l_] := cPoint[chrom[l]]
ChromDiagram := Show[ParametricPlot[{x[l], y[l]}, {l,380,750}, Frame->True, AspectRatio->Automatic, GridLines->Automatic, DisplayFunction->Identity], Graphics[{Line[{{x[380], y[380]}, {x[750], y[750]}}]}, DisplayFunction->Identity]]
PlotChromList[list_] := Show[ChromDiagram, Map[cPoint, list], DisplayFunction->$DisplayFunction]
(* generic rgb or xyz functions based on transforms *)
X[l_] := XyzList[[1]][l] Y[l_] := XyzList[[2]][l] Z[l_] := XyzList[[3]][l]
R[l_] := txr[[1]] . {XyzList[[1]][l], XyzList[[2]][l], XyzList[[3]][l]} G[l_] := txr[[2]] . {XyzList[[1]][l], XyzList[[2]][l], XyzList[[3]][l]} B[l_] := txr[[3]] . {XyzList[[1]][l], XyzList[[2]][l], XyzList[[3]][l]}
(* positive lobes only *)
RP[l_] := If[R[l]<0, 0, R[l]] GP[l_] := If[G[l]<0, 0, G[l]] BP[l_] := If[B[l]<0, 0, B[l]]
RPM[l_] := If[l<525, 0, RP[l]] (* major positive lobe only *)
(* inverses of xyz to rgb *)
XI[l_] := trx[[1]] . {R[l], G[l], B[l]} YI[l_] := trx[[2]] . {R[l], G[l], B[l]} ZI[l_] := trx[[3]] . {R[l], G[l], B[l]}
plotRGB := Plot[{R[l], G[l], B[l]}, {l,380,760}] plotXYZ := Plot[{X[l], Y[l], Z[l]}, {l,380,760}] plotXYZI := Plot[{XI[l], YI[l], ZI[l]}, {l,380,760}]
(* integrals of fundamentals *)
intXYZ := {NIntegrate[X[l], {l,380,760}], NIntegrate[Y[l], {l,380,760}], NIntegrate[Z[l], {l,380,760}]} intXYZI := {NIntegrate[XI[l], {l,380,760}], NIntegrate[YI[l], {l,380,760}], NIntegrate[ZI[l], {l,380,760}]} intRGB := {NIntegrate[R[l], {l,380,760}], NIntegrate[G[l], {l,380,760}], NIntegrate[B[l], {l,380,760}]} intRGBP := {NIntegrate[RP[l], {l,380,760}], NIntegrate[GP[l], {l,380,760}], NIntegrate[BP[l], {l,380,760}]} intRGBPM := {NIntegrate[RPM[l], {l,380,760}], NIntegrate[GP[l], {l,380,760}], NIntegrate[BP[l], {l,380,760}]}
(* compute integrals of primaries specified as an XYZ->RGB transform *)
NIntTxr[m_, CieFns_:XyzList] := ( txr = m; XyzList = CieFns; intRGB )
NIntTxrP[m_, CieFns_:XyzList] := ( txr = m; XyzList = CieFns; intRGBP )
NIntTxrPM[m_, CieFns_:XyzList] := ( txr = m; XyzList = CieFns; intRGBPM )
(* The following function computes a transformation from Cie XYZ color matching functions to RGB color matching functions, given a set of primaries specified by chromaticity coordinates only, and the chromaticities of the alignment white. This function has been verified against the results given in Rogers 85; the only difference is the normalization factor. This method seems a lot easier than Rogers'. *)
TxrFromChrom[{{Rx_, Ry_}, {Gx_, Gy_}, {Bx_, By_}}, {Wx_,Wy_}, CieFns_:XyzList] := (* Compute an XYZ to RGB transform from RGB chromaticity coordinates, equalized on the white point given by chromaticity coordinates. The result is also normalized for max[rgb]=1. *) Module[{Rz, Gz, Bz, Wz, m, sf}, ( Rz = 1 - Rx - Ry; Gz = 1 - Gx - Gy; Bz = 1 - Bx - By; Wz = 1 - Wx - Wy; m := Inverse[{{Rx, Gx, Bx}, {Ry, Gy, By}, {Rz, Gz, Bz}}]; (* equalize on white *) m = (1 / (m . {Wx,Wy,Wz})) m; (* scale for max[rgb]=1 *) m = RgbNormFactorFromTxr[m, CieFns] m; m )]
(* This function computes chromaticity coordinates of primaries, given a transform from Cie XYZ color matching functions to the color matching functions of those primaries. Note color matching functions are NOT the same as SPD's, hence you cannot use the technique of function ChromFromSPD here. In fact the computation is a lot simpler. Use this function as the inverse of function TxrFromChrom. This function has been verified against the results given in Rogers85, p. 395 ff. *)
ChromFromTxr[txr_] := Module[ {trx = Inverse[txr]}, ( Map[{#[[1]] / (#[[1]] + #[[2]] + #[[3]]), #[[2]] / (#[[1]] + #[[2]] + #[[3]])}&, {trx . {1,0,0}, trx . {0,1,0}, trx . {0,0,1}}] )]
(* Compute chromaticity coordinates from SPD's as fn of wavelength. This function has been verified on Cie sources A, B, C, D65, and E. The results are virtually the same as the chromaticities found in the literature. *)
ChromFromSPD[plist_, CieFns_:XyzList] := Module[{F, result, sum}, ( XyzList = CieFns; Map[ Function[ el, ( result = Map[ Function[x, NIntegrate[x, {l,380,760}]], {F[l] XyzList[[1]][l], F[l] XyzList[[2]][l], F[l] XyzList[[3]][l]} /. F->el ]; sum = Plus @@ result; Take[result / sum, 2] ) ], plist ] ) ]
(* some typical wavelengths for spectral colors *)
{LBlue, LGreen, LYellow, LRed} = {475, 500, 580, 700}
(* for displaying fundamentals *)
PlotTxr[m_,label_:"",style_:rgbColors,CieFns_:XyzList] := ( txr = m; XyzList = CieFns; Plot[{R[l], G[l], B[l]}, {l,380,760}, PlotRange->All, PlotLabel->label, PlotStyle->style] )
PlotTxrP[m_,label_:"",style_:rgbColors,CieFns_:XyzList] := ( txr = m; XyzList = CieFns; Plot[{RP[l], GP[l], BP[l]}, {l,380,760}, PlotRange->All, PlotLabel->label, PlotStyle->style] )
PlotTxrPM[m_,label_:"",style_:rgbColors,CieFns_:XyzList] := ( txr = m; XyzList = CieFns; Plot[{RPM[l], GP[l], BP[l]}, {l,380,760}, PlotRange->All, PlotLabel->label, PlotStyle->style] )
PlotTxo[m_,label_:"",style_:rgbColors,CieFns_:XyzList] := ( txr = m; XyzList = CieFns; Plot[{R[l], G[l], B[l]}, {l,380,760}, PlotRange->All, PlotLabel->label, PlotStyle->style] )
(* displaying Cie diagram with colored outline *)
(* this function assigns RGB colors to wavelengths, scaling them up to maximum intensity while preserving rgb hue *)
lcolor[l_] := Module[ {r,g,b}, ( {r,g,b} = {RPM[l], GP[l], BP[l]}; {r,g,b} = (1 / Max[r,g,b]) {r,g,b}; RGBColor[r,g,b] )]
(* only ParametricPlot3D seems to take a color argument *) (* the purple line does not show for some reason *) (* delayed evaluation because XyzToRgbCRT has to be defined first *)
ChromDiagramColor := ( txr = transforms`XyzToRgbCRT; XyzList = Cie31List; Show[ParametricPlot3D[{x[l], y[l], 0, lcolor[l]}, {l,380,750}, PlotRange->{{0,.80}, {0,.85}, {-.001,.001}}, ViewPoint->{0,0,25}, Axes->{True,True,False}, DisplayFunction->Identity], Graphics3D[{RGBColor[1,0,1], Line[{{x[380], y[380],0}, {x[750], y[750],0}}]}], Background->RGBColor[.5,.5,.5], DisplayFunction->Identity] )
cPointColor[{x_, y_}] := Graphics3D[{GrayLevel[0], PointSize[0.02], Point[{x,y,0}]}] cPointLColor[l_] := cPointColor[chrom[l]]
PlotChromListColor[list_] := Show[ChromDiagramColor, Map[cPointColor, list], DisplayFunction->$DisplayFunction]
(* Plot RGB gamuts, given a transformation from XYZ coordinates to RGB coordinates. Note that when displaying gamuts for transforms other than XyzToRgbCRT, the acutally displayed colors are not correct since they are limited to what the monitor can display (of which CRT is an approximation). The extent of the gamut is ok, however. The gamuts are displayed in the Cie chromaticity space (xy coordinates), so colors are normalized for maximum intensity. To plot CRT and LMS gamuts respectively, use these expressions: Show[Txr2ChromGamut[transforms`XyzToRgbCRT, 0.85, 64, "CRT"], DisplayFunction->$DisplayFunction] Show[Txr2ChromGamut[transforms`XyzToLmsVW, 0.85, 64, "LMS"], DisplayFunction->$DisplayFunction] *)
(* The function below uses Plot3D only because that allows colors of plotted points to be specified, and Plot doesn't. It plots a plane at z coordinate 0, and x,y determined by chromaticity coordinates. The range for the regular Cie chromaticity diagram is 0.85 *)
Txr2ChromGamut[txr_,range_,points_,label_,bg_:{.5,.5,.5},CieFns_:XyzList] := (XyzList = CieFns; Plot3D[{0, Chrom2RGBColorClip[{x,y}, txr, bg, CieFns]}, {x,0,range}, {y,0,range}, ViewPoint->{0,0,25}, Background->ReplacePart[bg, RGBColor, 0], PlotPoints->{points,points}, Lighting->False, Mesh->False, PlotLabel->label, Axes->{True,True,False}, DisplayFunction->Identity] )
(* The following functions compute RGB values from chromaticity coordinates and an XYZ->RGB transform (txr). The idea is that since chrom coords are related to XYZ coordinates by a constant 1/(X+Y+Z), we can use the chrom coords in stead of XYZ coords as input to the transform, and then scale the result so that the maximum component is 1. Since chrom coords carry no luminance info anyway, that is fine. The result should be gamma corrected for display of course, but Mathematica takes care of that. *)
Chrom2RGBColor[{x_,y_}, txr_, bg_:{.5,.5,.5}, CieFns_:XyzList] := ( XyzList = CieFns; RGBColor[#[[1]],#[[2]],#[[3]]]& [ (1/Max[#])#&[ Map[If[#<0,0,#]&, txr . {x, y, 1-x-y}]]] )
(* Same, but turns any point with a negative coordinate into bg gray. Displaying points computed in this way on the same gray bg makes only the valid points stand out on a gray bg. *)
Chrom2RGBColorClip[{x_,y_}, txr_, bg_:{.5,.5,.5}, CieFns_:XyzList] := ( XyzList = CieFns; RGBColor[#[[1]],#[[2]],#[[3]]]& [ If[Min[#]<0, bg, #]& [ (1/Max[#])#&[ txr . {x, y, 1-x-y}]]] )
(* These are some examples of gamuts in chromaticity coordinates *)
GamutNTSC := Txr2ChromGamut[transforms`XyzToRgbNTSC, 64, "NTSC"] GamutCRT := Txr2ChromGamut[transforms`XyzToRgbCRT, 64, "CRT"] GamutLMS := Txr2ChromGamut[transforms`XyzToLmsVW, 64, "LMS"]
(* Plot Opponent gamuts, given a transformation from opponent coordinates to RGB coordinates. The same restriction to displayed colors applies as above. Opponent gamuts are displayed at constant brightness planes through the opponent color space, with or without intensity normalization. *)
(* normalized *)
Tor2GamutN[tor_, brightness_, hrange_, vrange_, points_, label_] := Plot3D[{brightness, Tor2RGBColorClipN[{gr, by, brightness}, tor]}, {gr, -hrange, hrange}, {by, -hrange, hrange}, PlotPoints -> {points, points}, PlotRange -> {{-hrange, hrange}, {-hrange, hrange}, {0, vrange}}, Background -> GrayLevel[0.5], Lighting -> False, Mesh -> False, PlotLabel -> label, BoxRatios -> {1, 1, 2}, SphericalRegion -> True, DisplayFunction -> Identity]
(* The function for converting to rgb values has to limit display to possible points only; since the rgb transforms are scaled for a maximum value of one, any point which has an RGB component that is negative or exceeds 1 is not possible. Whatever transform is used, e.g. H81ToLms (human cone sensitivities), results are displayed on a monitor using its own gamut of course. Using Lms, good ranges are [0,1.8] for brightness and [-3.5,3.5] for gr/by. *)
Tor2RGBColorClipN[{gr_,by_,wb_}, tor_] := RGBColor[#[[1]],#[[2]],#[[3]]]& [ Which[ Min[#]<0, {.5,.5,.5}, Max[#]>1, {.5,.5,.5}, True, (1/Max[#])# ]& [ tor . {gr,by,wb} ]]
(* unnormalized *)
Tor2Gamut[tor_, brightness_, hrange_, vrange_, points_, label_] := Plot3D[{brightness, Tor2RGBColorClip[{gr, by, brightness}, tor]}, {gr, -hrange, hrange}, {by, -hrange, hrange}, PlotPoints -> {points, points}, PlotRange -> {{-hrange, hrange}, {-hrange, hrange}, {0, vrange}}, Background -> GrayLevel[0.5], Lighting -> False, Mesh -> False, PlotLabel -> label, BoxRatios -> {1, 1, 2}, SphericalRegion -> True, DisplayFunction -> Identity]
Tor2RGBColorClip[{gr_,by_,wb_}, tor_] := RGBColor[#[[1]],#[[2]],#[[3]]]& [ Which[ Min[#]<0, {.5,.5,.5}, Max[#]>1, {.5,.5,.5}, True, # ]& [ tor . {gr,by,wb} ]]
(* use these expressions to generate an animated plot of constant-brightness planes through H81 and WWAW opponent space, mapped into RgbCRT space. Note that using Lms results in different shapes. *)
AnimH81CRT := ShowAnimation[ Table[ Tor2Gamut[H81ToRgbCRT, b, 1, 2.5, 16, ""], {b,0,2.3,.1}]]
AnimH81Lms := ShowAnimation[ Table[ Tor2Gamut[H81ToLmsVW, b, 3, 2, 16, ""], {b,0,1.84,.08}]]
AnimWWAWCRT := ShowAnimation[ Table[ Tor2Gamut[WWAWToRgbCRT, b, 1, 1.5, 16, ""], {b,0,1.38,.06}]]
End[] (* private *) EndPackage[] (* package *)