Analogy of Proportion

Chromatic adaptation can be viewed as an Aristotelian analogy of proportion, see [1]. A general analogy of this type is usually written A : B = C : D and read as “A is to B as C is to D”. In our case the expression A : I is interpreted as “the color appearance of A in a viewing environment with illuminant I after the viewer is adapted to I”, or more simply “the appearance of A under illuminant I”. It is better to think of A not as an object color, but as a self-luminous color. The analogy A : I = B : J can be read as “the appearance of A under illuminant I is the same as the appearance of B under illuminant J”. This phenomenon is called chromatic adaptation and takes place in the human eye and brain. Of course, the exact values of A and B vary depending on the individual, so it only makes sense to think of averages, for a standard observer. All of these variables are typically XYZ tristimulus vectors, and that is what we will assume from now on. We require that all components of I and J are positive, i.e. the illuminants are in the positive octant.

We will adopt these axioms throughout this vignette: In A2) the first equality is equality of appearance, and the second is equality of vectors.

In most cases A, I, and J are known, and B is not. Solving the chromatic adaptation problem is solving the analogy A : I = X : J for X, where A is the source color, I is the the source illuminant (also called the test illuminant), and J is the target illuminant (also called the reference illuminant). By axiom A2) the target color solution X is unique.

Chromatic Adaptation Transforms

In most practical cases, illuminants I and J are fixed, and the source color A is allowed to vary widely. Recall that X is uniquely determined by A, so in the analogy A : I = X : J it makes sense to think of X as a transform T(A) of A. Define the ideal Chromatic Adaption Transform T from I to J by the analogy:

If T(A) : J = T(B) : J then A : I = B : I and so A = B, by (1.1). This shows that T is injective.

It is convenient to abbreviate Chromatic Adaption Transform by the acronym CAT. In theory the exactly accurate ideal CAT from I to J is unique, but in practice there are so many confounding variables and complexities that we cannot hope to find it. So from now on we open things up and look for approximations to T using various methods (see the next section).

There are 2 special cases where we can assert something exactly. Since I : I = J : J by axiom A1), we can say that T(I) = J exactly. And since A : I = A : I, we can say that if J = I, then T is the identity transform. In both theory and practice, we obviously want these 2 properties to be true, so we require that all approximations have them.

CAT Methods

1. identity. If J = I, T_II is the identity transform
2. inverse. T_JI is the inverse of T_IJ
3. commutative triangle. Adapting from I to J, and then from J to K, is the same as adapting from I to K. In symbols: T_JK ∘ T_IJ = T_IK (where ∘ denotes transform composition).

Property 1 (identity) is required by the last paragraph of the previous section. We now argue that properties 2 and 3 are desirable.

Property 2 (inverse) is true for an ideal CAT method, by the following argument: Swapping I and J gives T_IJ ∘ T_JI = identity, and this proves property 2. According to Hunt [4] p. 591, the non-linear version of the CMCCAT97 method does NOT have property 2. The degree to which property 2 fails is therefore a measure of the accuracy of the appearance matching of the CMCCAT97 CATs. In Burns [2] is a CAT method which in the original and direct form does not have the inverse property, but can be slightly modified to a “symmetric” version that does enjoy the property.

Property 3 (commutative triangle) is desirable because it makes it possible to use the intermediate Profile Connection Space (PCS with Illuminant D50) in ICC color profiles with no loss of accuracy. In property 3, J represents the PCS in the middle, and the endpoints I and K are viewing illuminants for devices. If property 3 fails, then the 2-step adaptation may be sub-optimal.

Remark 3.1 Assuming that each T_II is invertible (a very mild condition), properties 1. and 2. follow from 3. Letting J = K = I in property 3 we get T_II ∘ T_II = T_II. Since T_II is invertible we conclude T_II = identity and this proves property 1. Letting K = I in property 3 we get Swapping I and J gives and this proves property 2. Note that property 3 is similar to the cocycle condition in the construction of fiber bundles, see [9] and p. 14 of [8]. That property 3 implies properties 1 and 2 is on p. 8 of [8].

Remark 3.2 Property 1 does NOT follow from property 2. Here is a CAT method that is a counter-example. First suppose that the vectors I and J have the same length. Define T_IJ to be the transform that first rotates space π radians around I and then rotates by the most direct rotation from I to J. If I and J do not have the same length, then rotate I to the line spanned by J and then follow this by a uniform scaling so the result expands or shrinks to J. By “most direct rotation” we mean the rotation around the axis I × J (the 3D cross product). This T_IJ is perfectly well-defined (unless J = −I which is impossible since both are in the positive octant), and maps I to J. It is not hard to show that T_JI ∘ T_IJ = T_IJ ∘ T_JI = identity, so property 2 is true. But T_II is rotation of π radians around I, so property 1 is false.

The Simplest CAT Method - XYZ Scaling

A CAT method is called linear if and only if every T_IJ is a linear map (from ℝ³ to ℝ³). From this point on in the vignette all CAT methods are linear. T_IJ now denotes a 3x3 matrix. In ICC v4.0 color profiles this matrix is saved in the chromaticAdaptationTag or chad tag, and converts from device white to PCS white, see [5]. To avoid confusion between matrices and vectors, we shift notation and replace I, J, and K by 3-vectors u, v, and w (since lower-case i, j, and k look too much like integers). We also use the bold 1 := (1, 1, 1) for the 3-vector of all 1s, which is also a valid whitepoint.

1. identity. T_uu = I
2. inverse. T_vu = T_vu⁻¹
3. commutative triangle. T_vwT_uv = T_uw

The rest of this section follows Lindbloom [7]. Given 3-vectors u and v with all entries positive, define the XYZ scaling CAT method by: Here and through the rest of this vignette, the vector division v/u is component by component, i.e. the Hadamard division, see [10]. One easily verifies that: so this CAT method satisfies property 3. Since T_uu is trivially invertible, it also satisfies properties 1 and 2 by Remark 3.1. Each channel in XYZ is scaled independently. When this technique is used in electronic RGB cameras, and applied to the RGB channels independently, it is often called white-balancing. This CAT method is what one would get by transforming to XYZ to Lab using u as the whitepoint, and then from Lab to XYZ using v as the whitepoint (though this view hides the linearity of T_uv). This CAT method is sometimes called the “Wrong von Kries”, see [3].

Von-Kries-Based CAT Methods

All the linear CAT methods in common use are von-Kries-based which is now defined, loosely following Lindbloom [7]. Let M_a be an invertible 3x3 matrix that does not move 1 too much, i.e. M_a1 ≈ 1. First define One easily checks that T_u1u = 1 and T_1u1 = u and these matrices are inverses of each other. Now define a von-Kries-based CAT method in general by: One easily checks that so the method satisfies property 3. By Remark 3.1 it also satisfies properties 1 and 2.

Figure 5.1 - Commutative Diagram of CATs

Since the 3 inner triangles commute, the outer triangle commutes too. The symbol ℝ_u³ denotes the space of all possible XYZs under illuminant u. We are being sloppy here because a CAT makes no sense when one of the XYZs is negative; but since the CATs in the figure are linear maps, they can still be defined mathematically even if they make no physical sense.

Note that the dependence of the transforms on M_a is suppressed in this notation. In the von Kries theory, M_a transforms from XYZ to the cone response domain. In practice M_a is calculated from a large number of pairs of experimentally measured corresponding colors. For example the popular Bradford M_a was calculated from 58 pairs, see [6]. M_a is called the cone response matrix for the von-Kries-based CAT method. In ArgyllCMS ICC color profiles, M_a is saved in the SigAbsToRelTransSpace or arts private ICC tag, see [3].

In the special case M_a = I the transforms become: which is just the XYZ scaling method in the previous section. There is more going on. In words, equation (5.1) says: transform u and v by M_a and form the diagonal matrix one would get from XYZ scaling. To adapt color A from u to v, transform A by M_a, and then by the diagonal matrix, and then transform back again by M_a⁻¹. A fancy way to say the same thing is: a CAT matrix is von-Kries-based if and only if it is linearly conjugate to an XYZ scaling.

Since M_a has 9 parameters, it appears that these von-Kries-based CAT methods have 9 degrees of freedom. However, some changes of M_a do not change the transforms. Let D be a diagonal 3x3 matrix with positive entries on the diagonal, and let M′_a := DM_a. Then so changing M_a to DM_a leaves the CAT unchanged. Setting D = diag (M_a1)⁻¹ implies that M′_a1 = 1. The matrix diag (M_a1) is invertible because M_a1 ≈ 1 as we assumed above. Thus, any acceptable cone response matrix can be normalized so that M_a1 = 1 exactly, i.e. the row sums are all 1. So these von-Kries-based CAT methods have 9-3 = 6 degrees of freedom. Geometrically, M_a leaves the line generated by 1 pointwise fixed, so we think of it as a “twist” around the line (more general than a rotation around the line). If {1, b₂, b₃} is a basis of ℝ³, then M_a can map b₂ and b₃ to ℝ³ almost arbitrarily, and this gives the 3+3 = 6 degrees of freedom.

Other Linear CAT Methods

This section is a detour that asks: “Are there are linear CAT methods that are NOT von-Kries-based ?” Note from (5.1) that if T_uv is a von-Kries-based matrix, then T_uv is diagonalizable. So if we can construct a linear CAT method whose matrices are not diagonalizable, then it is not von-Kries-based.

As an example, define T_uv to be the most direct rotation of u to the ray generated by v, followed by uniform scaling to take u to v. Since all the non-trivial T_uv are products of rotations and scalings, they are not diagonalizable. Note that this CAT method has properties 1 and 2, but NOT property 3. It shows that properties 1 and 2 do NOT imply property 3.

Taking an idea from the previous section, we can modify this method to have property 3. Define T_u1 to be the most direct rotation of u to the ray generated by 1, then followed by a uniform scaling that makes u map to 1. Define T_1u to be the inverse of T_u1. Now define T_uv := T_1vT_u1 as before. All three properties in section CAT Methods are satisfied by the same arguments from the previous section. This CAT method may perform poorly. In Burns [2] the problem of negative tristimulus values is discussed, and this type of CAT method might make make the problem worse. I am sure that many other interesting (but impractical) examples can be constructed.

CAT Methods in Package spacesXYZ

Finally, we explore the CAT methods available in software.

library( spacesXYZ )

There is an S3 class CAT with constructor CAT(). The constructor takes arguments the source and target illuminant XYZ (denoted u and v above), and the CAT method. All the available non-trivial methods - Bradford, MCAT02, vonKries, and Bianco-Schettini - are von-Kries-based. For MCAT02 only the simple linear variant is used. The CAT object is a list with cone response matrix Ma (denoted M_a above), the adaptation matrix M (denoted T_uv above), and other things, see the CAT() man page.

Ma = CAT( source='A', target='D65', method='bradford' )$Ma ;  Ma

##         X       Y       Z
## L  0.8951  0.2664 -0.1614
## M -0.7502  1.7135  0.0367
## S  0.0389 -0.0685  1.0296

This is the famous Bradford cone-response-matrix, appearing in Lam [6] p. 3-46.

rowSums( Ma )

##      L      M      S 
## 1.0001 1.0000 1.0000

It appears that an attempt was made to normalize the row sums to 1, but roundoff made the last digit in the first row off by 1. There is no practical effect. The actual adaptation matrix is easily inspected and tested too:

theCAT  = CAT( source='A', target='D65', method='bradford' )
A  =  standardXYZ('A')
A %*% t(theCAT$M) - standardXYZ('D65')

##              X Y Z
## A 1.110223e-16 0 0

So M maps the XYZ of Illuminant A to that of D65 as required. Using explicit matrix multiplication is OK, but the function adaptXYZ() is preferred:

identical(  adaptXYZ( theCAT, A ), A %*% t(theCAT$M) )

## [1] FALSE

We can also inspect the row sums for method MCAT02.

rowSums( CAT( source='A', target='D65', method='MCAT02' )$Ma )

## L M S 
## 1 1 1

So for MCAT02 the normalization was more careful about roundoff. And for vonKries we have:

rowSums( CAT( source='A', target='D65', method='vonKries' )$Ma )

##       L       M       S 
## 1.02703 0.98472 0.91822

So for vonKries there was no normalization.

References

Appendix A - Recovery of M_a

Suppose we know that a chromatic adaption matrix T_uv is von-Kries-based and we also know the white points u and v. Can we recover the cone response matrix M_a ? This appendix presents a recipe to do that. We saw earlier that M_a is not unique; its rows are only defined up to a constant. So our solution will follow common practice and scale the rows to have sum 1.

Rearrange equation (5.1) to the form: We see that the columns of M_a^⊺ - equivalently the rows of M_a - are the eigenvectors of T_uv^⊺. Amazingly, the eigenvectors of T_uv^⊺ do not depend on u and v ! The eigenvectors are also only defined up to a constant, so rescaling them to have sum 1 is easy. But the eigenvectors (and eigenvalues) are also only defined up to a permutation; finding the right permutation is a little harder and is where u and v are used.

Here is a worked out example.

whiteA = standardXYZ("A")[1, ]  ;  whiteB = standardXYZ("B")[1, ]
theCAT = CAT( whiteA, whiteB, method='MCAT02' )
T   = theCAT$M ; Ma = theCAT$Ma
res = eigen( t(T) )
X   = t(res$vectors)  ;  X = diag( 1 / rowSums(X) ) %*% X  # X is 'first cut' at the unknown Ma

Compare Ma and X

Ma ; X

##         X      Y       Z
## L  0.7328 0.4296 -0.1624
## M -0.7036 1.6975  0.0061
## S  0.0030 0.0136  0.9834

##         [,1]   [,2]    [,3]
## [1,]  0.0030 0.0136  0.9834
## [2,] -0.7036 1.6975  0.0061
## [3,]  0.7328 0.4296 -0.1624

One can check that the row sums of Ma are 1. Now compare the white ratios and the eigenvalues:

as.numeric(Ma %*% whiteB / Ma %*% whiteA) ; res$values

## [1] 0.8643826 1.0850937 2.3297865

## [1] 2.3297865 1.0850937 0.8643826

So the row orders are reversed, and the eigenvalues too. The two permutations are the same: c(3,2,1).
In practice Ma is the unknown so we do not have it, but we do have the eigenvalues res$values. Since we know that res$values are returned in decreasing order, we can compute the desired permutation like this:

perm = order( Ma %*% whiteB / Ma %*% whiteA,  decreasing=TRUE )  ; perm

## [1] 3 2 1

The permutation perm now takes the ratios to the eigenvalues, but we want the inverse, so invert perm:

perm = order(perm) ; perm   

## [1] 3 2 1

res$values[perm]

## [1] 0.8643826 1.0850937 2.3297865

X = X[perm, ]  ;  X  ;  max( abs(X - Ma) )

##         [,1]   [,2]    [,3]
## [1,]  0.7328 0.4296 -0.1624
## [2,] -0.7036 1.6975  0.0061
## [3,]  0.0030 0.0136  0.9834

## [1] 2.220446e-15

We have recovered M_a with good accuracy.

Now suppose u and v are not known, but M_a is known to be in a small list of candidate matrices. Since there only 6 possible permutations that take matrix X to the right candidate, it will be easy to spot the right one.

Session Information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=C               LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] spacesXYZ_1.3-0 rmarkdown_2.28 

loaded via a namespace (and not attached):
 [1] digest_0.6.37        R6_2.5.1             microbenchmark_1.5.0 fastmap_1.2.0       
 [5] xfun_0.49            maketools_1.3.1      cachem_1.1.0         knitr_1.48          
 [9] htmltools_0.5.8.1    buildtools_1.0.0     lifecycle_1.0.4      cli_3.6.3           
[13] sass_0.4.9           jquerylib_0.1.4      compiler_4.4.1       sys_3.4.3           
[17] tools_4.4.1          evaluate_1.0.1       bslib_0.8.0          yaml_2.3.10         
[21] jsonlite_1.8.9       rlang_1.1.4