On the conversion from OSA UCS to CIEXYZ

Nico Schlömer

On the conversion from OSA-UCS to CIEXYZ

Abstract. This article revisits Kobayasi's and Yosiki's algorithm for conversion of OSA-UCS into XYZ cooordinates. It corrects some mistakes on the involved functions and initial guesses and shows that that hundreds of thousands of coordinates can be converted in less than a second with full accuracy.

In 1974, MacAdam published the definition of the OSA-UCS color space that tries to adhere particularly well to experimentally measured color distances. It combines work that had been going on since the late 1940s. One aspect of OSA-UCS is that, while the conversion from CIEXYZ coordinates into OSA-UCS $Lgj$ coordinates is straightforward, the conversion the other way around is not. In fact, there is no conversion method that works solely in elementary functions. Apparently, this had not been a design goal of OSA-UCS although is severely limits the usability of OSA-UCS.

In 2002, Kobayasi and Yosiki presented an algorithm for conversion from $Lgj$ to $XYZ$ coordinates that leverages Newton's method for solving nonlinear equation systems. Unfortunately, the article remains vague at important points and also contains false assertions about the nature of the involved functions.

In 2013, Cao et al. compared Kobayasi's and Yosiki's approach with some other, more complex methods based on artificial neural networks and found the latter to be superior.

In the present note, the author aims to iron out the inaccuracies in Kobayasi's article and improves the efficiency of the algorithm.

The forward conversion

The conversion from CIEXYZ-100 coordinates to OSA-UCS $Lgj$ coordinates is defined as follows:

• Compute $x$, $y$ coordinates via

  x = \frac{X}{X + Y + Z},\quad y = \frac{Y}{X + Y + Z}.
  

• Compute $K$ and $Y_0$ as

  \label{0}\tag{0}
  \begin{split}
    K &= 4.4934 x^2 + 4.3034 y^2 - 4.276 x y - 1.3744 x - 2.5643 y + 1.8103,\\
    Y_0 &= Y K.
  \end{split}
  

• Compute $L'$ and $C$ as

  \begin{split}
  \label{1}\tag{1}
  L' &= 5.9 \left(\sqrt[3]{Y_0} - \frac{2}{3} + 0.042 \sqrt[3]{Y_0 - 30}\right),\\
  C &= \frac{L'}{5.9 \left(\sqrt[3]{Y_0} - \frac{2}{3}\right)}.
  \end{split}
  

  (Note that $L'$ is $L$ in the original article.)

• Compute RGB as

  \label{2}\tag{2}
  \begin{bmatrix}
    R\\
    G\\
    B
  \end{bmatrix}
  =
  M
  \begin{bmatrix}
    X\\
    Y\\
    Z
  \end{bmatrix}
  \quad\text{with}\quad
  M=\begin{bmatrix}
    +0.7990 & 0.4194 & -0.1648\\
    -0.4493 & 1.3265 & +0.0927\\
    -0.1149 & 0.3394 & +0.7170
  \end{bmatrix}.
  

• Compute $a$, $b$ as

  \label{3}\tag{3}
  \begin{bmatrix}
    a\\
    b
  \end{bmatrix}
  =
  A
  \begin{bmatrix}
    \sqrt[3]{R}\\
    \sqrt[3]{G}\\
    \sqrt[3]{B}
  \end{bmatrix}
  \quad\text{with}\quad
  A = \begin{bmatrix}
    -13.7 & +17.7 & -4\\
    1.7 & +8 & -9.7
  \end{bmatrix}.
  

• Compute $L$, $g$, $j$ as

  L = \frac{L' - 14.3993}{\sqrt{2}},\quad g = Ca,\quad j = Cb.
  

The backward conversion

This section describes the conversion from the $Lgj$ to the $XYZ$ coordinates.

Given $L$, we can first compute

$$L' = L \sqrt{2} + 14.3993.$$

Equation $\eqref{1}$ gives the nonlinear relationship between $L'$ and $Y_0$ from which we will retrieve $Y_0$. First set $t\coloneqq \sqrt[3]{Y_0}$ and consider

$$\label{4}\tag{4} 0 = f(t) \coloneqq {\left(\frac{L'}{5.9} + \frac{2}{3} - t\right)}^3 - 0.042^3 (t^3 - 30).$$

$f$ is a monotonically decreasing cubic polynomial (see figure \ref{fig:singularity}).

Hence, it has exactly one root that can be found using the classical Cardano formula:

• Expand $f(t) = at^3 + bt^2 + ct + d$ with

  \begin{split}
    &u = \frac{L'}{5.9} + \frac{2}{3},\quad v = 0.042^3,\\
    &a = -(v + 1),\quad  b = 3u,\quad  c = -3u^2, \quad d = u^3 + 30v.
  \end{split}
  

• Compute the depressed form $\tilde{f}(x)=a(x^3 + px + q)$:

  p = \frac{3ac - b^2}{3a^2},\quad q = \frac{2b^3 - 9abc + 27a^2d}{27a^3}.
  

• Compute the root as

  t = -\frac{b}{3a} + \sqrt[3]{
    -\frac{q}{2} + \sqrt{{\left(\frac{q}{2}\right)}^2 + {\left(\frac{p}{3}\right)}^3}
  }
  + \sqrt[3]{
    -\frac{q}{2} - \sqrt{{\left(\frac{q}{2}\right)}^2 + {\left(\frac{p}{3}\right)}^3}
  }.
  

  Note that the expression in the square root, ${\left(q/2\right)}^2 + {\left(p/3\right)}^3$ is always positive since, as argued above, $f$ has exactly one root.

Remark. Kobayasi and Yosiki find the root of $f$ using Newton's method. A good initial guess here is $t = \frac{L'}{5.9} + \frac{2}{3}$ since the second term in $f(t)$, containing $0.042^3$, is comparatively small. Indeed it typically only takes around 10 iterations to converge to machine precision. Cardano's method finds the root at once at the expense of computing one square root and two cube roots. This approach is found to be about 15 times faster.

{ "layout": {"width": 400, "height": 300, "margin": {"l": 30, "r": 30, "b": 30, "t": 30, "pad": 4}}, "data": [{ "x": [4.0, 4.02, 4.04, 4.06, 4.08, 4.1, 4.12, 4.14, 4.16, 4.18, 4.2, 4.22, 4.24, 4.26, 4.28, 4.3, 4.32, 4.34, 4.36, 4.38, 4.4, 4.42, 4.44, 4.46, 4.48, 4.5, 4.52, 4.54, 4.56, 4.58, 4.6, 4.62, 4.64, 4.66, 4.68, 4.7, 4.72, 4.74, 4.76, 4.78, 4.8, 4.82, 4.84, 4.86, 4.88, 4.9, 4.92, 4.94, 4.96, 4.98, 5.0, 5.02, 5.04, 5.06, 5.08, 5.1, 5.12, 5.14, 5.16, 5.18, 5.2, 5.22, 5.24, 5.26, 5.28, 5.3, 5.32, 5.34, 5.36, 5.38, 5.4, 5.42, 5.44, 5.46, 5.48, 5.5, 5.52, 5.54, 5.56, 5.58, 5.6, 5.62, 5.64, 5.66, 5.68, 5.7, 5.72, 5.74, 5.76, 5.78, 5.8, 5.82, 5.84, 5.86, 5.88, 5.9, 5.92, 5.94, 5.96, 5.98, 6.0], "y": [0.785235946109928, 0.735087280230532, 0.687106436386218, 0.641245411020767, 0.59745620057795, 0.555690801501548, 0.515901210235332, 0.478039423223082, 0.442057436908571, 0.407907247735577, 0.375540852147875, 0.344910246589242, 0.315967427503453, 0.288664391334285, 0.262953134525512, 0.238785653520913, 0.216113944764262, 0.194890004699335, 0.175065829769908, 0.156593416419758, 0.13942476109266, 0.123511860232391, 0.108806710282726, 0.0952613076874412, 0.0828276488903126, 0.0714577303351169, 0.0611035484656296, 0.0517170997256262, 0.0432503805588832, 0.0356553874091769, 0.028884116720283, 0.0228885649359772, 0.017620728500036, 0.0130326038562349, 0.00907618744835042, 0.00570347572015815, 0.00286646511543439, 0.000517152077954889, -0.00139246694850414, -0.00291039552016687, -0.00408463719325717, -0.00496319552399913, -0.00559407406861669, -0.0060252763833339, -0.00630480602437471, -0.00648066654796317, -0.00660086151032324, -0.00671339446767894, -0.00686626897625427, -0.00710748859227323, -0.0074850568719598, -0.00804697737153799, -0.00884125364723183, -0.00991588925526531, -0.0113188877518624, -0.013098252693247, -0.0153019876356434, -0.0179780961352754, -0.021174581748367, -0.0249394480311421, -0.029320698539825, -0.0343663368306394, -0.0401243664598096, -0.0466427909835592, -0.0539696139581127, -0.0621528389396935, -0.0712404694845263, -0.0812805091488343, -0.0923209614888424, -0.104409830060774, -0.117595118420853, -0.131924830125303, -0.147446968730349, -0.164209537792216, -0.182260540867126, -0.201647981511303, -0.222419863280971, -0.244624189732356, -0.268308964421681, -0.293522190905168, -0.320311872739043, -0.34872601347953, -0.378812616682854, -0.410619685905235, -0.4441952247029, -0.479587236632074, -0.516843725248979, -0.55601269410984, -0.59714214677088, -0.640280086788325, -0.685474517718396, -0.73277344311732, -0.782224866541318, -0.833876791546617, -0.887777221689438, -0.943974160526009, -1.00251561161255, -1.06344957850529, -1.12682406476044, -1.19268707393425, -1.26108660958291], "type": "scatter" }] }

{ "layout": {"width": 400, "height": 300, "margin": {"l": 30, "r": 30, "b": 30, "t": 30, "pad": 4}, "yaxis": { "range": [-100, 100] }}, "data": [{ "x": [-1.0, -0.939393939393939, -0.878787878787879, -0.818181818181818, -0.757575757575758, -0.696969696969697, -0.636363636363636, -0.575757575757576, -0.515151515151515, -0.454545454545455, -0.393939393939394, -0.333333333333333, -0.272727272727273, -0.212121212121212, -0.151515151515151, -0.0909090909090908, -0.0303030303030303, 0.0303030303030303, 0.0909090909090908, 0.151515151515152, 1.18181818181818, 1.24242424242424, 1.3030303030303, 1.36363636363636, 1.42424242424242, 1.48484848484848, 1.54545454545455, 1.60606060606061, 1.66666666666667, 1.72727272727273, 1.78787878787879, 1.84848484848485, 1.90909090909091, 1.96969696969697, 2.03030303030303, 2.09090909090909, 2.15151515151515, 2.21212121212121, 2.27272727272727, 2.33333333333333, 2.39393939393939, 2.45454545454545, 2.51515151515152, 2.57575757575758, 2.63636363636364, 2.6969696969697, 2.75757575757576, 2.81818181818182, 2.87878787878788, 2.93939393939394, 3.0, 3.06060606060606, 3.12121212121212, 3.18181818181818, 3.24242424242424, 3.3030303030303, 3.36363636363636, 3.42424242424242, 3.48484848484849, 3.54545454545455, 3.60606060606061, 3.66666666666667, 3.72727272727273, 3.78787878787879, 3.84848484848485, 3.90909090909091, 3.96969696969697, 4.03030303030303, 4.09090909090909, 4.15151515151515, 4.21212121212121, 4.27272727272727, 4.33333333333333, 4.39393939393939, 4.45454545454546, 4.51515151515152, 4.57575757575758, 4.63636363636364, 4.6969696969697, 4.75757575757576, 4.81818181818182, 4.87878787878788, 4.93939393939394, 5.0], "y": [-107.413027646396, -107.309786831062, -107.215237122039, -107.127190104587, -107.042359529441, -106.954298360616, -106.847025918721, -106.673796503393, -106.283050417142, -105.164709350817, -101.697626835742, -91.7446608040379, -69.1567235939574, -33.4644220071798, 5.1529978523228, 38.1524502469182, 66.9186705474708, 97.8541264114677, 138.357047423149, 197.483926762583, 162.475784128249, 125.967621812498, 98.2246601475964, 76.6635115820041, 59.5986674360327, 45.8913206375379, 34.748884678301, 25.6055084607501, 18.048190278627, 11.7696062648161, 6.53714074322924, 2.17204573916972, -1.46489446641638, -4.48341990771192, -6.97097581197619, -8.99792296863166, -10.6213738820869, -11.8880550669845, -12.836469002897, -13.4985456088261, -13.9009169725394, -14.0659108181188, -14.0123317432564, -13.7560807176327, -13.310650174996, -12.6875225826085, -11.8964935122321, -10.9459352060274, -9.84301290556648, -8.59386342878862, -7.20374338230765, -5.6771528032405, -4.01793880526886, -2.22938286355019, -0.314274643424397, 1.72502529216217, 3.88653200963904, 6.16859043355976, 8.56983602610926, 11.0891608914696, 13.7256844599801, 16.4787280530134, 19.3477927478875, 22.3325400585856, 25.4327750269647, 28.6484313839621, 31.9795584937718, 35.4263098382213, 38.9889328353565, 42.6677598169073, 46.4632000149567, 50.375732429673, 54.405899468091, 58.5543012592435, 62.8215905639128, 67.2084682082878, 71.7156789801974, 76.3440079346035, 81.0942770619015, 85.9673422784606, 90.9640907039071, 96.0854381940253, 101.332327101924, 106.705724243391], "type": "scatter" }] }

Graph of $f(t)$ $\eqref{4}$ for $L'=25$. Note that the root is not in the turning point, but close to it. This is because of the small second term in $f$.

Graph of the function $\phi$ for $L$, $g$, $j$ computed from $X=12$, $Y=67$, $Z=20$. The singularity is at $w\approx 0.59652046418$. Note that the function has three roots only the largest of which is of interest.

From here, one can compute

$$\label{6}\tag{6} Y_0 = t^3,\quad C = \frac{L'}{5.9 \left(t - \frac{2}{3}\right)},\quad a = \frac{g}{C},\quad b = \frac{j}{C}.$$

With $a$ and $b$ at hand, it is now possible via equation $\eqref{3}$ to pin down $(\sqrt[3]{R}, \sqrt[3]{G}, \sqrt[3]{B})$ to only one degree of freedom, $w$. The exact value of $w$ will be found by Newton iteration. The function $\phi(w)$ of which a root needs to be found is defined as follows.

Append the matrix $A$ $\eqref{3}$ with a row such that the new $3\times3$-matrix $\tilde{A}$ is nonsingular and solve

$$\begin{bmatrix} a\\\ b\\\ w \end{bmatrix} = \tilde{A} \begin{bmatrix} \sqrt[3]{R}\\\ \sqrt[3]{G}\\\ \sqrt[3]{B} \end{bmatrix}$$

(Kobayasi, for instance, appends $[1, 0, 0]$ which corresponds to setting $w=\sqrt[3]{R}$.) Then compute the tentative $\tilde{X}$, $\tilde{Y}$, $\tilde{Z}$ via $\eqref{2}$ and further get the corresponding tentative $\tilde{Y}_0$ from $\eqref{0}$. Then $\phi(w) = \tilde{Y}_0(w) - Y_0$.

If the difference between $\tilde{Y}_0(w)$ and $Y_0$ from $\eqref{6}$ is 0, the correct $w$ has been found. Kobayasi states the function $\phi$ is "monotone increasing, convex downward, and smooth". Unfortunately, none of this is true (see figure \ref{fig:singularity}). In fact, the function has a singularity at $w$ chosen such that the computed tentative $\tilde{X}$, $\tilde{Y}$, $\tilde{Z}$ sum up to 0 while the individual values of $|\tilde{X}|, |\tilde{Y}|, |\tilde{Z}| > 0$. This happens if the tentative $[R, G, B]$ is orthogonal on $[1,1,1] M^{-1}$.

Fortunately, it seems that the function is indeed convex to the right of the singularity. Newton's method will hence find the correct (largest) root if the initial guess $w_0$ is chosen larger than the root. Since $w$ corresponds to $\sqrt[3]{R}$, it is reasonable to chose $w_0$ to be the maximum possible value that $\sqrt[3]{R}$ can take, namely that corresponding to $X=Y=100$, $Z=0$ (see $\eqref{2}$), $w_0=\sqrt[3]{79.9 + 41.94}\approx 4.9575$.

Cao et al. found that the conversion to from $Lgj$ to $XYZ$ takes so long that alternative methods need to be researched. They even find that the Newton iterations sometimes do not converge, or find the correct result only to a few digits of accuracy. The author cannot confirm these observations. The computation of hundreds of thousands of coordinates at once merely takes a second of computation time on a recent computer (figure \ref{fig:speed}).

To achieve this speed, it is important to vectorize all computation, i.e., not to perform the conversion for each $Lgj$-tuple individually one after another, but to perform all steps on the array. This also means to perform the Newton iteration on all tuples until the last one of them has converged successfully, even if some already converge in the first step. The redundant work inflicted by this approach is far outweighed by the advantages of vectorization.

All code is published as open-source in colorio.

{ "layout": { "width": 400, "height": 300, "margin": {"l": 30, "r": 30, "b": 30, "t": 30, "pad": 4}, "xaxis": {"type": "log", "title": {"text": "Size of the array"}}, "yaxis": {"type": "log", "title": {"text": "Runtime [s]"}} }, "data": [ { "name": "OSA-UCS", "x": [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, 512.0, 1024.0, 2048.0, 4096.0, 8192.0, 16384.0, 32768.0, 65536.0, 131072.0, 262144.0, 524288.0, 1048576.0, 2097152.0, 4194304.0], "y": [0.000471044, 0.000493364, 0.000493436, 0.000496653, 0.000569121, 0.000455073, 0.000551353, 0.000530194, 0.000632414, 0.001044303, 0.001434569, 0.001990027, 0.003218893, 0.005218656, 0.013776633, 0.044884756, 0.089765479, 0.162647182, 0.408672322, 1.01199472, 1.894782626, 3.468689603, 7.81700846], "type": "scatter" }, { "name": "CIECAM02", "x": [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, 512.0, 1024.0, 2048.0, 4096.0, 8192.0, 16384.0, 32768.0, 65536.0, 131072.0, 262144.0, 524288.0, 1048576.0, 2097152.0, 4194304.0], "y": [8.9266e-05, 9.2652e-05, 9.3239e-05, 9.3739e-05, 9.9035e-05, 9.8968e-05, 0.000105689, 0.000118602, 0.00014335, 0.000189389, 0.000271935, 0.00042539, 0.000735595, 0.001374689, 0.002681221, 0.008891838, 0.02621592, 0.049708766, 0.095939109, 0.194949605, 0.396554787, 0.851998144, 1.165202175], "type": "scatter" }, { "name": "CIELAB", "x": [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, 512.0, 1024.0, 2048.0, 4096.0, 8192.0, 16384.0, 32768.0, 65536.0, 131072.0, 262144.0, 524288.0, 1048576.0, 2097152.0, 4194304.0], "y": [1.9605e-05, 2.0169e-05, 2.0567e-05, 2.0203e-05, 2.262e-05, 2.065e-05, 2.195e-05, 2.8066e-05, 3.9437e-05, 3.0937e-05, 3.8563e-05, 5.544e-05, 0.000155565, 0.00028523, 0.000265324, 0.000501102, 0.002711648, 0.007986971, 0.013183137, 0.025345891, 0.067077266, 0.08853523, 0.170746303], "type": "scatter" } ] }

{ "layout": { "width": 400, "height": 300, "margin": {"l": 30, "r": 30, "b": 30, "t": 30, "pad": 4}, "xaxis": {"type": "log", "title": {"text": "Size of the OSA-UCS array"}} }, "data": [ { "x": [1.0, 2.0, 4.0, 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, 512.0, 1024.0, 2048.0, 4096.0, 8192.0, 16384.0, 32768.0, 65536.0, 131072.0, 262144.0, 524288.0, 1048576.0, 2097152.0, 4194304.0], "y": [787.066433566434, 748.043062200957, 631.69647696477, 540.053088803089, 430.216724738676, 197.28831378892, 135.024555079973, 69.6830099941211, 42.7417648170309, 37.6977113374785, 25.8176709590747, 20.1362440723177, 14.446842284113, 11.6819040387925, 16.5864959757117, 29.5641483554035, 31.3942138250337, 24.4242208227254, 30.4130064747132, 37.476925997706, 34.5039382592277, 31.8295473370856, 39.7625721082598], "type": "scatter" } ] }

Comparison with CIELAB and CIECAM02. The conversion of several hundred thousand $Lgj$ values takes about 1 second.

Computation speed relative to the evaluation of the cubic root. For large arrays, the conversion to $XYZ$ is about as costly as the evaluation of 35 cubic roots.