Color: CIELab and Tristimulus XYZ

Let’s talk about CIELab. I’m not going to provide any background or discussion. This wiki article is an excellent introduction.

And there’s a very nice article about HSB. Furthermore, it includes an example with many more parameters than just hue, saturation, and chroma computed, which is outstanding. It also has a lot more information about HSI versus HSB; the differences seem to be more than I was led to believe from my previous reading.

CIELab is specified as a transformation from XYZ tristimulus coordinates. What motivated this post is that the only inverse transformation I’d seen published – until I read the wiki article – is incomplete. That’s a polite way of saying it doesn’t always work correctly. And that’s a polite way of saying it’s wrong.

In fact, the inverse transformation in the wiki post is perfecty straight-forward, and I’ll show it to you in more detail than they do.

So, I want to show you the standard transformation which defines CIELab… and two inverse transformations, one of which is – at best – incomplete. In addition, there is more than one way to write the numbers in the transformations, and I’ll show you what I prefer.

Here’s the definition of CIELab out of Fairchild’s “Color Appearance Models”.

L = 116 f(Y/Yn) - 16

a = 500 ( f(X/Xn) - f(Y/Yn) )

b = 200 (f(Y/Yn) - f(Z/Zn) )

f(x) = x^{1/3}, \text{if }x > .008856

f(x) = 7.787 x + 16/116, \text{if } x \le .008856

Here's the code I wrote for it:

You should notice that my code computes and returns two additional results… \chi\ for "chroma" and h for hue angle (and I chose to use degrees rather than radians). I'll talk about them down the road (but not today).

I should point out, however, that h and \chi\ are nothing more than polar coordinates based on cartesian coordinates a,b.

Oh, properly speaking the CIELab coordinates are denoted with asterisks, L*, a*, b* – but since I’m not going to discuss any other Lab model in this post, I will use L, a, b.

Thanks to the wiki article, I know what the numbers .008856 and 7.787 "really" are: one is the cube of 6/29, and the other is one-third of the square of 29/6:

Oh, and 16/116, in lowest terms, is

As a result, I prefer to write the transformation using exact rational numbers rather than approximate decimal numbers. Here's my current code. It uses "If" to define the function f, and declares f local in addition to the internal variables.

Recall my most recent artists' color wheel…

We need the tristimulus values for a reference white. I chose to measure the values of the white background – boy, that was handy! – and got

w1 = \{96.42,100.003,82.489\}

(We'll see those numbers again in a new post about transforming between RGB and XYZ on my monitor.)

Here are the tristimulus values XYZ for those 12 colors on my monitor, measured by the DigitalColor Meter, which comes with OS X.

measTri = \left(\begin{array}{ccc} 36.151 & 20.801 & 2.655 \\ 38.409 & 24.292 & 3.259 \\ 45.496 & 35.254 & 5.157 \\ 59.079 & 56.259 & 8.792 \\ 80.151 & 88.846 & 14.438 \\ 51.678 & 72.464 & 12.347 \\ 44. & 68.045 & 11.783 \\ 60.269 & 79.199 & 79.834 \\ 16.269 & 11.154 & 68.051 \\ 23.947 & 15.573 & 68.616 \\ 52.42 & 31.955 & 70.706 \\ 39.606 & 23.169 & 17.108\end{array}\right)

Let me apply the XYZtoLab function to those 12 XYZ sets. Call the result "fullLab"; in addition to the L,a,b values, it has returned \chi\ and h, chroma and hue angle. It's been a while since I displayed matrices; the "TeXForm" command gives me what I need for the LaTeX.

fullLab = \left(\begin{array}{ccccc} 52.7303 & 64.2894 & 54.8817 & 84.5288 & 40.4863 \\ 56.3782 & 55.921 & 56.6725 & 79.6174 & 45.3824 \\ 65.9451 & 36.0473 & 61.9073 & 71.6374 & 59.7887 \\ 79.76 & 11.9181 & 70.2776 & 71.281 & 80.375 \\ 95.5149 & -10.5385 & 80.3919 & 81.0797 & 97.4682 \\ 88.1904 & -42.95 & 73.4483 & 85.0844 & 120.318 \\ 86.0279 & -54.8307 & 71.362 & 89.9941 & 127.537 \\ 91.3232 & -35.0906 & -12.7908 & 37.3491 & -159.973 \\ 39.8381 & 35.6103 & -91.303 & 98.0017 & -68.693 \\ 46.4086 & 45.288 & -80.4922 & 92.358 & -60.6362 \\ 63.3049 & 66.2499 & -53.2513 & 84.9985 & -38.7922 \\ 55.2453 & 64.5872 & 4.4514 & 64.7404 & 3.94263\end{array}\right)

Let me cut those results down so we just have L, a, b. The following line says "drop none of the rows, and drop the last two columns" – and then give me TeX after displaying the matrix…

compLab = \left(\begin{array}{ccc} 52.7303 & 64.2894 & 54.8817 \\ 56.3782 & 55.921 & 56.6725 \\ 65.9451 & 36.0473 & 61.9073 \\ 79.76 & 11.9181 & 70.2776 \\ 95.5149 & -10.5385 & 80.3919 \\ 88.1904 & -42.95 & 73.4483 \\ 86.0279 & -54.8307 & 71.362 \\ 91.3232 & -35.0906 & -12.7908 \\ 39.8381 & 35.6103 & -91.303 \\ 46.4086 & 45.288 & -80.4922 \\ 63.3049 & 66.2499 & -53.2513 \\ 55.2453 & 64.5872 & 4.4514\end{array}\right)

Now, let me use the DigitalColor Meter to measure the Lab values directly:

measLab = \left(\begin{array}{ccc} 52.734 & 63.973 & 54.605 \\ 56.382 & 55.633 & 56.387 \\ 65.943 & 35.855 & 61.602 \\ 79.76 & 11.812 & 69.93 \\ 95.518 & -10.504 & 80.02 \\ 88.194 & -42.789 & 73.102 \\ 86.032 & -54.625 & 71.023 \\ 91.327 & -34.961 & -12.742 \\ 39.84 & 35.406 & -90.949 \\ 46.407 & 45.047 & -80.184 \\ 63.308 & 65.926 & -53.043 \\ 55.248 & 64.266 & 4.398\end{array}\right)

Look at the differences between computed Lab values and measured Lab values:

\left(\begin{array}{ccc} -0.00367709 & 0.316393 & 0.276656 \\ -0.0037513 & 0.287984 & 0.285469 \\ 0.00212331 & 0.192295 & 0.305267 \\ 0.0000393356 & 0.106144 & 0.347565 \\ -0.00308963 & -0.0344622 & 0.371929 \\ -0.0035899 & -0.160959 & 0.346269 \\ -0.00407541 & -0.205681 & 0.339003 \\ -0.00381671 & -0.129572 & -0.0488465 \\ -0.0019136 & 0.204342 & -0.353955 \\ 0.0015953 & 0.240992 & -0.308182 \\ -0.00311854 & 0.323858 & -0.20832 \\ -0.00272351 & 0.321233 & 0.053397\end{array}\right)

Hmm. The differences in a and b are two orders of magnitude greater than the differences in L. I don't like that.

Here are the % differences:

\left(\begin{array}{ccc} -0.00697289 & 0.494573 & 0.50665 \\ -0.00665336 & 0.51765 & 0.506267 \\ 0.00321991 & 0.536314 & 0.495547 \\ 0.0000493174 & 0.898614 & 0.497018 \\ -0.00323461 & 0.328087 & 0.464795 \\ -0.00407045 & 0.376168 & 0.473679 \\ -0.00473709 & 0.376532 & 0.477315 \\ -0.00417917 & 0.370618 & 0.38335 \\ -0.0048032 & 0.577139 & 0.389179 \\ 0.00343763 & 0.534978 & 0.384344 \\ -0.00492598 & 0.491245 & 0.392737 \\ -0.00492962 & 0.499849 & 1.21412\end{array}\right)

Yikes! The discrepancy in the last measured versus computed b value is more than 1%.

I don't know why. (And I’ve checked the measured tristimulus and Lab values.) I'm hoping it has to do with our having only 8 bits for each of R,G,B — which is how the colors are specified physically.

In other words, I'm not going to worry about the discrepancy between computed and measured. Not until and unless I have reason to. The code is simple enough, and I've checked the numerical values used in it. My computed Lab values should be correct, given the measured tristimulus values.

(I’m confident that I’m right, and it will be embarrassing if I find I’m wrong somehow.)

Now let's look at the inverse, going from Lab to XYZ. Kang ("Computational Color Technology") offers the following – and I’ve seen that Glassner (bibliography) does, too… which was easy enough to code up:

There's just one little problem: the inversion is conditional only on the value of L. Equivalently, this inversion is conditional only on the value of Y/Yn; this inverse does not depend on X/Xn or Z/Zn.

But they affected the forward computation; how can they not affect the reverse?

Was Fairchild wrong?

I don't think so. I don't have the CIE documents which define CIELab — but Wyszecki and Stiles are quite explicit: if the value of X/Xn is less than that cutoff, use the linear form of f, rather than the cube root, to compute f(X/Xn); and if the value of Z/Zn is less than that cutoff, use the linear form of f, rather than the cube root, to compute f(Z/Zn). And ditto for Y/Yn, of course.

That's what Fairchild wrote, in different language.

Kang's inversion is correct if all three of X/Xn, Y/Yn, and Z/Zn are on the same side of the cutoff – that is, if the same form of f was used for all three ratios. But the definition of CIELab  – as far as I can tell from secondary sources – can require us to use the cube root three times, twice, once, or not at all – so Kang's inversion works if we used the cube root three times or never, but it must fail to invert correctly if we used the cube root only once or twice.

Here's an example. The first line is a set of XYZ… the second line is X/Xn, Y/Yn, Z/Zn, and we see that X/Xn is less than .008… the third line converts to CIELab and drops the last two answers (hue and chroma)… the fourth line converts the Lab coordinates back to XYZ – incorrectly: X/Xn is not right.

We did not recover the original X, which used the linear function instead of the cube root: we got .54 instead of .5 . We recovered Y and Z just fine, because the forward calculation used the cube roots for both of them — so we see that the inversion used the cube on all three, rather than just on two.

What to do? I was simply going to use the inversion, but check it. If I found that I needed to fix it up, I would.

But it turns out that the specification is almost trivial to invert correctly. Let's look at f(x), for extremely small values for x.

Here's the function from the XYZ to CIELab code… and evaluate it at the break point (6/29)^3:

That value of f should not be surprising: we're taking the cube root of (6/29)^3. Here's a picture, with the two pieces in different colors.

Oh, if you've never seen the "With" command, here's a nice example of what it does for us: the value "a" occurs in the code 5 times, but I only had to write (6/29)^3 once. In this case, I was just saving myself having to repeat – and possibly mess up – 5 identical constants. Another excellent use is for parametrizing some aspect of a plot, say something that is drawn from -a to a, when you're not sure what would be a good value of a.

Perhaps I should also remind you that I am using Dave Park's Presentations package: "Draw" a function rather than "Plot" a function; it makes it very easy to combine things.

That f(x) is defined piecewise is not a problem: we simply invert the two pieces of the function separately. if f(x) = y < 6/29, we invert the linear portion; otherwise, we invert the cube root.

Oh, I'm jumping ahead. Let's see why we only need the inverse of f.

The equation for L is

L = 116 f(Y/Yn) - 16

so

f(Y/Yn) = \frac{L+16}{116}

and then we apply f^{-1}\ to both sides,

f^{-1}(f(Y/Yn)) = Y/Yn = f^{-1} (\frac{L+16}{116})\ ,

which gave us Y/Yn on the LHS.

Finally,

Y = Yn f^{-1} (\frac{L+16}{116})\ .

The reason f^{-1}\ will suffice is that this first equation only involved one term in f.

Similarly,

a = 500 ( f(X/Xn) - f(Y/Yn) ) = 500 f(X/Xn) - 500 f(Y/Yn)

500 f(X/Xn) = a + + 500 f(Y/Yn)

f(X/Xn) = a/500 + f(Y/Yn)

but we have f(Y/Yn) from the equation for L, so we can eliminate it (this is crucial):

f(X/Xn) = a/500 + \frac{L+16}{116}\ .

Now apply f^{-1}\ and get

f^{-1}(f(X/Xn)) = X/Xn = f^{-1} (a/500 + \frac{L+16}{116})\ .

Getting Z/Zn and then Z is just like getting X/Xn; the third equation inverts as

Z = Zn f^{-1} (\frac{L+16}{116} - b/200)\ .

Now we just need to work out f^{-1}\ .

To invert the cube root is trivial: cube the argument. To invert the linear portion…

That is equivalent to what the wiki article has, namely

Here's my code for the inversion. I should add two optional parameters, so I can apply the code to a list that includes hue and chroma, but I'm not going to yet.

Let's run the inversion on the computed Lab values.

revLab = \left(\begin{array}{ccc} 36.151 & 20.801 & 2.655 \\ 38.409 & 24.292 & 3.259 \\ 45.496 & 35.254 & 5.157 \\ 59.079 & 56.259 & 8.792 \\ 80.151 & 88.846 & 14.438 \\ 51.678 & 72.464 & 12.347 \\ 44. & 68.045 & 11.783 \\ 60.269 & 79.199 & 79.834 \\ 16.269 & 11.154 & 68.051 \\ 23.947 & 15.573 & 68.616 \\ 52.42 & 31.955 & 70.706 \\ 39.606 & 23.169 & 17.108\end{array}\right)

Are those the XYZ tristimulus values we started with? Let's compute the differences between those and the original XYZ values:

I see a lot of very small numbers… are they all small?

Yes.

So. I am comfortable with displaying color combinations from books, specified in CMYK. Unfortunately, we learned that while Mathematica can draw a CMYK spec correctly on my monitor, it cannot tell me what the RGB is… I have to measure the RGB. Many of my books also specify an RGB as well as CMYK, but I don't trust their translation; I prefer to let Mathematica use my monitor profile to get RGB for my device.

Once I have the RGB, I can convert to HSB using my code – which was written to confirm that I understood the algorithm – or I can let Mathematica do it. I am comfortable with HSB for choosing colors on my computer.

Now I've added CIELab to the mix. It moves between XYZ tristimulus values and Lab. I can always measure the tristimulus values – I did that for my artists' color wheel.

But it would be nice if I could move between RGB and XYZ. And I used to be able to do that. As I said in a recent diary post, my monitor profile has been changed… I have a copy of the original, and I can see the changes. I found them, of course, when I discovered that my current measured tristimulus XYZ do not agree with my old ones.

No problem. I know what the changes are and I can write new transformations to move between RGB and XYZ – and that in turn will let me move between RGB and CIELab, and between HSB and CIELab.

Posted in color. Tags: , . 4 Comments »

4 Responses to “Color: CIELab and Tristimulus XYZ”

  1. Bob Says:

    Excellent!

  2. Martin Says:

    Your inverse transformation is not quite correct, as far as I can see. For LabToXYZ you wrote:
    fi[y_] = If [y > 6 / 29, y^3, 3 (6 / 29)^2 (y – 4/29)];
    The condition on the Wikipedia article is different: It should be “y > (6 / 29)^3”. And not “y> 6/29”, as you wrote.
    See also: https://en.wikipedia.org/wiki/Lab_color_space

  3. Martin Says:

    I want to add: very nice blog you have here! Especially the CIE stuff.

  4. Martin Says:

    Hi again.
    Just realized, that your calculations actually are right and Wikipedia is wrong!
    So please disregard my previous comment about your transformation not beeing correct.


Leave a comment