Calculating RGB↔XYZ matrix

Posted by Michał ‘mina86’ Nazarewicz on 3rd of February 2019

I’ve recently found myself in need of an RGB↔XYZ transformation matrix expressed to the maximum possible precision. Sources on the Internet typically limit the precision to a handful decimal places so I’ve performed do the calculations myself.

What we’re looking for is a 3-by-3 matrix $M$ which, when multiplied by red, green and blue coordinates of a colour, produces its XYZ coordinates. In other words, a change of basis matrix from a space whose basis vectors are RGB’s primary colours: $$ M = \begin{bmatrix} X_r & X_g & X_b \\ Y_r & Y_g & Y_b \\ Z_r & Z_g & Z_b \end{bmatrix} $$

Derivation

RGB primary colours are typically provided as chromaticity usually expressed as x and y coordinates. For sRGB they are defined in IEC 61966 standard (and also Rec. 709 document which is 170 francs cheaper). Converting them to XYZ space is simple: $X = x Y / y$ and $Z = (1 - x - y) Y / y,$ but leaves luminosity (the Y value) unknown.

	$\langle x, y\rangle$	$\langle X, Y, Z\rangle$
Red	$\langle 0.64, 0.33\rangle$	$\langle {64 \over 33} Y_r, \; Y_r, \; {1 \over 11} Y_r\rangle$
Green	$\langle 0.30, 0.60\rangle$	$\langle {1 \over 2} Y_g, \; Y_g, \; {1 \over 6} Y_g\rangle$
Blue	$\langle 0.15, 0.06\rangle$	$\langle {5 \over 2} Y_b, \; Y_b, \; {79 \over 6} Y_b\rangle$
White (D65)	$\langle 0.312713, 0.329016\rangle$	$\langle {312713 \over 329016} Y_w, \; Y_w, \; {358271 \over 329016} Y_w\rangle$

That’s where reference white point comes into play. Its coordinates in linear RGB space, $\langle 1, 1, 1 \rangle,$ can be plugged into the change of basis formula to yield the following equation: $$ \begin{bmatrix} X_w \\ Y_w \\ Z_w \end{bmatrix} = \begin{bmatrix} X_r & X_g & X_b \\ Y_r & Y_g & Y_b \\ Z_r & Z_g & Z_b \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} $$

For each colour $c$ (including white), $X_c$ and $Z_c$ can be expressed as a product of a known quantity and $Y_c$ (see table above). Furthermore, by definition of a white point, $Y_w = 1.$ At this point, luminosities of the primary colours are the only unknowns. To isolate them, let’s define $X'_c = X_c / Y_c$ and $Z'_c = Z_c / Y_c$ and see where that leads us: $$ \begin{align} \begin{bmatrix} X_w \\ Y_w \\ Z_w \end{bmatrix} &= \begin{bmatrix} X_r & X_g & X_b \\ Y_r & Y_g & Y_b \\ Z_r & Z_g & Z_b \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \\ &= \begin{bmatrix} X'_r Y_r & X'_g Y_g & X'_b Y_b \\ Y_r & Y_g & Y_b \\ Z'_r Y_r & Z'_g Y_g & Z'_b Y_b \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \\ &= \begin{bmatrix} X'_r & X'_g & X'_b \\ 1 & 1 & 1 \\ Z'_r & Z'_g & Z'_b \end{bmatrix} \begin{bmatrix} Y_r & 0 & 0 \\ 0 & Y_g & 0 \\ 0 & 0 & Y_b \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix} \\ &= \begin{bmatrix} X'_r & X'_g & X'_b \\ 1 & 1 & 1 \\ Z'_r & Z'_g & Z'_b \end{bmatrix} \begin{bmatrix} Y_r \\ Y_g \\ Y_b \end{bmatrix} \\ \\ \begin{bmatrix} Y_r \\ Y_g \\ Y_b \end{bmatrix} &= \begin{bmatrix} X'_r & X'_g & X'_b \\ 1 & 1 & 1 \\ Z'_r & Z'_g & Z'_b \end{bmatrix}^{-1} \begin{bmatrix} X_w \\ Y_w \\ Z_w \end{bmatrix} \end{align} $$

All quantities on the right-hand side are known, therefore $[Y_r \; Y_g \; Y_b]^T$ can be computed. Let’s tidy things up into a final formula.

Final formula

Given chromaticity of primary colours of an RGB space ($\langle x_r, y_r \rangle,$ $\langle x_g, y_g \rangle$ and $\langle x_b, y_b \rangle$) and its reference white point ($\langle x_w, y_w \rangle$), the matrix for converting linear RGB coordinates to XYZ is: $$ M = \begin{bmatrix} X'_r Y_r & X'_g Y_g & X'_b Y_b \\ Y_r & Y_g & Y_b \\ Z'_r Y_r & Z'_g Y_g & Z'_b Y_b \end{bmatrix} $$

which can also be written as $M = M' \times \mathrm{diag}(Y_r, Y_g, Y_b)$ where: $$\begin{align} & M' = \begin{bmatrix} X'_r & X'_g & X'_b \\ 1 & 1 & 1 \\ Z'_r & Z'_g & Z'_b \end{bmatrix}\!\!, \\ \\ & \left. \begin{array}{l} X'_c = x_c / y_c \\ Z'_c = (1 - x_c - y_c) / y_c \end{array} \right\} \textrm{ for each colour } c, \\ \\ & \begin{bmatrix} Y_r \\ Y_g \\ Y_b \end{bmatrix} = M'^{-1} \begin{bmatrix} X_w \\ Y_w \\ Z_w \end{bmatrix} \textrm{ and} \\ \\ & \begin{bmatrix} X_w \\ Y_w \\ Z_w \end{bmatrix} = \begin{bmatrix} x_w / y_w \\ 1 \\ (1 - x_w - y_w) / y_w \end{bmatrix}\!\!. \end{align} $$

Matrix converting coordinates in the opposite direction is the inverse of $M,$ i.e. $M^{-1}$.

Implementation

Having the theoretical part covered, it’s time to put the equations into practice, which brings up a question of language to use. Rust programmers may simply reach for the rgb_derivation crate which provides all the necessary algorithms and supports rational arithmetic using big integers. However, for the sake of brevity, rather than presenting the Rust implementation, we’re going to use Python instead. The code begins with an overview of the algorithm:

Chromaticity = collections.namedtuple('Chromaticity', 'x y')

def calculate_rgb_matrix(primaries, white):
    M_prime = (tuple(c.x / c.y             for c in primaries),
               tuple(1                     for _ in primaries),
               tuple((1 - c.x - c.y) / c.y for c in primaries))
    W = (white.x / white.y, 1, (1 - white.x - white.y) / white.y)
    Y = mul_matrix_by_column(inverse_3x3_matrix(M_prime), W)
    return mul_matrix_by_diag(M_prime, Y)

The function first constructs $M'$ matrix and $W = [X_w\; Y_w\; Z_w]^T$ column which are used to calculate $[Y_r\; Y_g\; Y_b]^T = M'^{-1} W$ formula. With that computed, the function returns $M' \times \mathrm{diag}(Y_r, Y_g, Y_b)$ which is the transform matrix.

All operations on matrices are delegated to separate functions. Since the matrices are small there is no need for any sophisticated algorithms. Instead, the most straightforward multiplication algorithms are used:

def mul_matrix_by_column(matrix, column):
    return tuple(
        sum(row[i] * column[i] for i in range(len(row)))
        for row in matrix)

def mul_matrix_by_diag(matrix, column):
    return tuple(
        tuple(row[c] * column[c] for c in range(len(column)))
        for row in matrix)

Only the function inverting a 3-by-3 matrix is somewhat more complex:

def inverse_3x3_matrix(matrix):
    def cofactor(row, col):
        minor = [matrix[r][c]
                 for r in (0, 1, 2) if r != row
                 for c in (0, 1, 2) if c != col]
        a, b, c, d = minor
        det_minor = a * d - b * c
        return det_minor if (row ^ col) & 1 == 0 else -det_minor

    comatrix = tuple(
        tuple(cofactor(row, col) for col in (0, 1, 2))
        for row in (0, 1, 2))
    det = sum(matrix[0][col] * comatrix[0][col] for col in (0, 1, 2))
    return tuple(
        tuple(comatrix[col][row] / det for col in (0, 1, 2))
        for row in (0, 1, 2))

It first constructs a matrix of cofactors of the input (i.e. comatrix). Because function’s argument is always a 3-by-3 matrix, each minor’s determinant can be trivially calculated using $\bigl|\begin{smallmatrix} a & b \\ c & d \end{smallmatrix}\bigr| = a d - b c$ formula. Once the comatrix is constructed, calculating the determinant of the input matrix and its inverse becomes a matter of executing a few loops.

The above code works for any RGB system. To get result for sRGB its sRGB’s primaries and white point chromaticities need to be passed:

def calculate_srgb_matrix():
    primaries = (Chromaticity(fractions.Fraction(64, 100),
                              fractions.Fraction(33, 100)),
                 Chromaticity(fractions.Fraction(30, 100),
                              fractions.Fraction(60, 100)),
                 Chromaticity(fractions.Fraction(15, 100),
                              fractions.Fraction( 6, 100)))
    white = Chromaticity(fractions.Fraction(312713, 1000000),
                         fractions.Fraction(329016, 1000000))
    return calculate_rgb_matrix(primaries, white)

Full implementation with other bells and whistles can be found inside of the ansi_colours repository.

The matrix

Finally, there’s the matrix itself. Rust programmers are in luck again. The srgb crate provides the matrix along with other values and functions needed to work with sRGB colour space. If the values are needed in other programming languages, they can be copied from the following listing:

    ⎡ 0.4124108464885388   0.3575845678529519  0.18045380393360833 ⎤
M ≈ ⎢ 0.21264934272065283  0.7151691357059038  0.07218152157344333 ⎥
    ⎣ 0.019331758429150258 0.11919485595098397 0.9503900340503373  ⎦
    ⎡ 33786752 / 81924984   29295110 / 81924984    14783675 / 81924984  ⎤
M = ⎢  8710647 / 40962492   29295110 / 40962492     2956735 / 40962492  ⎥
    ⎣  4751262 / 245774952  29295110 / 245774952  233582065 / 245774952 ⎦

      ⎡  3.240812398895283    -1.5373084456298136  -0.4985865229069666  ⎤
M⁻¹ ≈ ⎢ -0.9692430170086407    1.8759663029085742   0.04155503085668564 ⎥
      ⎣  0.055638398436112804 -0.20400746093241362  1.0571295702861434  ⎦
      ⎡   4277208 / 1319795    -2028932 / 1319795   -658032 / 1319795  ⎤
M⁻¹ = ⎢ -70985202 / 73237775  137391598 / 73237775  3043398 / 73237775 ⎥
      ⎣    164508 / 2956735     -603196 / 2956735   3125652 / 2956735  ⎦

It’s important to remember this is not all that is needed to perform conversion between sRGB and XYZ colour spaces. The matrix assumes a linear RGB coordinates normalised to [0, 1] range which requires dealing with gamma correction. I’ve written a follow-up article which describes how to deal with that.

Update (March 2021): The article has been updated with more precise value for the D65 standard illuminant. This affected the resulting numbers slightly.

	\(\langle x, y\rangle\)	\(\langle X, Y, Z\rangle\)
Red	\(\langle 0.64, 0.33\rangle\)	\(\langle {64 \over 33} Y_r, \; Y_r, \; {1 \over 11} Y_r\rangle\)
Green	\(\langle 0.30, 0.60\rangle\)	\(\langle {1 \over 2} Y_g, \; Y_g, \; {1 \over 6} Y_g\rangle\)
Blue	\(\langle 0.15, 0.06\rangle\)	\(\langle {5 \over 2} Y_b, \; Y_b, \; {79 \over 6} Y_b\rangle\)
White (D65)	\(\langle 0.312713, 0.329016\rangle\)	\(\langle {312713 \over 329016} Y_w, \; Y_w, \; {358271 \over 329016} Y_w\rangle\)