// Title: CIE Color Conversion Formulas
// Authors: Michael Horvath, Christoph Lipka
// Website: http://isometricland.net
// Created: 2016-11-20
// Updated: 2017-03-08
// This file is licensed under the terms of the CC-GNU LGPL.
// http://www.gnu.org/licenses/lgpl-2.1.html
// Illuminant =  D65
// Observer   =  2° (1931)
// +kfi0 +kff15 +kc


// constants
#declare cie_XYZWhiteReference1 =  95.047;
#declare cie_XYZWhiteReference2 = 100.000;
#declare cie_XYZWhiteReference3 = 108.883;
#declare cie_XYZEpsilon = 0.008856;
#declare cie_XYZKappa = 903.3;
#declare cie_uPrime = 4 * cie_XYZWhiteReference1/(cie_XYZWhiteReference1 + 15 * cie_XYZWhiteReference2 + 3 * cie_XYZWhiteReference3);
#declare cie_vPrime = 9 * cie_XYZWhiteReference2/(cie_XYZWhiteReference1 + 15 * cie_XYZWhiteReference2 + 3 * cie_XYZWhiteReference3);


// utility
#declare cie_selectRGBa = function(C,D) {select(C,-1,D)}
#declare cie_selectRGBb = function(C,D) {select(C-1,D,D,-1)}
#declare cie_linearRGBa = function(C) {C/12.92}
#declare cie_linearRGBb = function(C) {pow((C + 0.055)/(1 + 0.055), 2.4)}
#declare cie_linearRGBc = function(C) {select(C-0.04045,cie_linearRGBa(C),cie_linearRGBa(C),cie_linearRGBb(C))}
#declare cie_fD = function(C) {abs(C-0.5)-0.5}
#declare cie_fDist = function(Dr,Dg,Db) {max(Dr,Dg,Db)}
//#declare cie_fClip1 = function(X,A) {select(X-A,A,X)}
//#declare cie_fClip2 = function(X,A) {select(-(X-A),A,X)}
//#declare cie_fD = function(C) {abs(fClip2(fClip1(C,-0.1),1.1)-0.5)-0.5}
//#declare cie_fD = function(C) {abs(fClip1(C,-0.1)-0.5)-0.5}


// input L = between 0 and 100
// input C = between 0 and 128
// input H = between 0 and 360
// output L = between 0 and 100
// output A = between -128 and +128
// output B = between -128 and +128
#declare cie_funcLCH2LABa1 = function(L,C,H) {L}
#declare cie_funcLCH2LABa2 = function(L,C,H) {cosd(H) * C}
#declare cie_funcLCH2LABa3 = function(L,C,H) {sind(H) * C}


// input H = between 0 and 360
// input C = between 0 and 200
// input L = between 0 and 100
// output L = between 0 and 100
// output U = between -128 and +128
// output V = between -128 and +128
#declare cie_funcHCL2LUVa1 = function(H,C,L) {L}
#declare cie_funcHCL2LUVa2 = function(H,C,L) {cosd(H) * C}
#declare cie_funcHCL2LUVa3 = function(H,C,L) {sind(H) * C}


// input L = between 0 and 100
// input A = between -128 and +128
// input B = between -128 and +128
// output X: between 0 and 100
// output Y: between 0 and 100
// output Z: between 0 and 100
// Note that some of these cie_functions are out of order.
#declare cie_funcLAB2XYZa2 = function(L,A,B) {(L + 16) / 116}
#declare cie_funcLAB2XYZa1 = function(L,A,B) {cie_funcLAB2XYZa2(L,A,B) + A / 500}
#declare cie_funcLAB2XYZa3 = function(L,A,B) {cie_funcLAB2XYZa2(L,A,B) - B / 200}
#declare cie_funcLAB2XYZb1 = function(L,A,B) {cie_XYZWhiteReference1 * select(pow(cie_funcLAB2XYZa1(L,A,B), 3) - cie_XYZEpsilon, (cie_funcLAB2XYZa1(L,A,B) - 16.0 / 116.0) / 7.787, (cie_funcLAB2XYZa1(L,A,B) - 16.0 / 116.0) / 7.787, pow(cie_funcLAB2XYZa1(L,A,B), 3))}
#declare cie_funcLAB2XYZb2 = function(L,A,B) {cie_XYZWhiteReference2 * select(L - cie_XYZKappa * cie_XYZEpsilon, L / cie_XYZKappa, L / cie_XYZKappa, pow(((L + 16) / 116), 3))}
#declare cie_funcLAB2XYZb3 = function(L,A,B) {cie_XYZWhiteReference3 * select(pow(cie_funcLAB2XYZa3(L,A,B), 3) - cie_XYZEpsilon, (cie_funcLAB2XYZa3(L,A,B) - 16.0 / 116.0) / 7.787, (cie_funcLAB2XYZa3(L,A,B) - 16.0 / 116.0) / 7.787, pow(cie_funcLAB2XYZa3(L,A,B), 3))}


// input L = between 0 and 100
// input U = between -128 and +128
// input V = between -128 and +128
// output X: between 0 and 100
// output Y: between 0 and 100
// output Z: between 0 and 100
#declare cie_funcLUV2XYZu1 = function(L,U,V) {U/13/L + cie_uPrime}
#declare cie_funcLUV2XYZv1 = function(L,U,V) {V/13/L + cie_vPrime}
#declare cie_funcLUV2XYZa2 = function(L,U,V) {cie_XYZWhiteReference2 * select(8 - L, pow((L + 16)/116, 3), L * pow(3/29, 3))}
#declare cie_funcLUV2XYZa1 = function(L,U,V) {cie_funcLUV2XYZa2(L,U,V) * 9/4 * cie_funcLUV2XYZu1(L,U,V)/cie_funcLUV2XYZv1(L,U,V)}
#declare cie_funcLUV2XYZa3 = function(L,U,V) {cie_funcLUV2XYZa2(L,U,V) * (12 - 3 * cie_funcLUV2XYZu1(L,U,V) - 20 * cie_funcLUV2XYZv1(L,U,V))/4/cie_funcLUV2XYZv1(L,U,V)}


// input X: between 0 and 100
// input Y: between 0 and 100
// input Z: between 0 and 100
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
// Note that out-of-range colors are *NOT* corrected. You will have to do this yourself.
#declare cie_funcXYZ2RGBa1 = function(X,Y,Z) {X/100 *  3.2406 + Y/100 * -1.5372 + Z/100 * -0.4986}
#declare cie_funcXYZ2RGBa2 = function(X,Y,Z) {X/100 * -0.9689 + Y/100 *  1.8758 + Z/100 *  0.0415}
#declare cie_funcXYZ2RGBa3 = function(X,Y,Z) {X/100 *  0.0557 + Y/100 * -0.2040 + Z/100 *  1.0570}
#declare cie_funcXYZ2RGBb1 = function(X,Y,Z) {select(cie_funcXYZ2RGBa1(X,Y,Z) - 0.0031308, 12.92 * cie_funcXYZ2RGBa1(X,Y,Z), 12.92 * cie_funcXYZ2RGBa1(X,Y,Z), 1.055 * pow(cie_funcXYZ2RGBa1(X,Y,Z), 1 / 2.4) - 0.055)}
#declare cie_funcXYZ2RGBb2 = function(X,Y,Z) {select(cie_funcXYZ2RGBa2(X,Y,Z) - 0.0031308, 12.92 * cie_funcXYZ2RGBa2(X,Y,Z), 12.92 * cie_funcXYZ2RGBa2(X,Y,Z), 1.055 * pow(cie_funcXYZ2RGBa2(X,Y,Z), 1 / 2.4) - 0.055)}
#declare cie_funcXYZ2RGBb3 = function(X,Y,Z) {select(cie_funcXYZ2RGBa3(X,Y,Z) - 0.0031308, 12.92 * cie_funcXYZ2RGBa3(X,Y,Z), 12.92 * cie_funcXYZ2RGBa3(X,Y,Z), 1.055 * pow(cie_funcXYZ2RGBa3(X,Y,Z), 1 / 2.4) - 0.055)}


// input R: between 0 and 1
// input G: between 0 and 1
// input B: between 0 and 1
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
// only uses these when converting individual values, not inside isosurfaces or cie_function patterns
#declare cie_correctRGBa1 = function(R,G,B) {cie_selectRGBa(R,cie_selectRGBa(G,cie_selectRGBa(B,R)))}
#declare cie_correctRGBa2 = function(R,G,B) {cie_selectRGBa(R,cie_selectRGBa(G,cie_selectRGBa(B,G)))}
#declare cie_correctRGBa3 = function(R,G,B) {cie_selectRGBa(R,cie_selectRGBa(G,cie_selectRGBa(B,B)))}
#declare cie_correctRGBb1 = function(R,G,B) {cie_selectRGBb(R,cie_selectRGBb(G,cie_selectRGBb(B,cie_correctRGBa1(R,G,B))))}
#declare cie_correctRGBb2 = function(R,G,B) {cie_selectRGBb(R,cie_selectRGBb(G,cie_selectRGBb(B,cie_correctRGBa2(R,G,B))))}
#declare cie_correctRGBb3 = function(R,G,B) {cie_selectRGBb(R,cie_selectRGBb(G,cie_selectRGBb(B,cie_correctRGBa3(R,G,B))))}


// input L = between 0 and 100
// input A = between -128 and +128
// input B = between -128 and +128
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare cie_funcLAB2RGBa1 = function(L,A,B) {cie_funcXYZ2RGBb1(cie_funcLAB2XYZb1(L,A,B),cie_funcLAB2XYZb2(L,A,B),cie_funcLAB2XYZb3(L,A,B))}
#declare cie_funcLAB2RGBa2 = function(L,A,B) {cie_funcXYZ2RGBb2(cie_funcLAB2XYZb1(L,A,B),cie_funcLAB2XYZb2(L,A,B),cie_funcLAB2XYZb3(L,A,B))}
#declare cie_funcLAB2RGBa3 = function(L,A,B) {cie_funcXYZ2RGBb3(cie_funcLAB2XYZb1(L,A,B),cie_funcLAB2XYZb2(L,A,B),cie_funcLAB2XYZb3(L,A,B))}


// input L = between 0 and 100
// input U = between -128 and +128
// input V = between -128 and +128
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare cie_funcLUV2RGBa1 = function(L,U,V) {cie_funcXYZ2RGBb1(cie_funcLUV2XYZa1(L,U,V),cie_funcLUV2XYZa2(L,U,V),cie_funcLUV2XYZa3(L,U,V))}
#declare cie_funcLUV2RGBa2 = function(L,U,V) {cie_funcXYZ2RGBb2(cie_funcLUV2XYZa1(L,U,V),cie_funcLUV2XYZa2(L,U,V),cie_funcLUV2XYZa3(L,U,V))}
#declare cie_funcLUV2RGBa3 = function(L,U,V) {cie_funcXYZ2RGBb3(cie_funcLUV2XYZa1(L,U,V),cie_funcLUV2XYZa2(L,U,V),cie_funcLUV2XYZa3(L,U,V))}


// input L = between 0 and 100
// input C = between 0 and 128
// input H = between 0 and 360
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare cie_funcLCH2RGBa1 = function(L,C,H) {cie_funcXYZ2RGBb1(cie_funcLAB2XYZb1(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)),cie_funcLAB2XYZb2(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)),cie_funcLAB2XYZb3(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)))}
#declare cie_funcLCH2RGBa2 = function(L,C,H) {cie_funcXYZ2RGBb2(cie_funcLAB2XYZb1(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)),cie_funcLAB2XYZb2(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)),cie_funcLAB2XYZb3(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)))}
#declare cie_funcLCH2RGBa3 = function(L,C,H) {cie_funcXYZ2RGBb3(cie_funcLAB2XYZb1(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)),cie_funcLAB2XYZb2(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)),cie_funcLAB2XYZb3(cie_funcLCH2LABa1(L,C,H),cie_funcLCH2LABa2(L,C,H),cie_funcLCH2LABa3(L,C,H)))}
#declare cie_funcLCH2RGBb1 = function(L,C,H) {cie_correctRGBb1(cie_funcLCH2RGBa1(L,C,H),cie_funcLCH2RGBa2(L,C,H),cie_funcLCH2RGBa3(L,C,H))}
#declare cie_funcLCH2RGBb2 = function(L,C,H) {cie_correctRGBb2(cie_funcLCH2RGBa1(L,C,H),cie_funcLCH2RGBa2(L,C,H),cie_funcLCH2RGBa3(L,C,H))}
#declare cie_funcLCH2RGBb3 = function(L,C,H) {cie_correctRGBb3(cie_funcLCH2RGBa1(L,C,H),cie_funcLCH2RGBa2(L,C,H),cie_funcLCH2RGBa3(L,C,H))}


// input H = between 0 and 360
// input C = between 0 and 200
// input L = between 0 and 100
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare cie_funcHCL2RGBa1 = function(H,C,L) {cie_funcXYZ2RGBb1(cie_funcLUV2XYZa1(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)),cie_funcLUV2XYZa2(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)),cie_funcLUV2XYZa3(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)))}
#declare cie_funcHCL2RGBa2 = function(H,C,L) {cie_funcXYZ2RGBb2(cie_funcLUV2XYZa1(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)),cie_funcLUV2XYZa2(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)),cie_funcLUV2XYZa3(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)))}
#declare cie_funcHCL2RGBa3 = function(H,C,L) {cie_funcXYZ2RGBb3(cie_funcLUV2XYZa1(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)),cie_funcLUV2XYZa2(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)),cie_funcLUV2XYZa3(cie_funcHCL2LUVa1(H,C,L),cie_funcHCL2LUVa2(H,C,L),cie_funcHCL2LUVa3(H,C,L)))}
#declare cie_funcHCL2RGBb1 = function(H,C,L) {cie_correctRGBb1(cie_funcHCL2RGBa1(H,C,L),cie_funcHCL2RGBa2(H,C,L),cie_funcHCL2RGBa3(H,C,L))}
#declare cie_funcHCL2RGBb2 = function(H,C,L) {cie_correctRGBb2(cie_funcHCL2RGBa1(H,C,L),cie_funcHCL2RGBa2(H,C,L),cie_funcHCL2RGBa3(H,C,L))}
#declare cie_funcHCL2RGBb3 = function(H,C,L) {cie_correctRGBb3(cie_funcHCL2RGBa1(H,C,L),cie_funcHCL2RGBa2(H,C,L),cie_funcHCL2RGBa3(H,C,L))}


// testing
#local VLCH = <25, 50, 270>;
#local VLAB = <cie_funcLCH2LABa1(VLCH.x,VLCH.y,VLCH.z), cie_funcLCH2LABa2(VLCH.x,VLCH.y,VLCH.z), cie_funcLCH2LABa3(VLCH.x,VLCH.y,VLCH.z)>;
#local VXYZ = <cie_funcLAB2XYZb1(VLAB.x,VLAB.y,VLAB.z), cie_funcLAB2XYZb2(VLAB.x,VLAB.y,VLAB.z), cie_funcLAB2XYZb3(VLAB.x,VLAB.y,VLAB.z)>;
#local VRGB = <cie_funcXYZ2RGBb1(VXYZ.x,VXYZ.y,VXYZ.z), cie_funcXYZ2RGBb2(VXYZ.x,VXYZ.y,VXYZ.z), cie_funcXYZ2RGBb3(VXYZ.x,VXYZ.y,VXYZ.z)>;
#local VCOR = < cie_correctRGBb1(VRGB.x,VRGB.y,VRGB.z),  cie_correctRGBb2(VRGB.x,VRGB.y,VRGB.z),  cie_correctRGBb3(VRGB.x,VRGB.y,VRGB.z)>;
#debug concat("\nVLCH = ", vstr(3, VLCH, ", ", 0, 5), "\n")
#debug concat("\nVLAB = ", vstr(3, VLAB, ", ", 0, 5), "\n")
#debug concat("\nVXYZ = ", vstr(3, VXYZ, ", ", 0, 5), "\n")
#debug concat("\nVRGB = ", vstr(3, VRGB, ", ", 0, 5), "\n")
#debug concat("\nVCOR = ", vstr(3, VCOR, ", ", 0, 5), "\n")