Don’t forget @pippin the master of babl and gegl.

# Did I correctly implement the math for converting from CIE L*a*b* to sRGB (or other RGB colorspaces) in this shader?

**Tynach**#5

I did check Bruce Lindbloom’s site, and his math is a *mess*. He seems stuck in the thought process of treating each channel as independent of the others, which is a huge mistake to do inside GLSL (and causes duplicate code).

While it was nice to see some proof that I did the math regarding the L* trc correctly (his `κ`

and `ϵ`

variables match my `slope`

and `cutoffToGamma`

variables, except his `κ`

is multiplied by 100), the fact that he implements it differently for each channel (unlike Wikipedia, which has its own issues that caused me to have to do my own math to begin with for calculating those variables’ values) caused a ton of issues for me early on.

On top of that, he has an extra step in his RGB↔XYZ matrix calculation math that is unnecessary and can simply be left out (dividing the chromaticity coordinates for each primary by their ‘y’ component). I’m seriously glad I learned that math from Ryan Juckett’s page on the subject instead of from Bruce Lindbloom’s page, since it gives a thorough and detailed explanation and reasoning for each step and doesn’t result in any unnecessary steps.

I have no doubt that Bruce Lindbloom is very intelligent and probably understands all of this more than most, but his pages are severely lacking in detail and explanation - and in at least one case that’s probably led to some unnecessary work for some people. After all, the math he shows isn’t *wrong*, it’s just got an unnecessary step in it.

All that said, part of my uncertainty over L*a*b* conversion is because I don’t know what is and isn’t unnecessary duplicate work in this context. It almost looks like the ‘+16’ thing is meant to be part of the standard L* trc, if one were doing it in the 0-100 range instead of 0-1. Since my code expects to apply a trc on a 0-1 scale, it adds 0.16 as the ‘offset’ parameter, so might be taken care of there.

But it’s really hard to determine if that’s the case or not, since Wikipedia also has ‘+16’ bits in various parts of the per-channel equations, even though it (like my code) also abstracts the trc part as a single equation to apply to all 3 channels identically. So I still feel like I might not be fully understanding the math, and might be doing it incorrectly.

**Elle**(Elle Stone) #6

I don’t have a clue what most of the above posts in this thread are actually talking about. My only programming language is “c” and I’m not all that good at reading or writing “c” code. But Lindbloom’s equations are correct. I used his equations to rewrite babl’s XYZ<->LAB equations quite awhile ago, using double floating point. See this file:

And specifically look at these two functions:

static inline void XYZ_to_LAB

static inline void LAB_to_XYZ

The same file has some floating point functions that are more efficiently though less transparently written.

The babl CIE.c code produces the same LAB values as ArgyllCMS xicclu utility, so I’m confident that it’s correct. But it’s correct for ICC profiles which are not the same as the color spaces from which they were made because ICC profiles are chromatically adapted to D50 from whatever white point is given in the color space specs (D65 for sRGB). Please note that Lindbloom’s equations are correct but his values for D65 and D50 are not - he uses ASTM values instead of the values given in the sRGB specs and the ICC specs.

**Jossie**(Hermann-Josef) #7

For the “16 thing” you might want to consult Charles Poynton’s book. He writes “the scale factor of 116 and the offset of -16 modify the pure power law function such that the best-fit pure power function has an exponent of 0.42, not 1/3.” if that is what you are talking about.

Hermann-Josef

**Tynach**#8

Ok, I think I should probably explain some of my code, since GLSL is weird in a few ways, but also since I’m doing things in ways that seem significantly differently than other code (such as your own) and formulas suggest. I’ve been studying the code you linked to, especially `XYZ_to_LAB`

as it applies the gamma/transfer curve identically to all channels (and apparently it’s like that in Bruce Lindbloom’s XYZ→L*a*b* page as well; wish I’d looked at it and not just the page for the other direction). It looks like we might be doing the same thing, but I’m still unsure, so I’ll try converting what my code does into something more similar.

Firstly, GLSL runs on GPUs, and as such, it’s apparently considered bad practice to have code that branches much. So within my `toGamma`

and `toLinear`

functions, instead of `if`

and ternary statements, I have it simply compute both ‘branches’ and then ‘mix’ the results based on a per-channel conditional. So, lets convert this GLSL code to C:

```
vec4 toGamma(vec4 color, const transfer trc)
{
bvec4 cutoff = lessThan(color, vec4(trc.cutoffToGamma));
vec4 higher = (1.0 + trc.off)*pow(color, vec4(1.0/trc.power)) - trc.off;
vec4 lower = color*trc.slope;
color = mix(higher, lower, cutoff);
return color;
}
```

In C, the definition of the `transfer`

struct (very easy to find as it’s the first thing I define in my shaders) would remain the same, except I might change `float`

to `double`

(in GLSL, doubles *are* floats, you just change the precision - which is literally the first few lines of code in my shader).

So when we rewrite `toGamma`

in such a way that the C version is similar to `XYZ_to_LAB`

, we’d end up with something like this:

```
void toGamma(double r, double g, double b, transfer trc, double *to_r, double *to_g, double *to_b)
{
*to_r = r >= trc.cutoffToGamma ? (1.0 + trc.off) * pow(r, 1.0 / trc.power) - trc.off : r * trc.slope;
*to_g = g >= trc.cutoffToGamma ? (1.0 + trc.off) * pow(g, 1.0 / trc.power) - trc.off : g * trc.slope;
*to_b = b >= trc.cutoffToGamma ? (1.0 + trc.off) * pow(b, 1.0 / trc.power) - trc.off : b * trc.slope;
}
```

So, prior to receiving these responses, I had made a ‘private’ shader on Shadertoy that applied my L*a*b*→XYZ code backwards. I first convert each pixel in an image from RGB to XYZ, then from there I worked through my code backwards to go from XYZ to L*a*b*. After all that, I slapped in my original L*a*b*→XYZ code and checked to see if it gave me the original image back.

I admit that at first it didn’t, but that was because I’d swapped around the order that 2 things were in (or rather, I kept the original order instead of having them occur in the opposite order). After fixing that, I indeed got the original image back. After that, the darkest, closest to black pixels were for some reason labeled as ‘out of gamut’ and were thus turned into that gray color from the original shader, but getting rid of that ‘out of gamut’ check fixed that.

So, since it’s easier to compare XYZ→L*a*b*, ever since your response I’ve been working to compare your version with mine. Here’s the basic GLSL code for just the XYZ→L*a*b* conversion (originally just inlined inside `mainImage`

, but here I’ve split it into a separate function and included the definitions (and their comments) of a few things for convenience):

```
// Parameters for transfer characteristics (gamma curves)
struct transfer {
// Exponent used to linearize the signal
float power;
// Offset from 0.0 for the exponential curve
float off;
// Slope of linear segment near 0
float slope;
// Values below this are divided by slope during linearization
float cutoffToLinear;
// Values below this are multiplied by slope during gamma correction
float cutoffToGamma;
};
// Creates a whitepoint's xyz chromaticity coordinates from the given xy coordinates
#define white(x, y)\
vec3(x, y, (1.0 - x - y))
#define Bright(w)\
((w)/w.y)
// Standard illuminant D50 according to ICC standards (they specify a hex value
// for the 16-bit integer representation, as well as a specific way to decode it,
// so I did some math to figure out exactly the values they expect)
const vec3 whiteD50I = white(31595.0/91393.0, 32768.0/91393.0);
// Gamma for the CIE L*a*b* Lightness scale
const transfer gamLab = transfer(3.0, 0.16, 243.89/27.0, 0.08, 216.0/24389.0);
// Color Appearance Model (or 'reference') white point
const vec3 whiteCam = whiteD50I;
vec4 XYZ_to_LAB(vec4 color)
{
color.xyz /= Bright(whiteCam);
color = toGamma(color, gamLab);
color.xyz = (color.xyz*100.0 + 16.0)/116.0;
color.xz = (color.xz - color.y)*vec2(500, -200);
color.y = color.y*116.0 - 16.0;
return color;
}
```

If I convert it to C (like earlier), I get something like this:

```
#define D50_WHITE_REF_x (31595.0 / 91393.0)
#define D50_WHITE_REF_y (32768.0 / 91393.0)
#define D50_WHITE_REF_X (D50_WHITE_REF_x / D50_WHITE_REF_y)
#define D50_WHITE_REF_Y (1.0)
#define D50_WHITE_REF_Z ((1.0 - D50_WHITE_REF_x - D50_WHITE_REF_y) / D50_WHITE_REF_y)
#define LAB_GAMMA (3.0)
#define LAB_OFFSET (0.16)
#define LAB_SLOPE (243.89 / 27.0)
#define LAB_CUT_GAMMA (216.0 / 24389.0)
void XYZ_to_LAB(double X, double Y, double Z, double *to_L, double *to_a, double *to_b)
{
double xr = X / D50_WHITE_REF_X;
double yr = Y / D50_WHITE_REF_Y;
double zr = Z / D50_WHITE_REF_Z;
// Doing this inline, but I could instead pack LAB_* into a 'transfer' struct and reuse 'to_gamma()'
double fx = xr >= LAB_CUT_GAMMA ? (1.0 + LAB_OFFSET) * pow(xr, 1.0 / LAB_GAMMA) - LAB_OFFSET : xr * LAB_SLOPE;
double fy = yr >= LAB_CUT_GAMMA ? (1.0 + LAB_OFFSET) * pow(yr, 1.0 / LAB_GAMMA) - LAB_OFFSET : yr * LAB_SLOPE;
double fz = zr >= LAB_CUT_GAMMA ? (1.0 + LAB_OFFSET) * pow(zr, 1.0 / LAB_GAMMA) - LAB_OFFSET : zr * LAB_SLOPE;
fx = (fx * 100.0 + 16.0) / 116.0;
fy = (fy * 100.0 + 16.0) / 116.0;
fz = (fz * 100.0 + 16.0) / 116.0;
// I could instead set '*to_L' before modifying 'fy' above by simply multiplying by 100
*to_L = fy * 116.0 - 16.0;
*to_a = (fx - fy) * 500.0;
*to_b = (fz - fy) * -200.0;
}
```

I think it’s safe to say that your code has fewer instructions than mine, so would run *much* faster! However, I’d have to write a special ‘gamma’ function specifically for L*a*b*, rather than reusing the one I use for sRGB, Rec. 709, etc. At this point I suppose it depends on how much faster it is overall; I know mobile devices choke on the 3D version of my gamut shader, so if it ends up being enough of a speed boost to raise my phone’s framerate up from 15fps (which is *after* setting `steps`

to 30, down from 120), that’d be a pretty good win.

I feel it’s worth pointing out that my `toGamma`

function takes values in the same range as it outputs (0 to 1), and so the offset and multiplier are scaled accordingly (to `0.16`

and `1.16`

(calculated as `1.0 + trc.off`

)). Also, I could shorten the gamma/transfer curve bit by just using the cube root function (`cbrt()`

, which your code uses and I had to look up because I wasn’t familiar with it) instead of taking the value to the 1/3rd power. However, while I did inline that part of the code instead of reusing it, I decided to keep its content intact. Shows what the other code would likely turn into if inlined, at the very least.

So, it seems as if my code does at least take the +16 stuff into account when using that function, but now to figure out if it does so in a way that leads to correct results!

**Note**: I have not tested any of the C code I wrote above. I’m assuming it’s correct as I basically performed a literal translation from GLSL to C (roughly in the style your code’s in), but I may have changed something while editing it and caused the rest of it to break due to forgetting to change all the other things to match. I’ll test things more properly after I go eat lunch.

**Tynach**#9

Yeah, in my code that’s taken care of by the `toGamma`

and `toLinear`

functions, with an offset of 0.16 and a multiplier calculated from that (which ends up being `1.16`

). I give more details about that in my response to Elle, along with a translation of what my GLSL code does to C.

I sadly don’t have a copy of that book, and don’t currently have any source of income, but it’s really neat to hear that it gives all the explanations and details for things like that! I have read Poynton’s Color and Gamma FAQs on his website, though. I’ll keep my eyes out for it, see if it’s on sale or anything.

**aurelienpierre**(Aurélien Pierre) #10

Lindbloom has errors in the functions he uses for the XYZ -> Lab (or the other way around, I don’t remember) conversions. When reviewing darktable color spaces conversions lib, at first I thought the lib was wrong, but after checking other sources, turns out there are some mistakes in Lindbloom’s website.

**Tynach**#11

It’s also easy to simply misinterpret his equations. There’s a reason why I described them as ‘a mess’. For example, he defines everything backwards - constants are usually at the bottom, variables directly using those constants are above them, and variables using those variables are above that.

At some point I think there was some mixed around and not strictly backwards, but I really can’t remember for sure. I have trouble with spatial relations *any*way, so it could be I *thought* they were not in strictly reverse order even though they were.

**Elle**(Elle Stone) #12

Can you give some specifics? Even an example set of, for example, sRGB RGB channel values that don’t convert properly using the Lindbloom equations? Just saying “oh, there are errors in Lindbloom’s equations” doesn’t really help, doesn’t say anything at all to be honest.

**Elle**(Elle Stone) #13

@Tynach - you might try the colour-science python code to check accuracy of your code. Example commands are found scattered through this thread:

**gwgill**(Graeme W. Gill) #14

There are several mathematically equivalent ways of doing XYZ to L*a*b*, some more convoluted than others!

For XYZ to L*a*b*, the non-linear x’ y’ z’ values are offset by 16 (i.e. range from 16/116 to 116/116), which is then scaled and shifted to the range 0 … 100 for the L* value, but not for the computation of the a* and b* values. Note that this offset is also implicit in result of the power curve applied to values above 0.008856451586 - i.e. the power curve is shifted to make it meet the linear segment used for small values.

From argyllcms/icc/icc.c icmXYZ2Lab():

```
void icmXYZ2Lab(icmXYZNumber *w, double *out, double *in) {
double X = in[0], Y = in[1], Z = in[2];
double x,y,z,fx,fy,fz;
x = X/w->X;
y = Y/w->Y;
z = Z/w->Z;
if (x > 0.008856451586)
fx = pow(x,1.0/3.0);
else
fx = 7.787036979 * x + 16.0/116.0;
if (y > 0.008856451586)
fy = pow(y,1.0/3.0);
else
fy = 7.787036979 * y + 16.0/116.0;
if (z > 0.008856451586)
fz = pow(z,1.0/3.0);
else
fz = 7.787036979 * z + 16.0/116.0;
out[0] = 116.0 * fy - 16.0;
out[1] = 500.0 * (fx - fy);
out[2] = 200.0 * (fy - fz);
}
```

**afre**#15

Since we are already on a pinging spree, I will also ping @jdc in case he would like to comment.

**Tynach**#16

Took me a bit of time to learn the API, but using this script (where I manually specify everything except the Bradford transformation - which I don’t specify because I checked to ensure it was identical to the transformation matrix I use in my code), I managed to mostly confirm my XYZ→L*a*b* code:

```
import colour
whiteD65S = [0.3127, 0.329]
whiteD50I = [31595.0/91393.0, 32768.0/91393.0]
toXyz = [[506752.0/1228815.0, 87881.0/245763.0, 12673.0/70218.0],
[87098.0/409605.0, 175762.0/245763.0, 12673.0/175545.0],
[7918.0/409605.0, 87881.0/737289.0, 1001167.0/1053270.0]]
rgb = [0, 0, 1]
xyz = colour.RGB_to_XYZ(rgb, whiteD65S, whiteD50I, toXyz, 'Bradford')
lab = colour.XYZ_to_Lab(xyz, illuminant=whiteD50I)
print("XYZ:")
print("X: {:.15f}".format(xyz[0]))
print("Y: {:.15f}".format(xyz[1]))
print("Z: {:.15f}".format(xyz[2]))
print()
print("L*a*b*:")
print("L*: {:.15f}".format(lab[0]))
print("a*: {:.15f}".format(lab[1]))
print("b*: {:.15f}".format(lab[2]))
```

This gives the following output:

```
XYZ:
X: 0.143044519701569
Y: 0.060609958905422
Z: 0.713903816952745
L*a*b*:
L*: 29.565833519650816
a*: 68.285587891529943
b*: -112.033082056164830
```

Feeding my C code from the other post some XYZ values that were precisely calculated as fractions (using Qalculate, which also was used to calculate the fractional form for the RGB→XYZ matrix used in the Python code), I got this output:

```
XYZ:
X: 0.143044519701569
Y: 0.060609958905422
Z: 0.713903816952745
L*a*b*:
L*: 29.565833519650823
a*: 68.285587891529858
b*: -112.033082056164801
```

Seems quite a few significant digits survived! Had to specify so much information into Colour-Science though… Otherwise each value along the way would be slightly different than expected, and by the time it got to L*a*b* things were off after the first digit following the decimal place.

If I used the same RGB→XYZ matrix used by most (rather than calculating my own from the chromaticities of the primaries/white point), it probably would have turned out the same from the start.

**anon11264400**#17

Not sure what you are doing exactly with the fractional nuttery here, but Colour has plenty of defines to get the values you are manually putting in. Also note that there are two methods to calculate the NPMs for RGB to XYZ, and Colour provides both specification and derived as required.

You’ve wasted a good deal of effort manually entering values that Colour already provides with tremendous accuracy.

**Tynach**#18

Well, it seems my L*a*b* code seems to agree with at least Colour-Science (not so sure about xicclu, which gives yet *different* results compared to either, even when Colour-Science is *not* using my pre-calculated matrix).

But, it’s interesting that the offset is applied to L*, but not a* and b*. That might explain a lot, and overall it’s likely that maybe some of the extra work I’m doing in my own code is adjusting for that or something. Or might not be, I really am not sure.

But I’m still somewhat confused - your post states, “For L*a*b* to XYZ, …”, but your code is for the other way around. And that paragraph of your post starts with talking about x’, y’, and z’, *then* goes into L*, a*, and b*.

Your code snippet also seems to apply the offset in the case that they are below the 0.008… value, but *not* to the gamma-corrected value - ~~which I think might cause continuity issues if it’s not corrected for elsewhere~~…

… Wait a minute, that’s not the normal ‘slope’ value. That’s… A close approximation for the slope of the cubic root curve at `216.0/24389.0`

? If I’m doing my math right (which I probably am not; I’m just messing around with offsets and division instead of dealing with derivatives), a more accurate slope would be `841.0/108.0`

. That gives `7.78703703703703...`

I’ve not seen that figure given in any of the equation variants, but it’s interesting nonetheless. It does remove the discontinuity between the power function and the linear fragment of the line, though, so scratch that part of this post out up there.

**Tynach**#19

Yeah, I saw some mentions in the API docs about that, but couldn’t figure out how to use that part of the API. Meanwhile I already knew how to calculate the matrices in Qalculate, so I just did that.

If you could show me how to compute them in Colour, or point me to the proper part of the documentation (which I probably accidentally skipped right over without realizing it), that’d be really helpful! The library is so far pretty fantastic, I just don’t really know how to use it!

In particular, I could not find where the chromaticities for different white points (like D50 or D65) were kept - only their spectral power distributions (which might result in much more precise chromaticity coordinates, but also would likely result in inaccuracies when comparing to how standards like sRGB and ICC v4 define white points like D65 and D50).

**aurelienpierre**(Aurélien Pierre) #20

here : http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html

the y_r equation is wrong, it’s supposed to be the same as x_r and z_r:

**Tynach**#21

I think it might actually be the same, now that I’m looking closer at it. `(L+16)/116`

is used both as the term to cube, *and* is the definition for `fy`

- and on top of that, `L/κ`

is equivalent to `(116*fy - 16)/κ`

since you’d literally just be doing `(116*(L + 16)/116 - 16)/κ`

, causing the two 116s to cancel out, then you’re just adding 16 immediately before subtracting 16.

The fact this wasn’t obvious to me before just shows I’m not really good at reading other people’s equations sometimes. It’s easy for me to get confused with those equations, which is why I didn’t mention Bruce’s website in my original post.

**aurelienpierre**(Aurélien Pierre) #22

I redid the computation on my own 1 month ago, I did not find the same result. Maybe an error on my side, but still… If it’s the same, write it the same, no need to obfuscate the equations. I served my time at the university, I have had my share of sneaky teachers playing this dumb game, I’m too old for this nonsense.

**Tynach**#23

Agreed, especially since I’m doing GLSL, where batching operations into vectors and matrices gives free performance benefits (as there’s dedicated hardware for those operations and the compiler can work faster; GLSL is compiled when it’s run, so fast compile times actually can matter).

Having a slightly more complex route for one of them than is strictly necessary, but as a result ensuring all of them can run simultaneously, can end up being better in the long run. At least in my case.

I’m still probably doing *way* more work than even doing it that way though, since I’m:

- Applying the offset to just L* (making the equivalent of
`fy`

) - Calculating the equivalent to
`fx`

and`fz`

- Removing the offset from all 3 of them (even though it was only indirectly applied to
`fx`

and`fz`

) - Adding the offset back on, sorta, inside
`toLinear()`

But since I’m getting at least 13 digits of accuracy when comparing to Colour-Science’s output (and I’ve since tested with sRGB values of `[0.0001, 0.0001, 0.0001]`

to check against L* values below 0.008, and confirmed accuracy is kept in that case), I’m assuming that my code applies and removes that offset correctly in the end.