Quantify color difference using CIE94 and CIEDE2000 algorithms

// input is 16 bit sRGB
function RGB2LAB(STRUCT RGBcolor &rgb, STRUCT LABcolor &Lab)
    // X, Y and Z output refer to a D65/2° standard illuminant.
    Make /free w={rgb.red, rgb.green, rgb.blue}
    w /= 0xFFFF
    w = w > 0.04045 ? ((w + 0.055) / 1.055)^2.4 : w / 12.92
    w *= 100
    Make /free/N=3 xyz
    xyz[0] = w[0] * 0.4124 + w[1] * 0.3576 + w[2] * 0.1805
    xyz[1] = w[0] * 0.2126 + w[1] * 0.7152 + w[2] * 0.0722
    xyz[2] = w[0] * 0.0193 + w[1] * 0.1192 + w[2] * 0.9505
    Make /free XYZref={95.047,100,108.883}
    w = xyz / XYZref
    w = (w > 0.008856) ? w^(1/3) : 7.787 * w + (16 / 116)
    Lab.L = (116 * w[1]) - 16
    Lab.a = 500 * (w[0] - w[1])
    Lab.b = 200 * (w[1] - w[2])
end

function deltaE94(STRUCT LABcolor &Lab1, STRUCT LABcolor &Lab2)
    variable kL = 1, K1 = 0.045, K2 = 0.015
    variable C1 = sqrt(Lab1.a^2 + Lab1.b^2)
    variable C2 = sqrt(Lab2.a^2 + Lab2.b^2)
    variable deltaL = Lab2.L - Lab1.L
    variable deltaC = C2 - C1
    variable deltaH = (Lab2.a - Lab1.a)^2 + (Lab2.b-Lab1.b)^2 - deltaC^2
    deltaH = deltaH > 0 ? sqrt(deltaH) : 0 
    variable SL = 1
    variable SC = 1 + ( K1 * C1 )
    variable SH = 1 + ( K2 * C1 )
    return sqrt( (deltaL/kL/SL)^2 + (deltaC/SC)^2 + (deltaH/SH)^2 )
end

function deltaE2000(STRUCT LABcolor &Lab1, STRUCT LABcolor &Lab2)
    variable kL = 1, kC = 1, kH = 1 // Weight factors  
    variable deltaL = Lab2.L - Lab1.L
    variable Lbar = (Lab1.L + Lab2.L) / 2
    variable C1a = sqrt(Lab1.a^2 + Lab1.b^2)
    variable C2a = sqrt(Lab2.a^2 + Lab2.b^2)
    variable Cbar = (C1a + C2a)/2
    variable a1prime = Lab1.a + Lab1.a/2 * (1 - sqrt(Cbar^7/(Cbar^7+25^7)))
    variable a2prime = Lab2.a + Lab2.a/2 * (1 - sqrt(Cbar^7/(Cbar^7+25^7)))
    variable C1prime = sqrt(a1prime^2 + Lab1.b^2)
    variable C2prime = sqrt(a2prime^2 + Lab2.b^2)
    variable Cprimebar = (C1prime + C2prime)/2
    variable deltaCprime = C2prime - C1prime
    variable h1prime = Lab1.b==a1prime ? 0 : atan2(Lab1.b, a1prime)
    h1prime += 2*Pi*(h1prime<0)
    variable h2prime = Lab2.b==a2prime ? 0 : atan2(Lab2.b, a2prime)
    h2prime += 2*Pi*(h2prime<0)
    variable deltaHprime = h2prime-h1prime // small h
    deltaHprime += (abs(h1prime-h2prime)) <= Pi ? 0 : 2*Pi*(1-2*(h2prime>h1prime)) // small h
    deltaHprime = 2*sqrt(C1prime*C2prime)*sin(deltahprime/2) // large H
    variable Hprimebar = (h1prime + h2prime)/2
    Hprimebar += (abs(h1prime-h2prime)) <= Pi ? 0 : Pi*(1-2*((h1prime+h2prime) >= (2*Pi)))
    variable T = 1 - 0.17 * cos(Hprimebar - Pi/6) + 0.24 * cos(2*Hprimebar) + 0.32 * cos(3*Hprimebar + 6*Pi/180) - 0.2 * cos(4*Hprimebar - 63*Pi/180)
    variable SL = 1 + 0.015 * (Lbar-50)^2 / sqrt(20 + (Lbar-50)^2)
    variable SC = 1 + 0.045 * Cprimebar
    variable SH = 1 + 0.015 * Cprimebar * T
    variable RT = -2 * sqrt(Cbar^7/(Cbar^7+25^7)) * sin(Pi/3 * exp(-((Hprimebar-275*Pi/180)/(25*Pi/180))^2))
    return sqrt( (deltaL/SL/kL)^2 + (deltaCprime/SC/kC)^2 + (deltaHprime/SH/kH)^2 + RT * deltaCprime/SC/kC * deltaHprime/SH/kH )
end

structure LABcolor
    float L, a, b
endstructure