color_palette_conversion.rs

use num_traits::cast::ToPrimitive;
use num_traits::{abs, Pow, Zero};
use num_traits::real::Real;

// Formula from https://en.wikipedia.org/wiki/SRGB#From_sRGB_to_CIE_XYZ
pub fn srgb_to_cie_xyz(r: u8, g: u8, b: u8) -> (f64, f64, f64) {
    fn value_to_linear(value: &f64) -> f64 {
        if *value < 0.04045 {
            *value / 12.92
        } else {
            ((*value + 0.055) / 1.055).pow(2.4)
        }
    }

    let r_norm: f64 = r.to_f64().unwrap() / 255.0;
    let g_norm: f64 = g.to_f64().unwrap() / 255.0;
    let b_norm: f64 = b.to_f64().unwrap() / 255.0;

    let r_lin = value_to_linear(&r_norm);
    let g_lin = value_to_linear(&g_norm);
    let b_lin = value_to_linear(&b_norm);

    let x = 0.4124 * r_lin + 0.3576 * g_lin + 0.1805 * b_lin;
    let y = 0.2126 * r_lin + 0.7152 * g_lin + 0.0722 * b_lin;
    let z = 0.0193 * r_lin + 0.1192 * g_lin + 0.9505 * b_lin;

    return (x, y, z);
}

// Formula from https://en.wikipedia.org/wiki/CIELAB_color_space#From_CIEXYZ_to_CIELAB
pub fn cie_xyz_to_cie_lab(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    fn f(t: f64) -> f64 {
        // Constants:
        //   delta^3 = (6/29)^3 = 0.008856
        //   (1/3)delta^(-2) = 7.787037
        //   4/29 = 0.137931
        if t > 0.008856 {
            t.pow(0.333333)
        } else {
            t * 7.787037 + 0.137931
        }
    }

    // Divide by reference luminant, as specified in Standard Illuminant D65
    // https://en.wikipedia.org/wiki/Illuminant_D65
    // (default for midday daylight)
    let x_norm = x / 0.950489;
    let y_norm = y; // / 1.0
    let z_norm = z / 1.088840;

    let fx = f(x_norm);
    let fy = f(y_norm);
    let fz = f(z_norm);

    let l = 116. * fy - 16.;
    let a = 500. * (fx - fy);
    let b = 200. * (fy - fz);

    return (l, a, b);
}

// Formula from https://en.wikipedia.org/wiki/Color_difference#CIEDE2000
pub fn cie_de2000(lab1: (f64, f64, f64), lab2: (f64, f64, f64)) -> f64 {
    fn deg2rad(deg: f64) -> f64 {
        deg * std::f64::consts::PI / 180.
    }

    fn rad2deg(rad: f64) -> f64 {
        rad * 180. / std::f64::consts::PI
    }

    fn wrap_angle(angle: f64) -> f64 {
        ((angle % 360.) + 360.) % 360.
    }

    let (l1, a1, b1) = lab1;
    let (l2, a2, b2) = lab2;

    let delta_l = l2 - l1;
    let l_mean = (l1 + l2) / 2.;

    let c1 = ((a1.pow(2) + b1.pow(2)) as f64).sqrt();
    let c2 = ((a2.pow(2) + b2.pow(2)) as f64).sqrt();
    let c_mean = (c1 + c2) / 2.;

    // 6103515625 == 25^7
    let c_mean_7 :f64 = c_mean.pow(7);
    let constant1 = 1. + (1. - (c_mean_7 / (c_mean_7 + 6103515625.)).sqrt()) / 2.;
    let a1_prime = a1 * constant1;
    let a2_prime = a2 * constant1;

    let c1_prime = ((a1_prime.pow(2) + b1.pow(2)) as f64).sqrt();
    let c2_prime = ((a2_prime.pow(2) + b2.pow(2)) as f64).sqrt();

    let c_prime_mean = (c1_prime + c2_prime) / 2.;
    let delta_c_prime = c2_prime - c1_prime;

    let h1 = wrap_angle(rad2deg(b1.atan2(a1_prime)));
    let h2 = wrap_angle(rad2deg(b2.atan2(a2_prime)));

    let (delta_h_big, h_big_mean) = if c1_prime == 0. || c2_prime == 0. {
        (0., h1 + h2)
    } else {
        let delta_h = if abs(h1 - h2) <= 180. {
            h2 - h1
        } else if h2 <= h1 {
            h2 - h1 + 360.
        } else {
            h2 - h1 - 360.
        };

        let d_h_big = 2. * (c1_prime * c2_prime).sqrt() * rad2deg((delta_h / 2.).sin());
        let h_b = if abs(h1 - h2) <= 180. {
            (h1 + h2) / 2.
        } else if h1 + h2 < 360. {
            (h1 + h2 + 360.) / 2.
        } else {
            (h1 + h2 - 360.) / 2.
        };

        (d_h_big, h_b)
    };

    let t = 1. - 0.17 * rad2deg((h_big_mean - 30.).cos())
        + 0.24 * rad2deg((2. * h_big_mean).cos())
        + 0.32 * rad2deg((3. * h_big_mean + 6.).cos())
        - 0.20 * rad2deg((4. * h_big_mean - 63.).cos());

    let s_l = 1. + ((0.015 * (l_mean - 50.).pow(2)) / ((20. + (l_mean - 50.).pow(2)) as f64).sqrt());
    let s_c = 1. + (0.045 * c_prime_mean);
    let s_h = 1. * (0.015 * c_prime_mean * t);

    let c_prime_mean_seven: f64 = c_prime_mean.pow(7);
    let r_t = -2.
        * ((c_prime_mean_seven / (c_prime_mean_seven + 25.0.pow(7))) as f64).sqrt()
        * rad2deg((60. * ((-1. * ((h_big_mean - 275.) / 25.).pow(2)) as f64).exp()).sin());

    let k_c = 1.;
    let k_l = 1.;
    let k_h = 1.;

    let ratio_l: f64 = delta_l / (k_l * s_l);
    let ratio_c: f64  = delta_c_prime / (k_c * s_c);
    let ratio_h: f64  = delta_h_big / (k_h * s_h);

    let delta_e: f64 = ((ratio_l.pow(2) + ratio_c.pow(2) + ratio_h.pow(2) + (r_t * ratio_c * ratio_h)) as f64).sqrt();

    return delta_e;
}