posts/srgb-xyz-matrix.en.html

   1 <!-- subject: Calculating {RGB↔XYZ} matrix -->
   2 <!-- date: 2019-02-03 07:42:06 -->
   3 <!-- update: 2021-03-21 01:21:33 -->
   4 <!-- tags: rgb, srgb, colorimetry -->
   5 <!-- categories: Articles, Techblog -->
   6 <!-- math: true -->
   7
   8 <p>I’ve recently found myself in need of an RGB↔XYZ transformation matrix
   9   expressed to the maximum possible precision.  Sources on the Internet
  10   typically limit the precision to a handful decimal places so I’ve performed
  11   do the calculations myself.
  12
  13 <p>What we’re looking for is a 3-by-3 matrix \(M\) which, when multiplied by
  14   red, green and blue coordinates of a colour, produces its XYZ coordinates.  In
  15   other words,
  16   a <a href="https://www.khanacademy.org/math/linear-algebra/alternate-bases/change-of-basis/v/linear-algebra-change-of-basis-matrix">change
  17   of basis matrix</a> from a space whose basis vectors are RGB’s primary
  18   colours:
  19
  20   $$ M = \begin{bmatrix}
  21     X_r & X_g & X_b \\
  22     Y_r & Y_g & Y_b \\
  23     Z_r & Z_g & Z_b
  24   \end{bmatrix} $$
  25
  26 <!-- FULL -->
  27
  28
  29 <h2>Derivation</h2>
  30
  31 <p>RGB primary colours are typically provided as chromaticity usually expressed
  32   as <a href="https://en.wikipedia.org/wiki/CIE_1931_color_space#CIE_xy_chromaticity_diagram_and_the_CIE_xyY_color_space">x
  33   and y</a> coordinates.  For sRGB they are defined in
  34   <a href="https://webstore.iec.ch/publication/6169">IEC 61966 standard</a> (and
  35   also <a href="https://www.itu.int/rec/R-REC-BT.709-6-201506-I/en">Rec. 709
  36   document</a> which is 170 francs cheaper).  Converting them to XYZ space is
  37   simple: \(X = x Y / y\) and \(Z = (1 - x - y) Y / y,\) but leaves luminosity
  38   (the Y value) unknown.
  39
  40 <table>
  41   <thead>
  42     <tr>
  43       <th>
  44       <th scope=col>\(\langle x, y\rangle\)
  45       <th scope=col>\(\langle X, Y, Z\rangle\)
  46   <tbody>
  47     <tr>
  48       <th scope=row>Red
  49       <td>\(\langle 0.64, 0.33\rangle\)
  50       <td>\(\langle {64 \over 33} Y_r, \; Y_r, \; {1 \over 11} Y_r\rangle\)
  51     <tr>
  52       <th scope=row>Green
  53       <td>\(\langle 0.30, 0.60\rangle\)
  54       <td>\(\langle {1 \over 2} Y_g, \; Y_g, \; {1 \over 6} Y_g\rangle\)
  55       <tr>
  56         <th scope=row>Blue
  57         <td>\(\langle 0.15, 0.06\rangle\)
  58         <td>\(\langle {5 \over 2} Y_b, \; Y_b, \; {79 \over 6} Y_b\rangle\)
  59       <tr>
  60         <th scope=row>White (D65)
  61         <td>\(\langle 0.312713, 0.329016\rangle\)
  62         <td>\(\langle {312713 \over 329016} Y_w, \; Y_w, \; {358271 \over 329016} Y_w\rangle\)
  63 </table>
  64
  65 <p>That’s where reference white point comes into play.  Its coordinates in
  66   linear RGB space, \(\langle 1, 1, 1 \rangle,\) can be plugged into the change
  67   of basis formula to yield the following equation:
  68
  69   $$ \begin{bmatrix}
  70     X_w \\ Y_w \\ Z_w
  71   \end{bmatrix} = \begin{bmatrix}
  72     X_r & X_g & X_b \\
  73     Y_r & Y_g & Y_b \\
  74     Z_r & Z_g & Z_b
  75   \end{bmatrix} \begin{bmatrix}
  76     1 \\ 1 \\ 1
  77   \end{bmatrix} $$
  78
  79 <p>For each colour \(c\) (including white), \(X_c\) and \(Z_c\) can be expressed
  80   as a product of a known quantity and \(Y_c\) (see table above).  Furthermore,
  81   by definition of a white point, \(Y_w = 1.\) At this point, luminosities of
  82   the primary colours are the only unknowns.  To isolate them, let’s define
  83   \(X'_c = X_c / Y_c\) and \(Z'_c = Z_c / Y_c\) and see where that leads us:
  84
  85   $$ \begin{align}
  86     \begin{bmatrix}
  87       X_w \\ Y_w \\ Z_w
  88     \end{bmatrix} &=
  89     \begin{bmatrix}
  90       X_r & X_g & X_b \\
  91       Y_r & Y_g & Y_b \\
  92       Z_r & Z_g & Z_b
  93     \end{bmatrix} \begin{bmatrix}
  94       1 \\ 1 \\ 1
  95     \end{bmatrix} \\ &= \begin{bmatrix}
  96       X'_r Y_r & X'_g Y_g & X'_b Y_b \\
  97       Y_r      & Y_g      & Y_b \\
  98       Z'_r Y_r & Z'_g Y_g & Z'_b Y_b
  99     \end{bmatrix} \begin{bmatrix}
 100       1 \\ 1 \\ 1
 101     \end{bmatrix} \\ &= \begin{bmatrix}
 102       X'_r & X'_g & X'_b \\
 103       1    & 1    & 1 \\
 104       Z'_r & Z'_g & Z'_b
 105     \end{bmatrix} \begin{bmatrix}
 106       Y_r & 0 & 0 \\
 107       0 & Y_g & 0 \\
 108       0 & 0 & Y_b
 109     \end{bmatrix} \begin{bmatrix}
 110       1 \\ 1 \\ 1
 111     \end{bmatrix} \\ &= \begin{bmatrix}
 112       X'_r & X'_g & X'_b \\
 113       1    & 1    & 1 \\
 114       Z'_r & Z'_g & Z'_b
 115     \end{bmatrix} \begin{bmatrix}
 116       Y_r \\ Y_g \\ Y_b
 117     \end{bmatrix} \\ \\
 118
 119     \begin{bmatrix}
 120       Y_r \\ Y_g \\ Y_b
 121     \end{bmatrix} &= \begin{bmatrix}
 122       X'_r & X'_g & X'_b \\
 123       1    & 1    & 1 \\
 124       Z'_r & Z'_g & Z'_b
 125     \end{bmatrix}^{-1} \begin{bmatrix}
 126       X_w \\ Y_w \\ Z_w
 127     \end{bmatrix}
 128   \end{align} $$
 129
 130 <p>All quantities on the right-hand side are known, therefore \([Y_r \; Y_g \;
 131   Y_b]^T\) can be computed.  Let’s tidy things up into a final formula.
 132
 133
 134 <h2>Final formula</h2>
 135
 136 <p>Given chromaticity of primary colours of an RGB space (\(\langle x_r, y_r
 137   \rangle,\) \(\langle x_g, y_g \rangle\) and \(\langle x_b, y_b \rangle\)) and
 138   its reference white point (\(\langle x_w, y_w \rangle\)), the matrix for
 139   converting linear RGB coordinates to XYZ is:
 140
 141   $$ M = \begin{bmatrix}
 142     X'_r Y_r & X'_g Y_g & X'_b Y_b \\
 143     Y_r      & Y_g      & Y_b \\
 144     Z'_r Y_r & Z'_g Y_g & Z'_b Y_b
 145   \end{bmatrix} $$
 146
 147 <p>which can also be written as \(M = M' \times \mathrm{diag}(Y_r, Y_g, Y_b)\)
 148 where:
 149
 150   $$\begin{align}
 151     & M' = \begin{bmatrix}
 152       X'_r & X'_g & X'_b \\
 153       1    & 1    & 1 \\
 154       Z'_r & Z'_g & Z'_b
 155     \end{bmatrix}\!\!, \\ \\
 156
 157     & \left.
 158       \begin{array}{l}
 159          X'_c = x_c / y_c \\
 160          Z'_c = (1 - x_c - y_c) / y_c
 161       \end{array}
 162     \right\} \textrm{ for each colour } c, \\ \\
 163
 164     & \begin{bmatrix}
 165       Y_r \\ Y_g \\ Y_b
 166     \end{bmatrix} = M'^{-1} \begin{bmatrix}
 167       X_w \\ Y_w \\ Z_w
 168     \end{bmatrix} \textrm{ and} \\ \\
 169
 170     & \begin{bmatrix}
 171       X_w \\ Y_w \\ Z_w
 172     \end{bmatrix} = \begin{bmatrix}
 173       x_w / y_w \\ 1 \\ (1 - x_w - y_w) / y_w
 174     \end{bmatrix}\!\!.
 175   \end{align} $$
 176
 177 <p>Matrix converting coordinates in the opposite direction is the inverse of
 178   \(M,\) i.e. \(M^{-1}\).
 179
 180
 181 <h2>Implementation</h2>
 182
 183 <p>Having the theoretical part covered, it’s time to put the equations into
 184   practice, which brings up a question of language to use.  Rust programmers may
 185   simply reach for <a href="https://crates.io/crates/rgb_derivation">the
 186   rgb_derivation crate</a> which provides all the necessary algorithms and
 187   supports rational arithmetic using big integers.  However, for the sake of
 188   brevity, rather than presenting the Rust implementation, we’re going to use
 189   Python instead.  The code begins with an overview of the algorithm:
 190
 191 <pre>
 192 Chromaticity = collections.namedtuple('Chromaticity', 'x y')
 193
 194 def calculate_rgb_matrix(primaries, white):
 195     M_prime = (tuple(c.x / c.y             for c in primaries),
 196                tuple(1                     for _ in primaries),
 197                tuple((1 - c.x - c.y) / c.y for c in primaries))
 198     W = (white.x / white.y, 1, (1 - white.x - white.y) / white.y)
 199     Y = mul_matrix_by_column(inverse_3x3_matrix(M_prime), W)
 200     return mul_matrix_by_diag(M_prime, Y)
 201 </pre>
 202
 203 <p>The function first constructs \(M'\) matrix and \(W = [X_w\; Y_w\; Z_w]^T\)
 204   column which are used to calculate \([Y_r\; Y_g\; Y_b]^T = M'^{-1} W\)
 205   formula.  With that computed, the function returns \(M' \times
 206   \mathrm{diag}(Y_r, Y_g, Y_b)\) which is the transform matrix.
 207
 208 <p>All operations on matrices are delegated to separate functions.  Since the
 209   matrices are small there is no need for any sophisticated algorithms.
 210   Instead, the most straightforward multiplication algorithms are used:
 211
 212 <pre>
 213 def mul_matrix_by_column(matrix, column):
 214     return tuple(
 215         sum(row[i] * column[i] for i in range(len(row)))
 216         for row in matrix)
 217
 218 def mul_matrix_by_diag(matrix, column):
 219     return tuple(
 220         tuple(row[c] * column[c] for c in range(len(column)))
 221         for row in matrix)
 222 </pre>
 223
 224 <p>Only the function inverting a 3-by-3 matrix is somewhat more complex:
 225
 226 <pre>
 227 def inverse_3x3_matrix(matrix):
 228     def cofactor(row, col):
 229         minor = [matrix[r][c]
 230                  for r in (0, 1, 2) if r != row
 231                  for c in (0, 1, 2) if c != col]
 232         a, b, c, d = minor
 233         det_minor = a * d - b * c
 234         return det_minor if (row ^ col) & 1 == 0 else -det_minor
 235
 236     comatrix = tuple(
 237         tuple(cofactor(row, col) for col in (0, 1, 2))
 238         for row in (0, 1, 2))
 239     det = sum(matrix[0][col] * comatrix[0][col] for col in (0, 1, 2))
 240     return tuple(
 241         tuple(comatrix[col][row] / det for col in (0, 1, 2))
 242         for row in (0, 1, 2))
 243 </pre>
 244
 245 <p>It first constructs
 246   a <a href="https://en.wikipedia.org/wiki/Minor_(linear_algebra)#First_minors">matrix
 247   of cofactors</a> of the input (i.e. <code>comatrix</code>).  Because
 248   function’s argument is always a 3-by-3 matrix, each minor’s determinant can be
 249   trivially calculated using \(\bigl|\begin{smallmatrix} a & b \\ c &
 250   d \end{smallmatrix}\bigr| = a d - b c\) formula.  Once the comatrix is
 251   constructed, <a href="https://en.wikipedia.org/wiki/Minor_(linear_algebra)#Applications_of_minors_and_cofactors">calculating
 252   the determinant of the input matrix and its inverse</a> becomes a matter
 253   of executing a few loops.
 254
 255 <p>The above code works for any RGB system.  To get result for sRGB its sRGB’s
 256   primaries and white point chromaticities need to be passed:
 257
 258 <pre>
 259 def calculate_srgb_matrix():
 260     primaries = (Chromaticity(fractions.Fraction(64, 100),
 261                               fractions.Fraction(33, 100)),
 262                  Chromaticity(fractions.Fraction(30, 100),
 263                               fractions.Fraction(60, 100)),
 264                  Chromaticity(fractions.Fraction(15, 100),
 265                               fractions.Fraction( 6, 100)))
 266     white = Chromaticity(fractions.Fraction(312713, 1000000),
 267                          fractions.Fraction(329016, 1000000))
 268     return calculate_rgb_matrix(primaries, white)
 269 </pre>
 270
 271 <p>Full implementation with other bells and whistles can be
 272   found <a href="https://github.com/mina86/ansi_colours/blob/master/tools/luminance-coefficients.py">inside
 273   of the <code>ansi_colours</code> repository</a>.
 274
 275
 276 <h2>The matrix</h2>
 277
 278 <p>Finally, there’s the matrix itself.  Rust programmers are in luck
 279   again.  <a href="https://crates.io/crates/srgb">The srgb crate</a> provides
 280   <a href="https://docs.rs/srgb/0.1.4/srgb/xyz/constant.XYZ_FROM_SRGB_MATRIX.html">the
 281   matrix</a> along with other values and functions needed to work with sRGB
 282   colour space.  If the values are needed in other programming languages, they
 283   can be copied from the following listing:
 284
 285 <pre>
 286     ⎡ 0.4124108464885388   0.3575845678529519  0.18045380393360833 ⎤
 287 M ≈ ⎢ 0.21264934272065283  0.7151691357059038  0.07218152157344333 ⎥
 288     ⎣ 0.019331758429150258 0.11919485595098397 0.9503900340503373  ⎦
 289     ⎡ 33786752 / 81924984   29295110 / 81924984    14783675 / 81924984  ⎤
 290 M = ⎢  8710647 / 40962492   29295110 / 40962492     2956735 / 40962492  ⎥
 291     ⎣  4751262 / 245774952  29295110 / 245774952  233582065 / 245774952 ⎦
 292
 293       ⎡  3.240812398895283    -1.5373084456298136  -0.4985865229069666  ⎤
 294 M⁻¹ ≈ ⎢ -0.9692430170086407    1.8759663029085742   0.04155503085668564 ⎥
 295       ⎣  0.055638398436112804 -0.20400746093241362  1.0571295702861434  ⎦
 296       ⎡   4277208 / 1319795    -2028932 / 1319795   -658032 / 1319795  ⎤
 297 M⁻¹ = ⎢ -70985202 / 73237775  137391598 / 73237775  3043398 / 73237775 ⎥
 298       ⎣    164508 / 2956735     -603196 / 2956735   3125652 / 2956735  ⎦
 299 </pre>
 300
 301
 302 <p>It’s important to remember this is not all that is needed to perform
 303   conversion between sRGB and XYZ colour spaces.  The matrix assumes a linear
 304   RGB coordinates normalised to [0, 1] range which requires dealing with gamma
 305   correction.  I’ve written <a href="/2019/srgb-xyz-conversion/">a follow-up
 306   article</a> which describes how to deal with that.
 307
 308 <p>Updated in March 2021 with more precise value for the D65 standard
 309   illuminant.  This affected the resulting numbers slightly.