colr 0.3.1

A general purpose, extensible color type unifying color models and their operations at the type level.
Documentation
//! Perceptual color difference (delta E) metrics.

use crate::Color;
use crate::illuminant::Illuminant;
use crate::model::{Lab, Oklab};

#[inline(always)]
fn euclidean(l1: f32, a1: f32, b1: f32, l2: f32, a2: f32, b2: f32) -> f32 {
    let dl = l1 - l2;
    let da = a1 - a2;
    let db = b1 - b2;
    (dl * dl + da * da + db * db).sqrt()
}

// 25^7 used in the CIEDE2000 chroma weighting factor.
const CIEDE2000_C25_7: f32 = 6103515625.0;

// Angle constants from CIE 142-2001, in radians.
const DEG_6: f32 = core::f32::consts::PI / 30.0;
const DEG_25: f32 = 25.0 * core::f32::consts::PI / 180.0;
const DEG_60: f32 = core::f32::consts::FRAC_PI_3;
const DEG_63: f32 = 63.0 * core::f32::consts::PI / 180.0;
const DEG_275: f32 = 275.0 * core::f32::consts::PI / 180.0;

#[allow(clippy::excessive_precision)]
fn ciede2000(l1: f32, a1: f32, b1: f32, l2: f32, a2: f32, b2: f32) -> f32 {
    let c1 = (a1 * a1 + b1 * b1).sqrt();
    let c2 = (a2 * a2 + b2 * b2).sqrt();
    let c_bar = (c1 + c2) * 0.5;
    let c_bar7 = c_bar * c_bar * c_bar * c_bar * c_bar * c_bar * c_bar;
    let g = 0.5 * (1.0 - (c_bar7 / (c_bar7 + CIEDE2000_C25_7))).sqrt();
    let a1p = a1 * (1.0 + g);
    let a2p = a2 * (1.0 + g);

    let c1p = (a1p * a1p + b1 * b1).sqrt();
    let c2p = (a2p * a2p + b2 * b2).sqrt();
    let h1p = (b1).atan2(a1p).rem_euclid(core::f32::consts::TAU);
    let h2p = (b2).atan2(a2p).rem_euclid(core::f32::consts::TAU);

    let dl_p = l2 - l1;
    let dc_p = c2p - c1p;
    let dh_p = if c1p * c2p == 0.0 {
        0.0
    } else {
        (h2p - h1p + core::f32::consts::PI).rem_euclid(core::f32::consts::TAU)
            - core::f32::consts::PI
    };
    let dh_big = 2.0 * (c1p * c2p).sqrt() * (dh_p * 0.5).sin();

    let l_bar = (l1 + l2) * 0.5;
    let c_barp = (c1p + c2p) * 0.5;
    let h_barp = if c1p * c2p == 0.0 {
        h1p + h2p
    } else if (h1p - h2p).abs() <= core::f32::consts::PI {
        (h1p + h2p) * 0.5
    } else if h1p + h2p < core::f32::consts::TAU {
        (h1p + h2p + core::f32::consts::TAU) * 0.5
    } else {
        (h1p + h2p - core::f32::consts::TAU) * 0.5
    };

    let t = 1.0
        - 0.17 * (h_barp - core::f32::consts::FRAC_PI_6).cos()
        + 0.24 * (2.0 * h_barp).cos()
        + 0.32 * (3.0 * h_barp + DEG_6).cos()
        - 0.20 * (4.0 * h_barp - DEG_63).cos();

    let l_bar_shifted = l_bar - 50.0;
    let sl = 1.0 + 0.015 * l_bar_shifted * l_bar_shifted
        / (20.0 + l_bar_shifted * l_bar_shifted).sqrt();
    let sc = 1.0 + 0.045 * c_barp;
    let sh = 1.0 + 0.015 * c_barp * t;

    let c_barp7 = { let c2 = c_barp * c_barp; let c4 = c2 * c2; c4 * c2 * c_barp };
    let rc = 2.0 * (c_barp7 / (c_barp7 + CIEDE2000_C25_7)).sqrt();
    let d_theta = DEG_60
        * (-((h_barp - DEG_275) / DEG_25).powi(2)).exp();
    let rt = -(2.0 * d_theta).sin() * rc;

    let term_l = dl_p / sl;
    let term_c = dc_p / sc;
    let term_h = dh_big / sh;

    (term_l * term_l + term_c * term_c + term_h * term_h + rt * term_c * term_h).sqrt()
}

macro_rules! impl_delta_e_lab {
    ($n:literal, $bound:tt) => {
        impl<$bound: Illuminant> Color<[f32; $n], Lab<$bound>> {
            /// CIE76 color difference. Euclidean distance in L*a*b*.
            ///
            /// Fast but overestimates differences in saturated regions.
            /// Prefer `delta_e2000` when perceptual accuracy matters.
            #[inline(always)]
            pub fn delta_e76(self, other: Self) -> f32 {
                let [l1, a1, b1, ..] = self.inner();
                let [l2, a2, b2, ..] = other.inner();
                euclidean(l1, a1, b1, l2, a2, b2)
            }

            /// CIEDE2000 color difference.
            ///
            /// Reference: CIE 142-2001, Luo, Cui, Rigg (2001).
            /// Parametric factors kL, kC, kH are all 1 (graphic arts default).
            pub fn delta_e2000(self, other: Self) -> f32 {
                let [l1, a1, b1, ..] = self.inner();
                let [l2, a2, b2, ..] = other.inner();
                ciede2000(l1, a1, b1, l2, a2, b2)
            }
        }
    };
}

macro_rules! impl_delta_e_oklab {
    ($n:literal) => {
        impl Color<[f32; $n], Oklab> {
            /// Oklab color difference. Euclidean distance in Oklab.
            ///
            /// Oklab is designed to be perceptually uniform, so this
            /// metric is effective for most use cases without the
            /// complexity of CIEDE2000.
            #[inline(always)]
            pub fn delta_e_ok(self, other: Self) -> f32 {
                let [l1, a1, b1, ..] = self.inner();
                let [l2, a2, b2, ..] = other.inner();
                euclidean(l1, a1, b1, l2, a2, b2)
            }
        }
    };
}

impl_delta_e_lab!(3, W);
impl_delta_e_lab!(4, W);
impl_delta_e_oklab!(3);
impl_delta_e_oklab!(4);