CIEDE2000

CIEDE2000.cs

// Implementation of "The CIEDE2000 Color-Difference Formula: Implementation Notes, Supplementary Test Data, and Mathematical Observations".
public static float CIEDE2000(Lab lab1, Lab lab2)
{
    float Pi = Mathf.PI;
    float Pi2 = 2.0f * Mathf.PI;

    float kL = 1.0f, kC = 1.0f, kH = 1.0f;

    var mCs = (Mathf.Sqrt( lab1.a * lab1.a + lab1.b * lab1.b) + Mathf.Sqrt(lab2.a * lab2.a + lab2.b * lab2.b)) / 2.0f;
    var G = 0.5f * (1.0f - Mathf.Sqrt(Mathf.Pow(mCs, 7) / (Mathf.Pow(mCs, 7) + Mathf.Pow(25.0f, 7))));
    var a1p = (1.0f + G) * lab1.a;
    var a2p = (1.0f + G) * lab2.a;
    var C1p = Mathf.Sqrt(a1p * a1p + lab1.b * lab1.b);
    var C2p = Mathf.Sqrt(a2p * a2p + lab2.b * lab2.b);

    var h1p = Mathf.Abs(a1p) + Mathf.Abs(lab1.b) > double.Epsilon ? Mathf.Atan2(lab1.b, a1p) : 0.0f;
    if (h1p < 0.0) h1p += Pi2;
    var h2p = Mathf.Abs(a2p) + Mathf.Abs(lab2.b) > double.Epsilon ? Mathf.Atan2(lab2.b, a2p) : 0.0f;
    if (h2p < 0.0) h2p += Pi2;

    var dLp = lab2.L - lab1.L;
    var dCp = C2p - C1p;

    var dhp = 0.0f;
    var cProdAbs = Mathf.Abs(C1p * C2p);
    if (cProdAbs > Mathf.Epsilon && Mathf.Abs(h1p - h2p) <= Pi)
    {
        dhp = h2p - h1p;
    }
    else if (cProdAbs > Mathf.Epsilon && h2p - h1p > Pi)
    {
        dhp = h2p - h1p - Pi2;
    }
    else if (cProdAbs > Mathf.Epsilon && h2p - h1p < -Pi)
    {
        dhp = h2p - h1p + Pi2;
    }

    var dHp = 2.0f * Mathf.Sqrt(C1p * C2p) * Mathf.Sin(dhp / 2.0f);

    var mLp =  (lab1.L + lab2.L) / 2.0f;
    var mCp = (C1p + C2p) / 2.0f;

    var mhp = 0.0f;
    if (cProdAbs > Mathf.Epsilon && Mathf.Abs(h1p - h2p) <= Pi)
    {
        mhp = (h1p + h2p) / 2.0f;
    }
    else if (cProdAbs > Mathf.Epsilon && Mathf.Abs(h1p - h2p) > Pi && h1p + h2p < Pi2)
    {
        mhp = (h1p + h2p + Pi2) / 2.0f;
    }
    else if (cProdAbs > Mathf.Epsilon && Mathf.Abs(h1p - h2p) > Pi && h1p + h2p >= Pi2)
    {
        mhp = (h1p + h2p - Pi2) / 2.0f;
    }
    else if (cProdAbs <= Mathf.Epsilon)
    {
        mhp = h1p + h2p;
    }

    var T = 1.0f - 0.17f * Mathf.Cos(mhp - Pi / 6.0f) + .24f * Mathf.Cos(2.0f * mhp) +
        0.32f * Mathf.Cos(3.0f * mhp + Pi / 30.0f) - 0.2f * Mathf.Cos(4.0f * mhp - 7.0f * Pi / 20.0f);
    var dTheta = Pi / 6.0f * Mathf.Exp(-Mathf.Pow((mhp / (2.0f * Pi) * 360.0f - 275.0f) / 25.0f, 2));
    var RC = 2.0f * Mathf.Sqrt(Mathf.Pow(mCp, 7) / (Mathf.Pow(mCp, 7) + Mathf.Pow(25.0f, 7)));
    var mlpSqr = (mLp - 50.0f) * (mLp - 50.0f);
    var SL = 1.0f + 0.015f * mlpSqr / Mathf.Sqrt(20.0f + mlpSqr);
    var SC = 1.0f + 0.045f * mCp;
    var SH = 1.0f + 0.015f * mCp * T;
    var RT = -Mathf.Sin(2.0f * dTheta) * RC;

    var de00 = Mathf.Sqrt(
        Mathf.Pow(dLp / (kL * SL), 2) + Mathf.Pow(dCp / (kC * SC), 2) + Mathf.Pow(dHp / (kH * SH), 2) +
        RT * dCp / (kC * SC) * dHp / (kH * SH)
    );
    return de00;
}