/** This file contains functions for converting between colorspaces.

    See:
    Converting light frequency to RGB:
    https://stackoverflow.com/questions/1472514/convert-light-frequency-to-rgb

    Simple Analytic Approximations to the CIE XYZ Color Matching Functions: 
    https://jcgt.org/published/0002/02/01/

    CIE 1931 color space:
    https://en.wikipedia.org/wiki/CIE_1931_color_space

    http://www.brucelindbloom.com/index.html?Eqn_RGB_to_XYZ.html
*/



/** A multi-lobe, piecewise Gaussian fit of CIE 1931 XYZ Color Matching
    Functions by Wyman el al. The code is adopted from the Listing 1 (eq. 4) of

    Chris Wyman, Peter-Pike Sloan, and Peter Shirley, Simple Analytic
    Approximations to the CIE XYZ Color Matching Functions,
    Journal of Computer Graphics Techniques (JCGT), vol. 2, no. 2, 1-11, 2013.

    https://jcgt.org/published/0002/02/01/

    returns:
      [X, Y, Z]
*/
wavelengthToXYZ[wavelength is length] :=
{
   wave = wavelength / nm

   t1x = (wave - 442.0) * ((wave < 442.0) ? 0.0624 : 0.0374)
   t2x = (wave - 599.8) * ((wave < 599.8) ? 0.0264 : 0.0323)
   t3x = (wave - 501.1) * ((wave < 501.1) ? 0.0490 : 0.0382)

   x =   0.362 * exp[-0.5 t1x^2] + 
         1.056 * exp[-0.5 t2x^2] -
         0.065 * exp[-0.5 t3x^2]

   t1y = (wave - 568.8) * ((wave < 568.8) ? 0.0213 : 0.0247)
   t2y = (wave - 530.9) * ((wave < 530.9) ? 0.0613 : 0.0322)

   y =   0.821 * exp[-0.5 t1y^2] + 0.286 * exp[-0.5 t2y^2]
   
   t1z = (wave - 437.0) * ((wave < 437.0) ? 0.0845 : 0.0278)
   t2z = (wave - 459.0) * ((wave < 459.0) ? 0.0385 : 0.0725)

   z =   1.217 * exp[-0.5 t1z^2] + 0.681 * exp[-0.5 t2z^2]

   return [x,y,z]
}


/** Convert a wavelength to an RGB color in the sRGB space.

    returns:
      [R, G, B]
*/
wavelengthToRGB[wavelength is length] :=
{
   XYZ = wavelengthToXYZ[wavelength]
   return XYZtoRGB[XYZ]
}


/** Convert a color in the XYZ space to the sRGB color space.  The conversion
    matrix and color component transfer function is taken from

    http://www.color.org/srgb.pdf
    and
    http://www.brucelindbloom.com/index.html?Eqn_RGB_to_XYZ.html

    which follows the International Electrotechnical Commission standard
    IEC 61966-2-1 "Multimedia systems and equipment - Colour measurement and
    management - Part 2-1: Colour management - Default RGB colour space - sRGB"

    arguments:
      XYZ: an array containing the XYZ components from 0 to 1 [X, Y, Z]

    returns an array:
      [R, G, B]
*/
XYZtoRGB[XYZ] :=
{
   [x, y, z] = XYZ
   rl =  3.2406255 x + -1.537208  y + -0.4986286 z
   gl = -0.9689307 x +  1.8757561 y +  0.0415175 z
   bl =  0.0557101 x + -0.2040211 y +  1.0569959 z

   // TODO:  Clamp output components from 0 to 1?
   return [postProcessXYZ[rl],
           postProcessXYZ[gl],
           postProcessXYZ[bl]]
}

/** Apply color component transfer function for the above function. */
postProcessXYZ[c] :=
{
   c = clamp[c, 0, 1]
   return (c <= 0.0031308 ? 12.92 c : 1.055 c^(1/2.4) - 0.055)
}

/** Convert a color in the sRGB space to the XYZ space.  The conversion
    matrix and color component transfer function is taken from

    http://www.color.org/srgb.pdf
    and
    http://www.brucelindbloom.com/index.html?Eqn_RGB_to_XYZ.html

    which follows the International Electrotechnical Commission standard
    IEC 61966-2-1 "Multimedia systems and equipment - Colour measurement and
    management - Part 2-1: Colour management - Default RGB colour space - sRGB"

    arguments:
       RGB: an array containing the RGB components from 0 to 1 [R, G, B]

    returns:
       [X, Y, Z]
*/
RGBtoXYZ[RGB] :=
{
   [R, G, B] = RGB
   rl = preprocessRGB[R]
   gl = preprocessRGB[G]
   bl = preprocessRGB[B]

   X = 0.4124564 rl + 0.3575761 gl + 0.1804375 bl
   Y = 0.2126729 rl + 0.7151522 gl + 0.0721750 bl
   Z = 0.0193339 rl + 0.1191920 gl + 0.9503041 bl

   // TODO:  Clamp each channel between 0 and 1?
   return [X, Y, Z]
}

/** Perform inverse sRGB companding.  See:
    http://www.brucelindbloom.com/index.html?Eqn_RGB_to_XYZ.html
*/
preprocessRGB[V] :=
{
   V = clamp[V, 0, 1]
   if V <= 0.04045
      return V / 12.92
   else
      return ((V + 0.055) / 1.055)^2.4
}


/** Test function to plot wavelengths */
plotWavelengths[] :=
{
   gr = new graphics
   gr.backgroundColor[0,0,0]
   rl = new polyline
   gl = new polyline
   bl = new polyline
   g2 = new graphics
   g2.backgroundColor[0,0,0]
   g2.antialiased[false]
   stepsize = 1 nm
   for lambda = 375 nm to 725 nm step stepsize
   {
      [r,g,b] = wavelengthToRGB[lambda]
      rl.addPoint[lambda/nm/100, -r]
      gl.addPoint[lambda/nm/100, -g]
      bl.addPoint[lambda/nm/100, -b]
      g2.color[clamp[r,0,1],clamp[g,0,1],clamp[b,0,1]]
      g2.fillRectCenter[lambda/nm/10, 0, stepsize/nm/10, 10]
      println[format[lambda, "nm", 0] + "\t" + format[wavelengthToRGB[lambda], 1, 5]]
   }

   gr.color[1,0,0]
   gr.add[rl]
   gr.color[0,1,0]
   gr.add[gl]
   gr.color[0,0,1]
   gr.add[bl]
   gr.show[]
   g2.show[]
}