/* Authors: Alexei Boronine, Bastien Dejean */

/* References:
https://gist.github.com/baskerville/8297428
http://www.brucelindbloom.com/ */
 
/* The following values are taken directly from matrix.wxm */

refX: 3127/3290;
refY: 1;
refZ: 3583/3290;

M_RGB_XYZ: matrix(
    [506752/1228815, 87881/245763,  12673/70218],
    [87098/409605, 175762/245763,  12673/175545],
    [7918/409605, 87881/737289, 1001167/1053270]);
M_XYZ_RGB: invert(M_RGB_XYZ);

/* Note that the standard defines these constants as decimals
and is incorrect. We are using the correct non-standard values.
See: http://www.brucelindbloom.com/index.html?LContinuity.html */
epsilon: (6 / 29) ^ 3;
kappa: (29 / 3) ^ 3;
refU: (4 * refX) / (refX + (15 * refY) + (3 * refZ));
refV: (9 * refY) / (refX + (15 * refY) + (3 * refZ));

/* Not used, see linearRGBChannelFromLUV */
fromLinear(c):= if (c<=0.0031308)
    then 12.92 * c
    else 1.055 * c ^ (1 / 2.4) - 0.055;

f_inv(t) := if (t ^ 3 > epsilon) 
    then t ^ 3
    else (116 * t - 16) / kappa;

/* Returns the value of a linear sRGB channel specified by its matrix 
values. Note that this is not the final sRGB value since gamma correction 
is not applied, however, the value happens to be correct in two cases: 
0 and 1. (Solve fromLinear above to see for yourself). As it happens, the 
only use for this function is to solve it for 0 and 1, which means 
fromLinear won't be necessary.
*/

linearRGBChannelFromLUV(L, U, V, m1, m2, m3) := block(
    [varY, varU, varV, Y, X, Z],
    varY : f_inv((L + 16) / 116),
    varU : U / (13 * L) + refU,
    varV : V / (13 * L) + refV,
    Y    : varY * refY,
    X    : 0 - (9 * Y * varU) / ((varU - 4) * varV - varU * varV),
    Z    : (9 * Y - (15 * varV * Y) - (varV * X)) / (3 * varV),
    [X, Y, Z] . [m1, m2, m3]
);

/* Given L, U, RGB channel and limit t (0 or 1), return the V value
that will push the CIE LUV color outside the sRGB gamut for the given
channel and the given limit. */

solve(linearRGBChannelFromLUV(L, U, V, m1, m2, m3) = t, V);

/* Bounding function generator */

bound(L, m1, m2, m3, t) := solve(linearRGBChannelFromLUV(L, x, y, m1, m2, m3) = t, y)[1];

/* Generate bounding function that draws a line on CIE LUV space, stepping
outside of which will push the red RGB value below 0 */

bound(50, M_XYZ_RGB[1][1], M_XYZ_RGB[1][2], M_XYZ_RGB[1][3], 0);

block(
    [L],
    L: 50,
    plot2d([
        rhs(bound(L, M_XYZ_RGB[1][1], M_XYZ_RGB[1][2], M_XYZ_RGB[1][3], 0)),
        rhs(bound(L, M_XYZ_RGB[1][1], M_XYZ_RGB[1][2], M_XYZ_RGB[1][3], 1)),
        rhs(bound(L, M_XYZ_RGB[2][1], M_XYZ_RGB[2][2], M_XYZ_RGB[2][3], 0)),
        rhs(bound(L, M_XYZ_RGB[2][1], M_XYZ_RGB[2][2], M_XYZ_RGB[2][3], 1)),
        rhs(bound(L, M_XYZ_RGB[3][1], M_XYZ_RGB[3][2], M_XYZ_RGB[3][3], 0)),
        rhs(bound(L, M_XYZ_RGB[3][1], M_XYZ_RGB[3][2], M_XYZ_RGB[3][3], 1))
    ], [x,-400,400], [y,-400,400], [
        legend,
        "red = 0",
        "red = 1",
        "green = 0",
        "green = 1",
        "blue = 0",
        "blue = 1"
    ], [
        xlabel,
        "hue"
    ], [
        ylabel,
        "chroma"
    ])
);

fpprec: 17;
bfloat(M_XYZ_RGB);
bfloat(M_RGB_XYZ);
bfloat(refX);
bfloat(refY);
bfloat(refZ);
bfloat(refU);
bfloat(refV);
bfloat(kappa);
bfloat(epsilon);