sRGB↔L*a*b*↔LChab conversions
Posted by Michał ‘mina86’ Nazarewicz on 21st of March 2021 | (cite)
After writing about conversion between sRGB and XYZ colour spaces I’ve been asked about a related process: moving between sRGB and CIELAB (perhaps better known as L*a*b*). As this may be of interest to others, I’ve decided to go ahead and make an article out of it. I’ll also touch on CIELChab which is a closely related colour representation.
The L*a*b* colour space was intended to be perceptually uniform. While it’s not truly uniform it’s nonetheless useful and widely used in the industry. For example, it’s the basis of the ΔE*00 colour difference metric. LChab aim to make L*a*b* easier to interpret by replacing a* and b* axes with more intuitive chroma and hue parameters.
Importantly, the conversion between sRGB and L*a*b* goes through XYZ colour space. As such, the full process has multiple steps with a round trip conversion being: sRGB→XYZ→L*a*b*→XYZ→sRGB. Because of that structure I will describe each of the steps separately.
sRGB to XYZ
For brevity I’ll assume herein that we’re working with 8-bit colour depth, i.e. that red, green and blue components are integers in the
When doing any conversions from sRGB, we first need to recall that it uses gamma correction. Each component of the colour is compressed and needs to be expanded to move to linear RGB space. Only then the conversion matrix
Expressing all that in code, specifically Rust, we get:
/// Converts sRGB colour given as a triple of 8-bit integers /// into XYZ colour space (with white’s Y value equal 1). fn xyz_from_srgb(srgb: [u8; 3]) -> [f32; 3] { // Gamma expansion for 8-bit sRGB components. fn gamma_expansion(value255: u8) -> f32 { return if value255 <= 10 { value255 as f32 / 3294.6 } else { ((value255 as f32 + 14.025) / 269.025).powf(2.4) } } // Convert all components to linear RGB. let r = gamma_expansion(srgb[0]); let g = gamma_expansion(srgb[1]); let b = gamma_expansion(srgb[2]); // Switch to XYZ. [ r*0.4124108464885388 + g*0.3575845678529519 + b*0.18045380393360833, r*0.21264934272065283 + g*0.7151691357059038 + b*0.07218152157344333, r*0.019331758429150258 + g*0.11919485595098397 + b*0.9503900340503373, ] }
As mentioned, when working with different colour depths a different gamma_expansion
function is necessary. It’s not enough to scale all numbers since the
XYZ to L*a*b*
While the sRGB standard defines a specific white point — which is baked into the conversion matrix — CIELAB can be computed in reference to an arbitrary illuminant. Since we’re starting with a colour defined using D65 white point it is what I’m going to use in L*a*b* calculations as well. With that assumption, the formulæ are as follows:
Naïvely one might change lab()
colours. Chromatic adaptation is beyond the scope of this article so for now the important thing to remember is that formulæ and implementations presented here result in L*a*b* colour calculated in reference to D65 illuminant while other sources may use D50 instead.
Speaking of implementation, here it is for the XYZ→CIELAB conversion:
/// Converts colour from XYZ space to CIELAB (a.k.a. L*a*b) using /// D65 reference white point. fn lab_from_xyz(xyz: [f32; 3]) -> [f32; 3] { const EPSILON: f32 = 216.0 / 24389.0; const KAPPA: f32 = 24389.0 / 27.0; fn f(v: f32) -> f32 { if v > EPSILON { v.powf(1.0 / 3.0) } else { (KAPPA * v + 16.0) / 116.0 } } let fx = f(xyz[0] / 0.9504492182750991); let fy = f(xyz[1]); let fz = f(xyz[2] / 1.0889166484304715); [ 116.0 * fy - 16.0, 500.0 * (fx - fy), 200.0 * (fy - fz), ] }
L*a*b* to XYZ
Conversion from L*a*b* to sRGB involves performing all the steps in reverse order. First going from CIELAB to XYZ, then to linear RGB and finally to sRGB. The first step of this process is governed by the following equations:
To optimise calculations it’s also worth to observe that:
Taking that optimisations into consideration, the equations translate into the following Rust function:
/// Converts colour from CIELAB (a.k.a. L*a*b*) space /// using D65 reference white point to XYZ. fn xyz_from_lab(lab: [f32; 3]) -> [f32; 3] { const CBRT_EPSILON: f32 = 6.0 / 29.0; const KAPPA: f32 = 24389.0 / 27.0; fn f_inv(v: f32) -> f32 { if v > CBRT_EPSILON: { v.powi(3) } else { (v * 116.0 - 16.0) / KAPPA }; } let fy = (lab[0] + 16.0) / 116.0; let fx = (lab[1] / 500.0) + fy; let fz = fy - (lab[2] / 200.0); let x = f_inv(fx) * 0.9504492182750991; let y = if lab[0] > 8 { fy.powi(3) } else { lab[0] / KAPPA } let z = f_inv(fz) * 1.0889166484304715; [x, y, z] }
XYZ to sRGB
The second step is moving to sRGB colour space. Similarly to conversion in the other direction, this involves switching from XYZ to linear RGB first (by using the inverse conversion matrix) and then applying gamma compression function to get final sRGB coordinates.
Like before, this assumes 8-bit colour depth output. For general form of the
/// Converts colour given in XYZ colour space (with white’s Y value /// equal 1) into an sRGB colour specified as a triple of 8-bit integers. fn srgb_from_xyz(xyz: [f32; 3]) -> [u8; 3] { // Gamma compression from linear [0, 1] value to 8-bit integer. fn gamma_compression(linear: f32): u8 { let v = if linear <= 0.00313066844250060782371 { 3294.6 * linear } else { 269.025 * linear.pow(5.0 / 12.0) - 14.025 }; v.round().min(255.0).max(0.0) as u8 } let x = xyz[0]; let y = xyz[1]; let z = xyz[2]; let r = x* 3.240812398895283 - y*1.5373084456298136 - z*0.4985865229069666; let g = x*-0.9692430170086407 + y*1.8759663029085742 + z*0.04155503085668564; let b = x* 0.055638398436112804 - y*0.20400746093241362 + z*1.0571295702861434; [ gamma_compression(r), gamma_compression(g), gamma_compression(b), ] }
And that’s it as far as converting between L*a*b* and sRGB colour spaces is concerned. With the presented code, lab_from_xyz(xyz_from_srgb([r, g, b]))
gets an sRGB colour converted into L*a*b* and srgb_from_xyz(xyz_from_lab([l, a, b]))
moves in the opposite direction. There’s still a closely related colour representation — CIELChab or LChab — which we’ll look at next.
L*a*b* and LChab
The a* and b* components of the L*a*b* colour are not easy to interpret. Even knowing that they component represents green-red and blue-yellow axes respectively, it’s tricky to visualise the colour space. To make it easier for humans to work with, CIE defined an LChab colour space as a conversion of L*a*b* into a cylindrical coordinate system. It replaces a* and b* components with more intuitive chroma and hue. The relations between all the components is as follows:
Depending on use case it may be better to keep hue expressed in radians or to convert it to degrees. The former is usually easier for computation since mathematical instructions typically operate in radians. On the other hand, the latter is more intuitive as it’s easier to think in degrees than radians.
Below are the formulæ expressed in code with conversion to degrees present:
/// Converts colour given in CIELAB (a.k.a. L*a*b) space to /// CIELCh(ab) with hue expressed in degrees. fn lch_from_lab(lab: [f32; 3]) -> [f32; 3] { [ lab[0], lab[1].hypot(lab[2]), lab[2].atan2(lab[1]) * 360.0 / std::f32::consts::TAU, ] } /// Converts colour given in CIELCh(ab) space with hue /// expressed in degrees to CIELAB (a.k.a. L*a*b*). fn lab_from_lch(lch: [f32; 3]) -> [f32; 3] { let h_rad = lch[2] * std::f32::consts::TAU / 360.0; [ lch[0], lch[1] * h_rad.cos(), lch[1] * h_rad.sin(), ] }
Accuracy benchmark and possible improvements
To finish up, let’s test how good the presented code is. It’s not easy to benchmark precision of the calculations in the conversion code. There is no authoritative set of colours with their coordinates given in RGB and L*a*b* colour spaces. Still, there is one thing that’s easy to check.
Colours with a* and b* components equal to zero represent shades of grey. In other words, for any
fn main() { let mut count: usize = 0; let mut error: f64 = 0.0; let mut max: f32 = 0.0; for i in 0..256 { let lab = lab_from_xyz(xyz_from_srgb([i as u8, i as u8, i as u8])); if lab[1] != 0.0 || lab[2] != 0.0 { eprintln!("{:3} {:20} {:20} {:20}", i, lab[0], lab[1], lab[2]); count += 1; error = (lab[1] as f64).mul_add(lab[1] as f64, error); error = (lab[2] as f64).mul_add(lab[2] as f64, error); max = max.max(lab[1].abs()); max = max.max(lab[2].abs()); } } println!("{:3} {}e-9 {}e-6", count, error * 1e9, max * 1e6); }
When executed the code shows that 48 colours weren’t converted correctly with the maximum error of less than 30 millionth and total square error of less than 21 billionth. With floating point numbers some imprecision is always to be expected so this isn’t a bad result by no means. However, on architectures with fused multiply-add instruction accuracy can be improved by improving the matrix multiplications at the end of xyz_from_srgb
and at the beginning of srgb_from_xyz
function. I’ll leave that as an exercise to the reader.