In an earlier post, I’ve shown how to calculate an RGB↔XYZ conversion matrix. It’s only natural to follow up with a code for converting between sRGB and XYZ colour spaces. While the matrix is a significant portion of the algorithm, there is one more step necessary: gamma correction.

What is gamma correction?

Human perception of light’s brightness approximates a power function of its intensity. This can be expressed as \(P = S^\alpha\) where \(P\) is the perceived brightness and \(S\) is linear intensity. \(\alpha\) has been experimentally measured to be less than one which means that people are more sensitive to changes to dark colours rather than to bright ones.

Based on that observation, colour space’s encoding can be made more efficient by using higher precision when encoding dark colours and lower when encoding bright ones. This is akin to precision of floating-point numbers scaling with value’s magnitude. In RGB systems, the role of precision scaling is done by gamma correction. When colour is captured (for example from a digital camera) it goes through gamma compression which spaces dark colours apart and packs lighter colours more densely. When displaying an image, the opposite happens and encoded value goes through gamma expansion.

Many RGB systems use a simple \(S = E^\gamma\) expansion formula, where \(E\) is the encoded (or non-linear) value. With decoding \(\gamma\) approximating \(1/\alpha\), equal steps in encoding space correspond roughly to equal steps in perceived brightness. Image on the right demonstrates this by comparing two colour gradients. The first one has been generated by increasing encoded value in equal steps and the second one has been created by doing the same to light intensity. The former includes many dark colours while the latter contains a sudden jump in brightness from black to the next colour.

sRGB uses slightly more complicated formula stitching together two functions: $$ \begin{align} E &= \begin{cases} 12.92 × S & \text{if } S ≤ S_0 \\ 1.055 × S^{1/2.4} - 0.055 & \text{otherwise} \end{cases} \\[.5em] S &= \begin{cases} {E \over 12.92} & \text{if } E ≤ E_0 \\ \left({E + 0.055 \over 1.055}\right)^{2.4} & \text{otherwise} \end{cases} \\[.5em] S_0 &= 0.00313066844250060782371 \\ E_0 &= 12.92 × S_0 \\ &= 0.04044823627710785308233 \end{align} $$

The formulæ assume values are normalised to [0, 1] range. This is not always how they are expressed so a scaling step might be necessary.

sRGB encoding

Most common sRGB encoding uses eight bits per channel which introduces a scaling step: \(E_8 = ⌊E × 255⌉\). In an actual implementation, to increase efficiency and accuracy of gamma operations, it’s best to fuse the multiplication into aforementioned formulæ. With that arguably obvious optimisation, the equations become: $$ \begin{align} E_8 &= \begin{cases} ⌊3294.6 × S⌉ & \text{if } S ≤ S_0 \\ ⌊269.025 × S^{1/2.4} - 14.025⌉ & \text{otherwise} \end{cases} \\[.5em] S &= \begin{cases} {E_8 \over 3294.6} & \text{if } E_8 ≤ 10 \\ \left({E_8 + 14.025 \over 269.025}\right)^{2.4} & \text{otherwise} \end{cases} \\[.5em] S_0 &= 0.00313066844250060782371 \\ \end{align} $$

This isn’t the only way to represent colours of course. For example, 10-bit colour depth changes the scaling factor to 1024; 16-bit high colour uses five bits for red and blue channels while five or six for green producing different scaling factors for different primaries; and HDTV caps the range to [16, 235]. Needless to say, correct formulæ need to be chosen based on the standard in question.

The implementation

And that’s it. Encoding, gamma correction and the conversion matrix are all the necessary pieces to get the conversion implemented. Like before, Rust programmers can take advantage of the srgb crate which implemented full conversion. However, to keep things interesting, in addition, here’s the conversion code written in TypeScript:

type Tripple = [number, number, number];
type Matrix = [Tripple, Tripple, Tripple];

/**
 * A conversion matrix from linear sRGB colour space with coordinates normalised
 * to [0, 1] range into an XYZ space.
 */
const xyzFromRgbMatrix: Matrix = [
	[0.4124108464885388,   0.3575845678529519,  0.18045380393360833],
	[0.21264934272065283,  0.7151691357059038,  0.07218152157344333],
	[0.019331758429150258, 0.11919485595098397, 0.9503900340503373]
];

/**
 * A conversion matrix from XYZ colour space to a linear sRGB space with
 * coordinates normalised to [0, 1] range.
 */
const rgbFromXyzMatrix: Matrix = [
	[ 3.240812398895283,    -1.5373084456298136,  -0.4985865229069666],
	[-0.9692430170086407,    1.8759663029085742,   0.04155503085668564],
	[ 0.055638398436112804, -0.20400746093241362,  1.0571295702861434]
];

/**
 * Performs an sRGB gamma expansion of an 8-bit value, i.e. an integer in [0,
 * 255] range, into a floating-point value in [0, 1] range.
 */
function gammaExpansion(value255: number): number {
	return value255 <= 10
		? value255 / 3294.6
		: Math.pow((value255 + 14.025) / 269.025, 2.4);
}

/**
 * Performs an sRGB gamma compression of a floating-point value in [0, 1] range
 * into an 8-bit value, i.e. an integer in [0, 255] range.
 */
function gammaCompression(linear: number): number {
	let nonLinear: number = linear <= 0.00313066844250060782371
		? 3294.6 * linear
		: (269.025 * Math.pow(linear, 5.0 / 12.0) - 14.025);
	return Math.round(nonLinear) | 0;
}

/**
 * Multiplies a 3✕3 matrix by a 3✕1 column matrix.  The result is another 3✕1
 * column matrix.  The column matrices are represented as single-dimensional
 * 3-element array.  The matrix is represented as a two-dimensional array of
 * rows.
 */
function matrixMultiplication3x3x1(matrix: Matrix, column: Tripple): Tripple {
	return matrix.map((row: Tripple) => (
		row[0] * column[0] + row[1] * column[1] + row[2] * column[2]
	)) as Tripple;
}

/**
 * Converts sRGB colour given as a triple of 8-bit integers into XYZ colour
 * space.
 */
function xyzFromRgb(rgb: Tripple): Tripple {
	return matrixMultiplication3x3x1(
		xyzFromRgbMatrix, rgb.map(gammaExpansion) as Tripple);
}

/**
 * Converts colour from XYZ space to sRGB colour represented as a triple of
 * 8-bit integers.
 */
function rgbFromXyz(xyz: Tripple): Tripple {
	return matrixMultiplication3x3x1(
		rgbFromXyzMatrix, xyz).map(gammaCompression) as Tripple;
}

Updated in March 2021 with more precise value for the D65 standard illuminant. This affected values in xyzFromRgbMatrix and rgbFromXyzMatrix matrices.