RGB2HUSL

husl.mac

/* References: */
/* http://www.brucelindbloom.com/ */
/* https://github.com/jrus/chromatist */
/* https://github.com/boronine/husl */

/* D65 white point */
refX: 0.95047;
refY: 1.0;
refZ: 1.08883;

labE: (6 / 29) ^ 3;
labK: (29 / 3) ^ 3;
refU: (4 * refX) / (refX + (15 * refY) + (3 * refZ));
refV: (9 * refY) / (refX + (15 * refY) + (3 * refZ));

/* Source: */
/* http://www.brucelindbloom.com/index.html?WorkingSpaceInfo.html */

x_r: 0.64;
y_r: 0.33;
x_g: 0.30;
y_g: 0.60;
x_b: 0.15;
y_b: 0.06;

X_r: x_r / y_r;
Y_r: 1;
Z_r: (1 - x_r - y_r) / y_r;

X_g: x_g / y_g;
Y_g: 1;
Z_g: (1 - x_g - y_g) / y_g;

X_b: x_b / y_b;
Y_b: 1;
Z_b: (1 - x_b - y_b) / y_b;

M_P: matrix([X_r, X_g, X_b],
            [Y_r, Y_g, Y_b],
            [Z_r, Z_g, Z_b]);

M_S: invert(M_P) . [refX, refY, refZ];

/* sRGB D65 matrices */
M_RGB_XYZ: M_P * transpose(addcol(M_S, M_S, M_S));
M_XYZ_RGB: invert(M_RGB_XYZ);


/* Gamma Companding */
/* linearValue(v) := v^2.2; */
/* sRGB Companding */
linearValue(v) := (if (v <= 0.04045) then v / 12.92 else ((v + 0.055) / 1.055)^2.4);
/* nonLinearValue(v) := v^(1 / 2.2); */
nonLinearValue(v) := v;
/* nonLinearValue(v) := (if (v <= 0.0031308) then 12.92 * v else 1.055 * (v ^ (1 / 2.4)) - 0.055); */

LUV2LCH(L, U, V) := block([C, H],
                          C: sqrt(U^2 + V^2),
                          H: 180 * atan2(V, U) / %pi,
                          [L, C, H]);

LCH2LUV(L, C, H) := block([H_Rad, U, V],
                          H_Rad: %pi * H / 180,
                          U : C * cos(H_Rad),
                          V : C * sin(H_Rad),
                          [L, U, V]);

XYZ2LUV(X, Y, Z) := block([y, u, v, L, U, V],
                          y: Y / refY,
                          L: (if (y <= labE) then (labK * y) else (116 * y^(1/3) - 16)),
                          u: 4 * X / (X + 15 * Y + 3 * Z),
                          v: 9 * Y / (X + 15 * Y + 3 * Z),
                          U: 13 * L * (u - refU),
                          V: 13 * L * (v - refV),
                          [L, U, V]);

LUV2XYZ(L, U, V) := block([u, v, X, Y, Z],
                          u: refU + U / (13 * L),
                          v: refV + V / (13 * L),
                          Y: refY * (if (L > (labK * labE)) then ((L + 16) / 116)^3 else L / labK),
                          X: Y * (9 * u) / (4 * v),
                          Z: Y * (12 - 3 * u - 20 * v) / (4 * v),
                          [X, Y, Z]);

RGB2XYZ(R, G, B) := M_RGB_XYZ . map(linearValue, [R, G, B]);

XYZ2RGB(X, Y, Z) := matrixmap(nonLinearValue, M_XYZ_RGB . [X, Y, Z]);


RGB2LCH(R, G, B) := block([XYZ, LUV],
                          XYZ: RGB2XYZ(R, G, B),
                          LUV: XYZ2LUV(XYZ[1][1], XYZ[2][1], XYZ[3][1]),
                          LUV2LCH(LUV[1], LUV[2], LUV[3]));

LCH2RGB(L, C, H) := block([LUV, XYZ],
                          LUV: LCH2LUV(L, C, H),
                          XYZ: LUV2XYZ(LUV[1], LUV[2], LUV[3]),
                          XYZ2RGB(XYZ[1], XYZ[2], XYZ[3]));

limitChromas(L, H) := block([M, i, j, results],
                            results: [],
                            M: LCH2RGB(L, C, H),
                            for i: 1 while i < 4 do
                               for j: 0 while j < 2 do
                                   results: cons(rhs(solve(M[i][1] = j, C)[1]), results),
                            reverse(results));

positivep(x) := x > 0; 

RGB2HUSL(R, G, B) := block([LCH, chromas, maxC],
                           LCH : RGB2LCH(R, G, B),
                           chromas: limitChromas(LCH[1], LCH[3]),
                           maxChroma: lmin(sublist(chromas, positivep)),
                           [LCH[3], 100 * LCH[2] / maxChroma, LCH[1]]);

ratprint: false;                                 

/* getGNUPlotPreamble(L) := printf(false, "set title 'L = ~d'; */
/*                                         set size square; */
/*                                         set xtics 60; */
/*                                         set ytics 60; */
/*                                         set style line 1 lc rgb '#FE8980'; */
/*                                         set style line 2 lc rgb '#5EBE7C'; */
/*                                         set style line 3 lc rgb '#59B0FF'", L); */

/* Generates a plot for L = 50 */

/* block([L], */
/*       L: 50, */
/*       plot2d(limitChromas(L, x), */
/*              [x, 0, 360], */
/*              [y, 0, 360], */
/*              [legend, false], */
/*              [style, [lines, 1.5], [lines, 3]], */
/*              [color, red, red, green, green, blue, blue], */
/*              [gnuplot_preamble, getGNUPlotPreamble(L)], */
/*              [xlabel, "H"], */
/*              [ylabel, "C"], */
/*              [plot_format, gnuplot])); */

/* Creates 99 PNG files for L=1 to 99 */
/* A GIF animation can be created afterwards with: */
/*     gm convert plot??.png GIF:- | \ */
/*     gifsicle --delay 10 --colors 32 --multifile --loop - > plot.gif */

/* for L:1 thru 99 step 1 do */
/*     plot2d(limitChromas(L, x), */
/*            [x, 0, 360], */
/*            [y, 0, 360], */
/*            [legend, false], */
/*            [style, [lines, 1.5], [lines, 3]], */
/*            [color, red, red, green, green, blue, blue], */
/*            [gnuplot_preamble, getGNUPlotPreamble(L)], */
/*            [xlabel, "H"], */
/*            [ylabel, "C"], */
/*            [gnuplot_term, "pngcairo size 500, 500"], */
/*            [gnuplot_out_file, sconcat("plot", printf(false, "~2,'0d", L), ".png")], */
/*            [plot_format, gnuplot]); */

/* vim: set ft=maxima: */