Supplementary Documentation


Published February 27, 2019 by Scott Allen Burns; last updated May 8, 2020

Change Log
2/27/2019 Original publication.
4/23/2019 Added section on symmetry.
8/1/2019 Added citation to published paper in Color Research & Application.
2/11/2020 Updated the symmetry section to expand on how it can be derived.
3/26/20 Added change log details.
5/8/20 Minor correction of location of update dates. Update of links to publications.


Introduction

This page contains supplementary documentation for the article “Chromatic Adaptation Transform by Spectral Reconstruction,” published in Color Research and Application in 2019 (Burns SA. Chromatic adaptation transform by spectral reconstruction. Color Res Appl. 2019;44(5):682-693). Any reference to “the Article” on this page refers to this document. The documentation that follows is organized according to the section numbers appearing in the Article. A preprint version of this paper is also available here.

Please note that the content of this page should not be considered peer reviewed by the journal, as it may be updated at any time. This page was last updated on the date shown above.


2. Examples of CATs Producing Negative Tristimulus Values

An instance of CAT02 producing a destination color with negative tristimulus values is presented here, specifically, for the case of Munsell 2.5G 4/8. The spectrophotometrically measured reflectance curve for this color is:

Munsell 2.5 4/8 reflectance (click to enlarge).


A purplish source illuminant is used, shown in the following table. Equal energy is used as the destination illuminant. Before applying CAT02, the source color tristimulus values must be obtained under the action of the illuminant. To do this, the 1931 CMFs are referenced to the source illuminant, W_S:

Computation of As (click to enlarge).


Next, the source tristimulus values, XYZ_S, are computed from the reflectance curve for Munsell 2.5G 4/8, \rho_{(2.5G4/8)}, and A_S:

(1)   \begin{equation*} {XYZ}_S = {A_S}'\;\rho_{(2.5G4/8)} = \begin{bmatrix} 0.0928 \\ 0.0784 \\ 0.2220 \end{bmatrix}. \end{equation*}

Forming the CAT02 matrix starts with the transformation from tristimulus values to sharpened cone response space, M_{CAT02}:

(2)   \begin{equation*} M_{CAT02} = \left[\begin{array}{rrr} 0.7328	&0.4296	&-0.1624 \\ -0.7036	&1.6975	&0.0061\\ 0.0030	&0.0136	&0.9834 \end{array}\right]. \end{equation*}

The source and destination white points are XYZ_{wps} = (2.145, 1, 4.729) and XYZ_{wpd} = (1.000, 1.000, 1.000), which when pre-multiplied by M_{CAT02} give white point cone responses RGB_{wps} = (1.2338, 0.2168, 4.6700) and RGB_{wpd} = (1.000, 1.000, 1.000). The diagonal von Kries coefficient matrix (using full adaptation, D=1) comprises ratios of destination to source white point cone responses:

(3)   \begin{equation*} \Lambda = \left[\begin{array}{ccc} 0.8104 & 0 & 0 \\ 0 & 4.6126 & 0 \\ 0 & 0 & 0.2140 \end{array}\right]. \end{equation*}

The final CAT02 matrix can now be computed as

(4)   \begin{equation*} \mathrm{CAT02} = M_{CAT02}^{-1}\;\Lambda\;M_{CAT02} = \left[\begin{array}{rrr} 1.5561&-1.8013&-0.1137 \\ -1.2669&3.8661&-0.0313 \\ 0.0134&-0.0450&0.2147 \end{array}\right]. \end{equation*}

The predicted destination color, XYZ_{Dp}, produced by CAT02 from the source tristimulus values, XYZ_S, ends up with negative tristimulus values:

(5)   \begin{equation*} XYZ_{Dp} = \mathrm{CAT02} \left[ \begin{array}{r} 0.0928 \\ 0.0784 \\ 0.2220 \end{array} \right] = \left[ \begin{array}{r} -0.0220 \\ 0.1785 \\ 0.0454 \end{array} \right]. \end{equation*}


4. Chromatic Adaptation Transform by Spectral Reconstruction

Table 2 in the Article presents \Delta E_{94}^\ast color differences between the experimental corresponding-colors dataset destination colors, XYZ_D, and the CAT-predicted destination colors, XYZ_{Dp}, for four methods: HPE, CAT02, CAT16, and the proposed CAT. In all cases, full adaptation (D=1) is used.

To confirm the accuracy of these computations, they were compared to instances of the same computation by other authors in the literature, and were found to match exactly the published D=1 values of:

\texttt{Table 4} of Reference: Susstrunk S, Holm J, Finlayson GD. Chromatic adaptation performance of different RGB sensors. Proc. SPIE 4300, Color Imaging: Device-Independent Color, Color Hardcopy, and Graphic Arts VI. Reiner Eschbach R, Marcu GG, Editors. 2001:172-183.

\texttt{Table AII} of Reference: Bianco S, Schettini R. Two new von Kries based chromatic adaptation transforms found by numerical optimization. Color Res Appl. 2010;35(3):184-192.

\texttt{Table II (CMCCAT-McC)} of Reference: Hunt RWG, Li CJ, Juan LY, Luo MR. Further improvements to CIECAM97s. Color Res Appl. 2002;27(3):164-170.

Published D=1 color difference values for CAT16 have not been found in the literature, for verification purposes.


5. Spectral Reconstruction Computation

Click here to open a text file containing a numerical example of the “spectral_recon” function shown in Section 5 of the Article, applied to XYZ = (0.7,0.6,0.5) using equal energy as the reference illuminant. It is best viewed by copying the text into an editor like Notepad and turning off word wrap. The values of each assignment statement on each line of the code are shown so it can be used to debug an implementation of the method. There may be small differences from one implementation to the next, mostly due to differences in floating point precision and different implementations of linear equation solving among the various utilities.

The author has also developed additional algorithms for reconstructing reflectance curves, including some that put an upper bound on the curve, \rho \le 1. See this web page for more information on these additional methods and comparisons to natural paints and pigments, including Munsell color chips. An application of these spectral reconstruction methods for subtractive color mixture can be found here.


6. Implementing the CAT

Here are two MATLAB/Octave functions that implement the proposed CAT. Note that the matrix C used here (the finite-differencing constants array) is normally referred to as D in other presentations concerning reflectance reconstruction. The switch to C in the CAT presentation is due to the convention of D being used as the degree of adaptation.

The first function precomputes the necessary arrays for a given choice of source and destination illuminant white point tristimulus values (Y should equal 1 in both).

    \[ \begin{verbatim} function [C,As,Ad,Asd] = precomputed_arrays(A,XYZwps,XYZwpd) % A is nx3 matrix of color matching functions (CMFs) by columns. % XYZwps is 3x1 vector of source illuminant white point tristimulus values. % XYZwpd is 3x1 vector of destination illuminant white point tristimulus values. %        Note: XYZwps and XYZwpd should both have Y=1. % C is nxn matrix of finite-differencing constants. % As is nx3 matrix of Ws-referenced CMFs. % Ad is nx3 matrix of Wd-referenced CMFs. % Asd is nx3 matrix of Ws and Wd-referenced CMFs. n=size(A,1); C=full(gallery('tridiag',n,-2,4,-2)); C(1,1)=2; C(n,n)=2; Ws = spectral_recon(A,C,XYZwps); % See the Article for this function source code. Wd = spectral_recon(A,C,XYZwpd); % See the Article for this function source code. As=diag(Ws)*A; Ad=diag(Wd)*A; Asd=diag(Ws)*Ad; \end{verbatim} \]

The second function produces a predicted destination color, XYZdp, given a triplet of source tristimulus values, XYZs.

    \[ \begin{verbatim} function XYZdp=CAT_sym(C,As,Ad,Asd,XYZs) % Symmetric version of proposed CAT. % C is nxn matrix of finite-differencing constants. % As, Ad, and Asd are precomputed for a particular choice of %     source and destination illuminants. % As is nx3 Ws-referenced CMFs, As = diag(Ws)*A. % Ad is nx3 Wd-referenced CMFs, As = diag(Wd)*A. % Asd is nx3 Ws and Wd-referenced CMFs, Asd = diag(Ws)*Ad. % XYZs is 3x1 vector of source tristimulus values. % XYZdp is 3x1 vector of predicted destination tristimulus values. rho=spectral_recon_sym(As,Asd,C,XYZs); % See the Article for this function source code. XYZd=Ad'*rho; XYZdp=XYZd*XYZs(2)/XYZd(2); \end{verbatim} \]

Here is an example:

    \[ \begin{verbatim} >> % A is nx3 matrix of CMFs (1931 2 deg in this case). >> >> % Source illuminant A: >> XYZwps = [1.0981; 1.0000; 0.3555]; >> >> % Destination illuminant D65: >> XYZwpd = [0.9501; 1.0000; 1.0882] >> >> % Precompute arrays for illum A and D65: >> [C,As,Ad,Asd] = precomputed_arrays(A,XYZwps,XYZwpd); >> >> % Three examples of CAT application: >> XYZdp=CAT_sym(C,As,Ad,Asd,[.2;.3;.1]) XYZdp =     0.1699     0.3000     0.2415 >> XYZdp=CAT_sym(C,As,Ad,Asd,[.4;.5;.6]) XYZdp =     0.4393     0.5000     1.4401 >> XYZdp=CAT_sym(C,As,Ad,Asd,[1.0981; 1.0000; 0.3555]) XYZdp =     0.9501     1.0000     1.0882 \end{verbatim} \]

The last example demonstrates that achromatic source colors map to achromatic colors with respect to the destination illuminant, as expected.


7. Symmetry

To show how substituting some instances of A_S with A_{SD} makes the proposed CAT symmetric, we start with the simpler case of the non-log version of reflectance reconstruction, since it has a closed-form solution. The NLP being solved in this case is

(6)   \begin{equation*} \begin{split} \underset{\rho}{\mathrm{min}}\;\;&\tfrac{1}{2}\,\rho' C \rho\\ \mathsf{s.t.}\;\;&{A_S}'\;\rho = XYZ_S. \end{split} \end{equation*}

See Section 3 of the Article for definitions of the symbols above. Again, note that the symbol C is used here instead of the customary D for the finite-differencing matrix because of the convention in CAT literature of using D to represent the degree of adaptation. The stationary conditions arising from the method of Lagrange multipliers are

(7)   \begin{equation*} \begin{split} \partial \mathcal{L}/\partial \rho &= C\,\rho + A_S\,\lambda = 0 \\ \partial \mathcal{L}/\partial \lambda &= {A_S}'\,\rho - XYZ_S = 0. \end{split} \end{equation*}

Solving the first stationary condition for \rho gives

(8)   \begin{equation*} \rho = C^{-1} {A_S}\,\lambda. \end{equation*}

Substituting this expression into the second stationary condition and solving for \lambda yields

(9)   \begin{equation*} \lambda = {({A_S}'\,C^{-1} A_S)}^{-1} XYZ_S. \end{equation*}

Finally, substituting the expression for \lambda into the previous equation gives us a closed-form solution for the spectral reconstruction:

(10)   \begin{equation*} \rho = C^{-1} A_S\,{({A_S}'\,C^{-1} A_S)}^{-1} XYZ_S. \end{equation*}

To compute how this reconstructed reflectance would appear under a destination illuminant, we simply pre-multiply \rho by a destination-illuminant-referenced set of color matching functions, A_D:

(11)   \begin{equation*} {XYZ}_D = {A_D}'\,\rho = {A_D}'\,C^{-1} A_S\,{({A_S}'\,C^{-1} A_S)}^{-1} XYZ_S. \end{equation*}

This gives us a 3 \times 3 CAT matrix of the form

(12)   \begin{equation*} CAT = {A_D}'\,C^{-1} A_S\,{({A_S}'\,C^{-1} A_S)}^{-1}. \end{equation*}

Some comments are in order at this point.

1. The C matrix in the form as presented in the Article is singular, and so it is not suitable for this version of spectral reconstruction. However, a slight change to it, such as changing the first and last elements of the main diagonal to 4 instead of 2, makes the matrix invertible. The effect this change has on the reconstructed reflectance curve is to pull the ends of the curve to zero. This can be beneficial in some cases, as colored objects found in nature (e.g., common pigments; see Figure 8 in the Article) often have reflectances that drop to zero on the ultraviolet side, so the reconstruction is actually a better match on that end. The singularity of C will not be an issue in the subsequent development, as it will not be inverted in the end result (Equations 8 and 9 in the Article).

2. The CAT derived above needs only to be computed once for a source/destination illuminant pair. Once this 3 \times 3 CAT matrix is formed, each destination color computed from a source color requires only simple matrix multiplication, just as is the case with existing CATs. It will in fact perform very well when tested against the corresponding-color datasets. Unfortunately, it is capable of producing destination colors with negative tristimulus values, since the reconstructed reflectances can have negative values. The propensity for generating negative destination colors with this CAT is even more pronounced than it is with CAT02 and CAT16, and so it is not recommended for this reason. Nevertheless, it is a useful stepping stone in explaining the symmetric version of the method proposed in the Article.

Moving on to the issue of symmetry of the method, one way to express symmetry of the CAT is to describe mathematically the composite action of first transforming from source to destination illuminants, and then doing the reverse, transforming from destination to source illuminants. This composition should yield the identity matrix:

(13)   \begin{equation*} I = ({A_S}'\,C^{-1} A_D\,{({A_D}'\,C^{-1} A_D)}^{-1})\;\;({A_D}'\,C^{-1} A_S\,{({A_S}'\,C^{-1} A_S)}^{-1}). \end{equation*}

By replacing two instances of A_S and two instances of A_D with the dual-referenced A_{SD} we get

(14)   \begin{equation*} I = {A_S}'\,C^{-1} A_{SD}\;\;\boldsymbol{{({A_D}'\,C^{-1} A_{SD})}^{-1}\;\;{A_D}'\,C^{-1} A_{SD}}\;\;{({A_S}'\,C^{-1} A_{SD})}^{-1}. \end{equation*}

The central pair of expressions (in bold) now cancel, and then the remaining two expressions also cancel, leaving the identity matrix, as desired. The author provides no other justification for the change to A_{SD} other than the observation that it simply works to make the CAT symmetric without any apparent negative side effects.

If we trace backwards to locate where the change to A_{SD} could first have occurred, we find that if the stationary conditions are altered to be

(15)   \begin{equation*} \begin{split} \partial \mathcal{L}/\partial \rho &= C\,\rho + A_{SD}\,\lambda = 0 \\ \partial \mathcal{L}/\partial \lambda &= {A_S}'\,\rho - XYZ_S = 0, \end{split} \end{equation*}

then the symmetric form of the CAT is produced. This change only affects the Lagrange multiplier term in the first set of stationary conditions and does not alter the constraint expression, {A_S}'\,\rho = XYZ_S, which is necessary to ensure that the reconstructed reflectance has the required W_S-referenced tristimulus values.

Now turning to the log version of the spectral reconstruction used in the Article, an analogous treatment is applied to the first set of stationary conditions:

(16)   \begin{equation*} \begin{split} \partial \mathcal{L}/\partial z &= C\,z + \mathrm{diag}(e^z)\,A_{SD}\,\lambda = 0 \\ \partial \mathcal{L}/\partial \lambda &= {A_S}'\,e^z - XYZ_S = 0. \end{split} \end{equation*}

This change will result in the symmetric version of the proposed method, as presented in Equations 12 and 13 of the Article.


UPDATE 2/11/2020
Back when the paper was published, I indicated that the earliest point at which the switch to A_{SD} could occur was at the stationary conditions stage. I have since come to realize that it can be traced all the way back to a modified form of the original optimization statement. Considering the non-log version again, we replace

(17)   \begin{equation*} \begin{split} \underset{\rho}{\mathrm{min}}\;\;&\tfrac{1}{2}\,\rho' C \rho\\ \mathsf{s.t.}\;\;&{A_S}'\;\rho = XYZ_S. \end{split} \end{equation*}

with

(18)   \begin{equation*} \begin{split} \underset{\rho}{\mathrm{min}}\;\;&\tfrac{1}{2}\,\rho' C \; \overline{W_D}^{\,-1} \rho\\ \mathsf{s.t.}\;\;&{A_S}'\;\rho = XYZ_S, \end{split} \end{equation*}

The new term in the objective function is \overline{W_D}^{\,-1}, which is a diagonal matrix with the inverse of the destination illuminant on the main diagonal (and zeros elsewhere). If we construct the stationary conditions, we get

(19)   \begin{equation*} \begin{split} \partial \mathcal{L}/\partial \rho &= C\,\overline{W_D}^{\,-1}\,\rho + A_S\,\lambda = 0 \\ \partial \mathcal{L}/\partial \lambda &= {A_S}'\,\rho - XYZ_S = 0. \end{split} \end{equation*}

Notice that if we multiply the first of these equations by \overline{W_D}, we have

(20)   \begin{equation*} \begin{split} \partial \mathcal{L}/\partial \rho &= C\,\rho + \overline{W_D}\,A_S\,\lambda = 0 \\ \partial \mathcal{L}/\partial \lambda &= {A_S}'\,\rho - XYZ_S = 0. \end{split} \end{equation*}

The final step is to note that A_{SD} = \overline{W_D}\,A_S, which transforms the equations above to the form shown in Equation 15. The remainder of the explanation of the symmetric CAT follows from Equation 15 above.

It is natural to ask how we should interpret the new objective function of Equation 18. One way to do this would be to absorb the illuminant term into the C term, yielding some sort of illuminant-weighted finite-differencing matrix. Perhaps a better way to express it would be to split the illuminant term and equally absorb it into the two \rho terms:

(21)   \begin{equation*} \begin{split} \underset{\rho}{\mathrm{min}}\;\;&\tfrac{1}{2}\,(\overline{W_D}^{\,-1/2}\,\rho)' \;C \; (\overline{W_D}^{\,-1/2}\,\rho)\\ \mathsf{s.t.}\;\;&{A_S}'\;\rho = XYZ_S, \end{split} \end{equation*}

The objective function can now be interpreted as a sum of squares of differences:

(22)   \begin{equation*} (w_1^{\,-1/2} \rho_1 -\, w_2^{\,-1/2} \rho_2)^2 + (w_2^{\,-1/2} \rho_2 \,-\, w_3^{\,-1/2} \rho_3)^2 + ... + (w_{n-1}^{\,-1/2} \rho_{n-1} -\, w_n^{\,-1/2} \rho_n)^2, \end{equation*}

where w_i is the i^{th} term of illuminant W_D. In the symmetric version of the CAT, we are seeking an illuminant-weighted smoothest (least slope squared) reflectance reconstruction to account for symmetry.

END OF UPDATE


Table 4 of the Article demonstrates the symmetric version in comparison to the original version:


See the Article for the meaning of the various columns. The entries in this table are produced by the following process.

First, the MATLAB/Octave function for reconstructing a reflectance curve for the original version is:

    \[ \begin{verbatim} function rho=spectral_recon(As, C, XYZ) % As is nx3 Ws-referenced CMFs, As = diag(Ws)*A % C is nxn matrix of finite differencing constants (Eqn 5) % XYZ is a 3x1 vector of target Ws-referenced tristimulus values % rho is a nx1 vector of reflectance values (or zeros if failure) n=size(As,1); rho = zeros(n,1); z=zeros(n,1); lambda=zeros(3,1); count=0; maxit=20; % max number of iterations ftol=1.0e-8; % solution tolerance while count <= maxit     r=exp(z);     dr=diag(r);     v=dr*As*lambda;     F=[C*z+v; As'*r-XYZ];     J=[C+diag(v), dr*As; As'*dr, zeros(3)];     delta=J\(-F); % solve system of equations J*delta = -F     z=z+delta(1:n);     lambda=lambda+delta(n+1:n+3);     if all(abs(F)<ftol)         rho=exp(z);         return     end     count=count+1; end \end{verbatim} \]

and the corresponding similar function for the symmetric case is:

    \[ \begin{verbatim} function rho=spectral_recon_sym(As, Asd, C, XYZ) % As is nx3 Ws-referenced CMFs, As = diag(Ws)*A % Asd is nx3 Ws- and Wd-referenced CMFs, Asd = diag(Wd)*As % C is nxn matrix of finite differencing constants % XYZ is a 3x1 vector of target Ws-referenced tristimulus values % rho is a nx1 vector of reflectance values (or zeros if failure) n=size(As,1); rho = zeros(n,1); z=zeros(n,1); lambda=zeros(3,1); count=0; maxit=20; % max number of iterations ftol=1.0e-8; % solution tolerance while count <= maxit     r=exp(z);     dr=diag(r);     v=dr*Asd*lambda;     F=[C*z+v; As'*r-XYZ];     J=[C+diag(v), dr*Asd; As'*dr, zeros(3)];     delta=J\(-F); % solve system of equations J*delta = -F     z=z+delta(1:n);     lambda=lambda+delta(n+1:n+3);     if all(abs(F)<ftol)         rho=exp(z);         return     end     count=count+1; end \end{verbatim} \]

The main difference is the use of a dual-referenced set of CMFs in the symmetric version, A_{sd}. See the Article for definitions of all quantities. See the Article for MATLAB/Octave code to create matrix C.

The following MATLAB/Octave commands create the values in Table 4 for the non-symmetric case:

    \[ \begin{verbatim} >> % Illuminant white points (A & D65) >> XYZwps=[1.0981;1;0.3555] XYZwps =     1.0981     1.0000     0.3555 >> >> XYZwpd=[0.9501;1;1.0882] XYZwpd =     0.9501     1.0000     1.0882 >> >> % spectral reconstruction of illum >> Ws=spectral_recon(A, C, XYZwps); >> Wd=spectral_recon(A, C, XYZwpd); >> >> % source & dest illum-weighted CMFs >> As=diag(Ws)*A; >> Ad=diag(Wd)*A; >> >> % source color >> XYZs=[0.2;0.3;0.1] XYZs =     0.2000     0.3000     0.1000 >> >> % source spectral reconstruction >> rho=spectral_recon(As, C, XYZs); \end{verbatim} \]

    \[ \begin{verbatim} >> >> % unadjusted destination color >> XYZd=Ad'*rho XYZd =     0.2038     0.3581     0.2896 >> >> % adjusted destination color >> % and output of proposed CAT >> XYZdp=XYZd*(XYZs(2)/XYZd(2)) XYZdp =     0.1707     0.3000     0.2426 >> >> % next apply CAT from dest->source >> rho2=spectral_recon(Ad, C, XYZdp); >> XYZs2=As'*rho2 XYZs2 =     0.1716     0.2500     0.0847 >> >> adjustment >> XYZsp=XYZs2*(XYZdp(2)/XYZs2(2)) XYZsp =     0.2059     0.3000     0.1016 >> >> % Note that XYZsp doesn't match XYZs \end{verbatim} \]

The symmetric case in the Table has values generated by the following MATLAB/Octave commands:

    \[ \begin{verbatim} >> % Illuminant white points (A & D65) >> XYZwps=[1.0981;1;0.3555] XYZwps =     1.0981     1.0000     0.3555 >> XYZwpd=[0.9501;1;1.0882] XYZwpd =     0.9501     1.0000     1.0882 >> >> % spectral reconstruction of illuminant >> Ws=spectral_recon(A, C, XYZwps); >> Wd=spectral_recon(A, C, XYZwpd); >> >> % source & dest illum-weighted CMFs >> As=diag(Ws)*A; >> Ad=diag(Wd)*A; >> >> % dual referenced CMFs >> Asd=diag(Wd)*As; >> >> % source color >> XYZs=[0.2;0.3;0.1] XYZs =     0.2000     0.3000     0.1000 \end{verbatim} \]

    \[ \begin{verbatim} >> >> % source spectral reconstruction >> rho=spectral_recon_sym(As, Asd, C, XYZs); >> >> % unadjusted destination color >> XYZd=Ad'*rho XYZd =     0.2033     0.3589     0.2890 >> >> % adjusted destination color >> % and output of proposed symmetric CAT >> XYZdp=XYZd*(XYZs(2)/XYZd(2)) XYZdp =     0.1699     0.3000     0.2416 >> >> % next apply CAT from dest->source >> rho2=spectral_recon_sym(Ad, Asd, C, XYZdp); >> XYZs2=As'*rho2 XYZs2 =     0.1672     0.2507     0.0836 >> >> adjustment >> XYZsp=XYZs2*(XYZdp(2)/XYZs2(2)) XYZsp =     0.2000     0.3000     0.1000 >> >> % Note XYZsp matches XYZs \end{verbatim} \]