PGFPlots: Mapping CIELab Color Space to a sphere by conversion to RGB color values

by Marius   Last Updated October 09, 2019 22:23 PM - source

After already posting a question about Drawing a sphere and mapping CIELab color space to it (i.e. generally asking if pgfplots allows an axis-dependent color mapping beforehand), I searched all the calculus and theory stuff about CIELab color space I needed, in order to get the equations I need for a (possible) color mapping.

My idea now is to take the x,y,z values of my sphere (see MWE below), where, for the CIELab Color Space, I consider my sphere coordinates (all calculated in spherical coordinates incl. sin and cos) to be x=a, y=b and z=L.

Starting with CIELab Color Space, the calculation should go as follows: L,a,b --> X,Y,Z --> R,G,B.

Using this, I first defined necessary functions with pgfmanual, p. 1032 / 1033 and

declare function = {<function definitions>}

plus ifthenelse (same package)

<statement> ? <yes> : <no> 

to calculate some coefficients xr,yr,zr which are required to perform L,a,b --> X,Y,Z as follows:

 declare function={% Coefficients xr,yr,zr for X,Y,Z calculation

    xr(\L,\a) = ( (\a/500) + ((\L+16)/116) )^3 > 0.008856 ?%
                ( ((\a/500) + ((\L + 16) / 116))^3 ) :%
                ( 116 * ((\a/500) + ((\L + 16) / 116)) - 16 ) / 903.3;

    yr(\L) = \L > ( 0.008856 * 903.3 ) ?%
             ((\L + 16)/116)^3 :%
             \L / 903.3;

    zr(\L,\b) = ( ((\L + 16) / 116) - (\b/200) )^3 > 0.008856 ?%
                ((((\L + 16) / 116) - (\b/200))^(3)) :%
                ( 116 *  (((\L + 16) / 116) - (\b/200)) - 16 ) / 903.3;
  }

Doing this, I acquire the X,Y,Z (tristimulus) values by multiplication of those coefficients xr,yr,zr with some illuminant constants Xn,Yn,Zn to achieve

  X = xr(\L,\a) * Xn
  Y = yr(\L)    * Yn
  Z = zr(\L,\b) * Zn

Following this, a matrix (matrix elements Mij)

       [M11 M12 M13]  -->  gives R value
  M =  [M21 M22 M23]  -->  gives G value
       [M31 M32 M33]  -->  gives B value

for XYZ --> RGB (both 3x1 vectors) conversion can be used to get the final results (functions R,G,B, depending on \L,\a,\b):

   R(\L,\a,\b) =   xr(\L,\a) * Xn * M11
                 + yr(\L)    * Yn * M12
                 + zr(\L,\b) * Zn * M13;

   G(\L,\a,\b) =   xr(\L,\a) * Xn * M21
                 + yr(\L)    * Yn * M22
                 + zr(\L,\b) * Zn * M23;

   B(\L,\a,\b) =   xr(\L,\a) * Xn * M31
                 + yr(\L)    * Yn * M32
                 + zr(\L,\b) * Zn * M33;

Having this, it should (theoretically / as I would think) be possible to simply plug in the coordinates x=\a [-100;100], y=\b [-100;100] and z=\L [0;100] of my sphere to get the results for each point meta in \addplot3:

 \addplot3[point meta = {symbolic = {<R>,<G>,<B>}}] (% Define sphere to be mapped on, incl. limited domains
            {100*cos(azimuth)*sin(polar)},% x
            {100*sin(azimuth)*sin(polar)},% y
            {50*cos(polar)+50}% z
 );

Unfortunately, the CIELab model does not exactly represent a sphere (rather an ellipsoid, since the three values of RGB (color gamut) can not fill the whole color space of XYZ or Lab completely, which is also a pretty nasty companion for every company producing colors in TVs, Smartphones etc.), which is why one also has to place exceptions when mapping RGB to a sphere.

Here, I simply defined the point meta = {symbolic = {<R>,<G>,<B>}} to replace values < 0 --> 0 and > 1 --> 1, in order to prevent the compiling process from crashing (using nested ifthenelse):

 point meta={%
   symbolic={%
      % R Values in [0;1]
      ifthenelse( R(z,x,y) < 0 , 0.0 , ifthenelse( R(z,x,y) > 1 , 1.0 , R(z,x,y) ) ),%
      % G Values [0;1]
      ifthenelse( G(z,x,y) < 0 , 0.0 , ifthenelse( G(z,x,y) > 1 , 1.0 , G(z,x,y) ) ),%
      % B Values [0;1]
      ifthenelse( B(z,x,y) < 0 , 0.0 , ifthenelse( B(z,x,y) > 1 , 1.0 , B(z,x,y) ) )%
   }%
 }

And here comes the thing: It's not compiling successfully, popping errors like Sorry, an internal routine of the floating point got an ill-formatted floating point number and even (for some reason) Unknown function 'Y' in 'xr(1,Y)' (might be the origin?) which leaves me helpless, since I am lacking the experience to fix it and to exactly know how pgfplots processes the data in mathparse.

Any ideas how to fix this / make it running successfully?

Here's my Code* so far (incl. concrete values plugged in and the commented mesh/colorspace explicit color input=rgb255-option, since I, though having the literature, am not sure (yet) whether R(\L,\a,\b), G(\L,\a,\b) and B(\L,\a,\b) will pop values in [0;1] or [0;255]; I would expect [0;1], though):

 \documentclass[tikz,border=3mm]{standalone}
   \usepackage{pgfplots}
     \pgfplotsset{compat=1.16}
       \usepgfplotslibrary{patchplots}

 \begin{document}
   \begin{tikzpicture}[%
     declare function={%
       % Calculation scheme:   Lab --> XYZ --> RGB [0;1]
       % Math, formulas and values based on
       %  - https://en.wikipedia.org/wiki/Illuminant_D65
       %  - https://en.wikipedia.org/wiki/CIELAB_color_space
       %  - https://en.wikipedia.org/wiki/Adobe_RGB_color_space
       %  - http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
       %  - http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
       %  - http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html (with XYZ -> sRGB)
       %
       % Declaring Functions to calculate CIE XYZ values --> RGB from the following variables:
       % \L = L-value of CIE Lab space
       % \a = a-value of CIE Lab space
       % \b = b-value of CIE Lab space
       % xr, yr, zr = Coefficients required to get XYZ (from L,a,b)
       %
       % D65 illuminant tristimulus values (T = 6504 K;)
       % Xn = 0.968774;
       % Yn = 1.0;
       % Zn = 1.121774;
       %
       % CIE constants (E = epsilon, K = kappa)
       % E = 0.008856;
       % K = 903.3;
       %
       % XYZ -> sRGB matrix:
       %       [3.2404542 -1.5371385 -0.4985314]   --> R Value
       %   M = [-0.9692660  1.8760108  0.0415560]  --> G Value
       %       [0.0556434 -0.2040259  1.0572252]   --> B Value
       %
       % Coefficients xr,yr,zr for X,Y,Z calculation
       xr(\L,\a) = ( (\a/500) + ((\L+16)/116) )^3 > 0.008856 ?% > E?
                   ( ( (\a/500) + ((\L+16) / 116) )^3 ) :%
                   ( 116 * ((\a/500) + ((\L+16)/116)) - 16 ) / 903.3;
       %
       yr(\L)    = \L > ( 0.008856 * 903.3 ) ?% < (E * K)?
                   ( (\L+16)/116 )^3 :%
                   \L / 903.3;
       %
       zr(\L,\b) = ( ((\L+16)/116) - (\b/200) )^3 > 0.008856 ?% > E?
                   ( ( ((\L+16)/116) - (\b/200) )^3 ) :%
                   ( 116 * (((\L+16)/116) - (\b/200)) - 16 ) / 903.3;%
       % Calculation of R,G,B via illuminant properties and matrix values (XYZ -> sRGB)
       R(\L,\a,\b) = xr(\L,\a) * 0.968774 * 3.2404542 + yr(\L) * 1.0 * (-1.5371385) + zr(\L,\b) * 1.121774 * (-0.4985314); % Including Xn, Yn, Zn, first matrix row
       G(\L,\a,\b) = xr(\L,\a) * 0.968774 * (-0.9692660) + yr(\L) * 1.0 * 1.8760108 + zr(\L,\b) * 1.121774 * 0.0415560;     % Including Xn, Yn, Zn, second matrix row
       B(\L,\a,\b) = xr(\L,\a) * 0.968774 * 0.0556434 + yr(\L) * 1.0 * (-0.2040259) + zr(\L,\b) * 1.121774 * 1.0572252;    % Including Xn, Yn, Zn, third matrix row
     }
   ]
   \begin{axis}[axis equal,
     width = 10cm,
     height = 10cm,
     axis lines = center,
     xmin = -120,
     xmax = 120,
     ymin = -120,
     ymax = 120,
     zmin = 0,
     zmax = 100,
     ticks = none,
     enlargelimits = 0.3,
     z buffer = sort,
     view/h = 45,
     scale uniformly strategy = units only]

     \addplot3 [%
       patch,
       patch type=bilinear,
       variable = \azimuth,
       variable y = \polar,
       domain = 0:360,
       y domain = 0:180,
       fill opacity = 0.5,
       draw opacity = 1,
       line width = 0.001 pt,
       samples = 10, % only for faster compilation
       mesh/color input=explicit mathparse,
       % mesh/colorspace explicit color input=rgb255, % if RGB values are calculated in [0 ; 255]
       point meta={%
         symbolic={%
           % R Values [0;1]
           ifthenelse( R(z,x,y) < 0 , 0.0 , ifthenelse( R(z,x,y) > 1 , 1.0 , R(z,x,y) ) ),% check r < 0 and r > 1
           % G Values [0;1]
           ifthenelse( G(z,x,y) < 0 , 0.0 , ifthenelse( G(z,x,y) > 1 , 1.0 , G(z,x,y) ) ),% check g < 0 and g > 1
           % B Values [0;1]
           ifthenelse( B(z,x,y) < 0 , 0.0 , ifthenelse( B(z,x,y) > 1 , 1.0 , B(z,x,y) ) )% check b < 0 and b > 1
         }%
       },
     ] (% Define sphere to be mapped on
          {100*cos(azimuth)*sin(polar)},%  x
          {100*sin(azimuth)*sin(polar)},%  y
          {50*cos(polar)+50}%              z
       );
    \end{axis}
   \end{tikzpicture}
 \end{document}

Thanks a lot for your ideas! :) - Marius.


*adapted from Schrödinger's cat's proposal in the previous question.



Related Questions


Pgfplots partial sphere

Updated March 21, 2017 17:23 PM

How to convert pgfplots to high quality PNG file?

Updated September 03, 2018 13:23 PM

Remapping German sharp s (Eszett)

Updated March 28, 2015 09:07 AM

Static map marking

Updated March 02, 2018 20:23 PM