Color: Determining a Forward Matrix for Your Camera

We understand from the previous article that rendering color during raw conversion essentially means mapping raw data in the form of rgb triplets into a standard color space via a Profile Connection Space in a two step process

    \[ Raw Data \rightarrow  XYZ_{D50} \rightarrow RGB_{standard} \]

The first step white balances and demosaics the raw data, which at that stage we will refer to as rgb, followed by converting it to XYZ_{D50} Profile Connection Space through linear projection by an unknown ‘Forward Matrix’ (as DNG calls it) of the form

(1)   \begin{equation*} \left[ \begin{array}{c} X_{D50} \\ Y_{D50} \\ Z_{D50} \end{array} \right] = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \times \left[ \begin{array}{c} r \\ g \\ b \end{array} \right] \end{equation*}

Determining the nine a coefficients of this matrix M is the main subject of this article[1].

The second step projects the resulting image information from XYZ_{D50} to the ‘output’ colorimetric color space chosen by the photographer, say sRGB or Adobe RGB.  The necessary linear matrices for this transformation are standardized and readily available online.

9 Equations and 9 Unknowns

So how do we determine the nine coefficients of the Forward Matrix in Equation 1?

The 9 unknown coefficients operate on white balanced and demosaiced rgb data to transform it linearly into XYZ_{D50} data.  It follows that if we had the rgb values from three pixels and we knew their corresponding XYZ_{D50} triplets we could solve for the a coefficients.  The results would be valid for the given hardware, scene and illumination.

For instance we could capture in the raw data 3 patches of uniform diffuse reflectance illuminated by a known light source with the camera whose matrix we want to determine, thus obtaining three sets of rgb values; then measure with a spectrophotometer or similar instrument the reflectance of the 3 patches and illuminant spectral power distribution; and calculate the  XYZ_{D50} values that the reflectance and illuminant imply.

All that would be left to do then is assemble the 3×3 pairs  of rgb and corresponding XYZ_{D50} data  in the form of 9 equations and 9 unknowns to solve for the coefficients of the relative Forward Matrix M.  Below the a‘s that make up M are the unknowns, the X,Y,Z’s and r,g,b‘s would be known:

(2)   \begin{equation*} \begin{align*} X_1 &= a_{12}r_1 + a_{12}g_1 + a_{13}b_1 \\ Y_1 &= a_{21}r_1 + a_{22}g_1 + a_{23}b_1 \\ Z_1 &= a_{31}r_1 + a_{32}g_1 + a_{33}b_1 \\ X_2 &= a_{12}r_2 + a_{12}g_2 + a_{13}b_2 \\ \vdots \\ Z_3 &= a_{31}r_3 + a_{32}g_3 + a_{33}b_3 \\ \end{align*} \end{equation*}

Looks easy?    In fact it is much easier than that.

1 Capture, 72 Equations

X-Rite produces a number of ColorChecker 24 patch standard targets whose reflectance information is published and well thumbed.  For this example I will use their handy Passport Photo version.

The 24 patches in the ColorChecker target carry reflectances that are supposed to be representative of everyday photographic subjects, as found in skin, foliage and sky colors.  The color scientists at BabelColor.com have measured a sample of 30 ColorChecker targets over the years and compared them to published specifications (pre-November 2014 formulations shown, as my unit is older than that)[2]:

babel-cc24
Figure 1.  Average measurements of 30 pre-2015 ColorChecker 24 targets  by BabelColor.com

Note that 1 \Delta E_{00} is supposed to represent a just noticeable color difference, so you can see that with the exception of purple and white the patches do seem to provide a reasonably stable reference.  This information, as well as average patch reflectance from 380 to 730nm at 10nm increments, is available in a spreadsheet from the relative link above.

Bingo!  the L*a*b* color space is a simple transformation away from XYZ_{D50}.  With the data above and a single raw capture of a ColorChecker Passport Photo target in clear noonish sunlight (call it D50) we’ve got 24 as opposed to just 3 sets of raw and reference data triplets needed to solve for the coefficients of the Forward Matrix in Equation 1.  And we can use the larger data set to make sure that the coefficients fit a wider number of potential photographic subjects – that’s why they call it a compromise color matrix.  Shutter Release.

wbdsc_3378
Figure 2. Pre-2015 X-Rite ColorChecker Passport Photo captured by a Nikon D610 mounting a 24-120mm/4 at 120mm, 1/800s, f/5.6, ISO100 in early September at 11:27 on a clear mountain day.  The Correlated Color Temperature is in the 5000-5200K range.

Computing the Coefficients by the Normal Equation

Ok, so now we have the cc24 target illuminated by a roughly D50 source captured in the raw file of a Nikon D610+24-120mm/4. Next we read the mean Raw values in the three color channels for each of the 24 patches with a tool like RawDigger [3].  Then we white balance them based on the third gray patch from the right in the neutral bottom row and demosaic them.  The result is now the white balanced and demosaiced raw data set that we need, one rgb triplet per patch.

We could then obtain the reference XYZ_{D50} corresponding to each patch by measuring them or via ColorChecker values published by X-Rite or BabelColor as above and solve the 72 equations for the coefficients in matrix M that best fit the available data using the Normal Equation as follows:

(3)   \begin{equation*} M = (rgb^T * rgb)^{-1} * rgb * XYZ_{D50} \end{equation*}

with ^T representing the transpose and ^{-1} the inverse of the relative arrays.  Voilà, by the magic of linear algebra Equation 3 will produce the 3×3 compromise color matrix M that will result in the smallest sum of square differences to reference values for the target at hand.  However that turns out not to be the best way to  solve for M because the method gives equal weight to all values, while the Human Visual System tends to be more sensitive in some parts of the XYZ color space than others.

A Better Color Difference Metric: dE2000

In fact in what follows I will use the standard color difference \Delta E_{00} (CIEDE2000) as the minimization criterion.  We can set up a spreadsheet with the 24 rgb triplets and 9 cells representing the coefficients in the forward matrix M as arrays.  Seed the coefficients to, say, 1/3 each.  Matrix multiply the  rgb triplets by the seed matrix, convert to L*a*b*, compute \Delta E_{00} differences to the reference values in Figure 1 above and let Excel Solver figure out what values of the 9 coefficients of M minimize the sum of the differences.  The resulting matrix will be the best compromise found for the 24 patches under that illuminant, taking HVS color sensitivity in consideration therefore giving each patch approximately equal perceptual importance.

The Compromise Forward Matrix

I effectively followed that last procedure but using Matlab/Octave instead of Excel.   Excellent toolkit OptProp and its built-in ColorChecker reference data lent a helping hand[4].  This is the Forward Matrix M obtained for my setup:

Figure 3.

Cool.  These are the \Delta E_{00} differences for each patch resulting from the transformation using OptProp’s ColorChecher reference data.

optprop-database-d2000
Figure 3.  CIEDE2000 difference between the measured raw values from Figure 2 transformed to L*a*b* via XYZD50 by the computed Forward Matrix and the ColorChecker reference values built into optprop that were used to obtain it.

The average \Delta E_{00} is 1.5 and maximum is 3.8 in the light skin patch.  Recall that 1 is a just noticeable difference.  I repeated the exercise this time using the BabelColor 30 database as reference and this is the resulting forward matrix then:

forwardmatrix-d610-babel

Looks like the differences to reference data from BabelColor’s database are more evenly distributed:  average \Delta E_{00} is still 1.5 but Maximum is a lower 3.1 in the ‘green’ patch.

babel30-database-d2000
Figure 4. CIEDE2000 difference between the measured raw values from Figure 2 transformed to L*a*b* via XYZD50 by the computed Forward Matrix and the BabelColor database ColorChecker reference values that were used to obtain it.

To do it properly I should have measured the Spectral Power Distribution of the illuminant and the reflectance of the patches with a spectrophotometer around the time that the capture was taken (didn’t have one, X-Rite ColorMunki Photo/i1 Studio on the way, see post scriptum).  Using published reference SPDs and reflectances pollute the data and contribute to degrade the results – but you get the idea.

Correcting Matrix ‘Errors’: Profiles

Note however that the compromise matrix is just that, a compromise, and even if the procedure had been perfect it could result in relatively large errors.  For better overall color rendering performance the preferred method is to correct such errors through optional nonlinear color profile adjustments that are typically introduced while in XYZ_{D50} via ProPhoto RGB HSV lookup tables.  The tables and application methods are usually referred to as camera ‘profiles’ (ICC and DCP for instance), a subject beyond the scope of this article[5].

Calculating Sensitivity Metamerism Index (SMI)

While the data is out we can calculate the Sensitivity Metamerism Index, which for us is equal to 100 minus 5.5 times the average \Delta E (not \Delta E_{00}) of just the 18 color patches[6].  The shown values (misnamed CRI above) are not necessarily maximized for the given setup because the matrix finding routine is built upon minimizing \Delta E_{00}. Still my D610 with its 24-120mm/4 around D50 shows SMIs of 80 in the first case and 83 in the second.   They jump to 82 and 86 respectively if the routine is setup to minimize \Delta E instead.  Not bad, although I know some people do not give much credence to this metric – and I can understand why.

Step 2: Matrix to Output Color Space

Now that we have done the hard work of determining the Forward Matrix to convert white balanced raw data to the PCS with this illuminant, we need a standard matrix to move it on to the colorimetric color space chosen by the photographer for output.  In this case we will map it to sRGB_{D65} by multiplication with the following transform matrix obtained from Bruce Lindbloom’s site[7]

lindbloom1

The product of this and the forward matrices calculated earlier produces the following combined white balanced raw rgb\rightarrow sRGB linear transforms for a D50ish illuminant.  My results are shown top and bottom, with DXOmark’s for the D610 at D50 in the center for reference[8]:

tosrgbdxo2
Figure 5. Forward Matrices from white balanced raw data to XYZD50 for a Nikon D610+24-120/4  in approximately D50 light. Top: matrix computed through the procedure in this article from the raw data in Figure 2 with optprop’s built in ColorChecker reference values. Middle: DXOmark.com’s D50 matrix. Bottom: matrix computed in this article from the raw data in Figure 2 with BabelColor 30 database ColorChecker reference values.

Pretty close, which confirms the illuminant of my captures was indeed near D50.

So that’s where color matrices come from and why we need them. Next, a closer look and putting them to work.

 

Post Scriptum

I got myself an X-Rite ColorMunki Photo spectrophotometer (it’s identical to the recently rebranded i1 Studio), it’s a blast.  It comes with its own mini ColorChecker 24 Classic, I assume with the new formulation.  I measured it with the Munki and captured it in the raw data in somewhat similar conditions to Figure 2.  Correlated Color Temperature was about 5050K.  These are the matrices that are produced from it (white balanced off my WhiBal card):

Pretty similar and SMI is now 85.  In fact if I set the routine up to minimize \Delta E instead of \Delta E_{00} SMI is 86.  The sum of rows in the XYZ matrix should represent the white point of the target illuminant, which in this case was just the reference data for D50.  The matrix white point is [0.9609  1  0.8214] which has a correlated color temperature of 5016K, according to Bruce.  Good, because D50 has a CCT of 5002K.

Here are the  \Delta E_{00} differences, measured vs captured:

Figure 6. CIEDE2000 obtained from ColorChecker mini measured with a ColorMunki spectrophotometer and captured in the raw data in late december at about 1700m, 2pm, sunny day with snow on the ground.  CCT is about 5050K.

Now average \Delta E_{00} is 1.17, there are only three patches above 2 and the rest seem to be better controlled overall. I guess pushing down the outliers  is the reason why linear matrices are not enough and profiles are necessary.  To get even closer I should have brought the Munki with me and measured the illuminant.  Well, next time.

Notes and References

1. Lots of provisos and simplifications for clarity as always.  I am not a color scientist, so if you spot any mistakes please let me know.
2. The ColorChecker pages at BabelColor.com can be found here. Careful that there was a change in formulations in November 2014.
3. RawDigger can be found here and dcraw here.
4. optprop can be found here.
5. For an example of profile implementation see for instance Adobe’s Digital Negative Specification Version 1.4.0.0.
6. See here for a description of the Sensitivity Metamerism Index and here for DXO’s take on it.
7. See this page at Bruce Lindbloom’s site for precise matrices from XYZD50 to many colorimetric color spaces.
8. DXOmark.com D610 color measurement can be found here by clicking on the ‘color response’ tab.