Implement CIEDE2000 color difference for Lab/Lch

palette/src/color_difference.rs

@@ -0,0 +1,129 @@
+use component::FloatComponent;
+use from_f64;
+
+///A trait for calculating the color difference between two colors.
+pub trait ColorDifference {
+    ///The type of the calculated color difference
+    type Scalar: FloatComponent;
+
+    ///Return the difference or distance between two colors
+    fn get_color_difference(&self, other: &Self) -> Self::Scalar;
+}
+
+/// Container of components necessary to calculate CIEDE color difference
+pub struct LabColorDiff<T: FloatComponent> {
+    /// Lab color lightness
+    pub l: T,
+    /// Lab color a* value
+    pub a: T,
+    /// Lab color b* value
+    pub b: T,
+    /// Lab color chroma value
+    pub chroma: T,
+}
+
+/// Calculate the CIEDE2000 color difference for two colors in Lab color space.
+/// There is a "just noticeable difference" between two colors when the delta E
+/// is roughly greater than 1. Thus, the color difference is more suited for
+/// calculating small distances between colors as opposed to large differences.
+#[rustfmt::skip]
+pub fn get_ciede_difference<T: FloatComponent>(this: &LabColorDiff<T>, other: &LabColorDiff<T>) -> T {
+    let c_bar = (this.chroma + other.chroma) / from_f64(2.0);
+    let c_bar_pow_seven = c_bar * c_bar * c_bar * c_bar * c_bar * c_bar * c_bar;
+    let twenty_five_pow_seven = from_f64(6103515625.0);
+    let pi_over_180 = from_f64::<T>(core::f64::consts::PI / 180.0);
+
+    let g = from_f64::<T>(0.5)
+        * (from_f64::<T>(1.0)
+            - (c_bar_pow_seven / (c_bar_pow_seven + twenty_five_pow_seven)).sqrt());
+    let a_one_prime = this.a * (from_f64::<T>(1.0) + g);
+    let a_two_prime = other.a * (from_f64::<T>(1.0) + g);
+    let c_one_prime = (a_one_prime * a_one_prime + this.b * this.b).sqrt();
+    let c_two_prime = (a_two_prime * a_two_prime + other.b * other.b).sqrt();
+
+    let calc_h_prime = |b: T, a_prime: T| -> T {
+        if b == T::zero() && a_prime == T::zero() {
+            from_f64(0.0)
+        } else {
+            let result = b.atan2(a_prime).to_degrees();
+            if result < T::zero() {
+                result + from_f64(360.0)
+            } else {
+                result
+            }
+        }
+    };
+    let h_one_prime = calc_h_prime(this.b, a_one_prime);
+    let h_two_prime = calc_h_prime(other.b, a_two_prime);
+
+    let h_prime_difference = (h_one_prime - h_two_prime).abs();
+
+    let delta_h_prime: T = if c_one_prime == T::zero() || c_two_prime == T::zero() {
+        from_f64(0.0)
+    } else {
+        if h_prime_difference <= from_f64(180.0) {
+            h_two_prime - h_one_prime
+        } else {
+            if h_two_prime <= h_one_prime {
+                h_two_prime - h_one_prime + from_f64(360.0)
+            } else {
+                h_two_prime - h_one_prime - from_f64(360.0)
+            }
+        }
+    };
+
+    let delta_big_h_prime = from_f64::<T>(2.0)
+        * (c_one_prime * c_two_prime).sqrt()
+        * (delta_h_prime / from_f64(2.0) * pi_over_180).sin();
+    let h_bar_prime = if c_one_prime == T::zero() || c_two_prime == T::zero() {
+        h_one_prime + h_two_prime
+    } else {
+        if h_prime_difference > from_f64(180.0) {
+            (h_one_prime + h_two_prime + from_f64(360.0)) / from_f64(2.0)
+        } else {
+            (h_one_prime + h_two_prime) / from_f64(2.0)
+        }
+    };
+
+    let l_bar = (this.l + other.l) / from_f64(2.0);
+    let c_bar_prime = (c_one_prime + c_two_prime) / from_f64(2.0);
+
+    let t: T = from_f64::<T>(1.0)
+        - from_f64::<T>(0.17) * ((h_bar_prime - from_f64(30.0)) * pi_over_180).cos()
+        + from_f64::<T>(0.24) * ((h_bar_prime * from_f64(2.0)) * pi_over_180).cos()
+        + from_f64::<T>(0.32) * ((h_bar_prime * from_f64(3.0) + from_f64(6.0)) * pi_over_180).cos()
+        - from_f64::<T>(0.20) * ((h_bar_prime * from_f64(4.0) - from_f64(63.0)) * pi_over_180).cos();
+    let s_l = from_f64::<T>(1.0)
+        + ((from_f64::<T>(0.015) * (l_bar - from_f64(50.0)) * (l_bar - from_f64(50.0)))
+            / ((l_bar - from_f64(50.0)) * (l_bar - from_f64(50.0)) + from_f64(20.0)).sqrt());
+    let s_c = from_f64::<T>(1.0) + from_f64::<T>(0.045) * c_bar_prime;
+    let s_h = from_f64::<T>(1.0) + from_f64::<T>(0.015) * c_bar_prime * t;
+
+    let delta_theta = from_f64::<T>(30.0)
+        * (-(((h_bar_prime - from_f64(275.0)) / from_f64(25.0))
+            * ((h_bar_prime - from_f64(275.0)) / from_f64(25.0))))
+        .exp();
+    let c_bar_prime_pow_seven = c_bar_prime
+        * c_bar_prime
+        * c_bar_prime
+        * c_bar_prime
+        * c_bar_prime
+        * c_bar_prime
+        * c_bar_prime;
+    let r_c: T = from_f64::<T>(2.0)
+        * (c_bar_prime_pow_seven / (c_bar_prime_pow_seven + twenty_five_pow_seven)).sqrt();
+    let r_t = -r_c * (from_f64::<T>(2.0) * delta_theta * pi_over_180).sin();
+
+    let one = from_f64::<T>(1.0);
+    let k_l = one;
+    let k_c = one;
+    let k_h = one;
+    let delta_l_prime = other.l - this.l;
+    let delta_c_prime = c_two_prime - c_one_prime;
+
+    ((delta_l_prime / (k_l * s_l)) * (delta_l_prime / (k_l * s_l))
+        + (delta_c_prime / (k_c * s_c)) * (delta_c_prime / (k_c * s_c))
+        + (delta_big_h_prime / (k_h * s_h)) * (delta_big_h_prime / (k_h * s_h))
+        + (r_t * delta_c_prime * delta_big_h_prime) / (k_c * s_c * k_h * s_h))
+        .sqrt()
+}

palette/src/lab.rs

@@ -7,6 +7,9 @@ use {clamp, from_f64};
 use {Alpha, LabHue, Lch, Xyz};
 use {Component, ComponentWise, FloatComponent, GetHue, Limited, Mix, Pixel, Shade};
 
+use color_difference::ColorDifference;
+use color_difference::{get_ciede_difference, LabColorDiff};
+
 /// CIE L\*a\*b\* (CIELAB) with an alpha component. See the [`Laba`
 /// implementation in `Alpha`](struct.Alpha.html#Laba).
 pub type Laba<Wp, T = f32> = Alpha<Lab<Wp, T>, T>;
@@ -303,6 +306,34 @@ where
     }
 }
 
+impl<Wp, T> ColorDifference for Lab<Wp, T>
+where
+    T: FloatComponent,
+    Wp: WhitePoint,
+{
+    type Scalar = T;
+
+    fn get_color_difference(&self, other: &Lab<Wp, T>) -> Self::Scalar {
+        // Color difference calculation requires Lab and chroma components. This
+        // function handles the conversion into those components which are then
+        // passed to `get_ciede_difference()` where calculation is completed.
+        let self_params = LabColorDiff {
+            l: self.l,
+            a: self.a,
+            b: self.b,
+            chroma: (self.a * self.a + self.b * self.b).sqrt(),
+        };
+        let other_params = LabColorDiff {
+            l: other.l,
+            a: other.a,
+            b: other.b,
+            chroma: (other.a * other.a + other.b * other.b).sqrt(),
+        };
+
+        get_ciede_difference(&self_params, &other_params)
+    }
+}
+
 impl<Wp, T> ComponentWise for Lab<Wp, T>
 where
     T: FloatComponent,

palette/src/lch.rs

@@ -1,6 +1,8 @@
 use core::marker::PhantomData;
 use core::ops::{Add, AddAssign, Sub, SubAssign};
 
+use color_difference::ColorDifference;
+use color_difference::{get_ciede_difference, LabColorDiff};
 use encoding::pixel::RawPixel;
 use white_point::{WhitePoint, D65};
 use {clamp, from_f64};
@@ -303,6 +305,53 @@ where
     }
 }
 
+/// CIEDE2000 distance metric for color difference.
+impl<Wp, T> ColorDifference for Lch<Wp, T>
+where
+    T: FloatComponent,
+    Wp: WhitePoint,
+{
+    type Scalar = T;
+
+    fn get_color_difference(&self, other: &Lch<Wp, T>) -> Self::Scalar {
+        // Prepare a* and b* from Lch components to calculate color difference
+        let self_a = clamp(
+            self.chroma.max(T::zero()) * self.hue.to_radians().cos(),
+            from_f64(-128.0),
+            from_f64(127.0),
+        );
+        let self_b = clamp(
+            self.chroma.max(T::zero()) * self.hue.to_radians().sin(),
+            from_f64(-128.0),
+            from_f64(127.0),
+        );
+        let other_a = clamp(
+            other.chroma.max(T::zero()) * other.hue.to_radians().cos(),
+            from_f64(-128.0),
+            from_f64(127.0),
+        );
+        let other_b = clamp(
+            other.chroma.max(T::zero()) * other.hue.to_radians().sin(),
+            from_f64(-128.0),
+            from_f64(127.0),
+        );
+        let self_params = LabColorDiff {
+            l: self.l,
+            a: self_a,
+            b: self_b,
+            chroma: self.chroma,
+        };
+        let other_params = LabColorDiff {
+            l: other.l,
+            a: other_a,
+            b: other_b,
+            chroma: other.chroma,
+        };
+
+        get_ciede_difference(&self_params, &other_params)
+    }
+}
+
 impl<Wp, T> Saturate for Lch<Wp, T>
 where
     T: FloatComponent,

palette/src/lib.rs

@@ -183,6 +183,7 @@ pub use rgb::{GammaSrgb, GammaSrgba, LinSrgb, LinSrgba, Srgb, Srgba};
 pub use xyz::{Xyz, Xyza};
 pub use yxy::{Yxy, Yxya};
 
+pub use color_difference::ColorDifference;
 pub use component::*;
 pub use convert::{ConvertFrom, ConvertInto, FromColor, IntoColor, OutOfBounds};
 pub use encoding::pixel::Pixel;
@@ -364,6 +365,7 @@ mod yxy;
 mod hues;
 
 pub mod chromatic_adaptation;
+mod color_difference;
 mod component;
 mod convert;
 pub mod encoding;

palette/tests/convert/data_ciede_2000.csv

@@ -0,0 +1,35 @@
+lab1_l,lab1_a,lab1_b,lab2_l,lab2_a,lab2_b,delta_e
+50.0000,2.6772,-79.7751,50.0000,0.0000,-82.7485,2.0425
+50.0000,3.1571,-77.2803,50.0000,0.0000,-82.7485,2.8615
+50.0000,2.8361,-74.0200,50.0000,0.0000,-82.7485,3.4412
+50.0000,-1.3802,-84.2814,50.0000,0.0000,-82.7485,1.0000
+50.0000,-1.1848,-84.8006,50.0000,0.0000,-82.7485,1.0000
+50.0000,-0.9009,-85.5211,50.0000,0.0000,-82.7485,1.0000
+50.0000,0.0000,0.0000,50.0000,-1.0000,2.0000,2.3669
+50.0000,-1.0000,2.0000,50.0000,0.0000,0.0000,2.3669
+50.0000,2.4900,-0.0010,50.0000,-2.4900,0.0009,7.1792
+50.0000,2.4900,-0.0010,50.0000,-2.4900,0.0010,7.1792
+50.0000,2.4900,-0.0010,50.0000,-2.4900,0.0011,7.2195
+50.0000,2.4900,-0.0010,50.0000,-2.4900,0.0012,7.2195
+50.0000,-0.0010,2.4900,50.0000,0.0009,-2.4900,4.8045
+50.0000,-0.0010,2.4900,50.0000,0.0010,-2.4900,4.8045
+50.0000,-0.0010,2.4900,50.0000,0.0011,-2.4900,4.7461
+50.0000,2.5000,0.0000,50.0000,0.0000,-2.5000,4.3065
+50.0000,2.5000,0.0000,73.0000,25.0000,-18.0000,27.1492
+50.0000,2.5000,0.0000,61.0000,-5.0000,29.0000,22.8977
+50.0000,2.5000,0.0000,56.0000,-27.0000,-3.0000,31.9030
+50.0000,2.5000,0.0000,58.0000,24.0000,15.0000,19.4535
+50.0000,2.5000,0.0000,50.0000,3.1736,0.5854,1.0000
+50.0000,2.5000,0.0000,50.0000,3.2972,0.0000,1.0000
+50.0000,2.5000,0.0000,50.0000,1.8634,0.5757,1.0000
+50.0000,2.5000,0.0000,50.0000,3.2592,0.3350,1.0000
+60.2574,-34.0099,36.2677,60.4626,-34.1751,39.4387,1.2644
+63.0109,-31.0961,-5.8663,62.8187,-29.7946,-4.0864,1.2630
+61.2901,3.7196,-5.3901,61.4292,2.2480,-4.9620,1.8731
+35.0831,-44.1164,3.7933,35.0232,-40.0716,1.5901,1.8645
+22.7233,20.0904,-46.6940,23.0331,14.9730,-42.5619,2.0373
+36.4612,47.8580,18.3852,36.2715,50.5065,21.2231,1.4146
+90.8027,-2.0831,1.4410,91.1528,-1.6435,0.0447,1.4441
+90.9257,-0.5406,-0.9208,88.6381,-0.8985,-0.7239,1.5381
+6.7747,-0.2908,-2.4247,5.8714,-0.0985,-2.2286,0.6377
+2.0776,0.0795,-1.1350,0.9033,-0.0636,-0.5514,0.9082

palette/tests/convert/data_ciede_2000.rs

@@ -0,0 +1,80 @@
+/*
+Data from http://www2.ece.rochester.edu/~gsharma/ciede2000/
+
+Tests Lab color differences with expected delta E*
+
+Note: Test uses `f64` because `f32` failed Travis CI builds on Linux for Lch on
+        case 13 or 14 which is noted in the paper as testing accuracy of hue
+        angle and atan calcuation (calculated: 4.7460666, expected: 4.8045).
+        MacOS and Windows passed the tests so be wary when using f32 on Linux.
+*/
+
+extern crate approx;
+
+use csv;
+use palette::white_point::D65;
+use palette::ColorDifference;
+use palette::{Lab, Lch};
+
+#[derive(Deserialize, PartialEq)]
+struct Cie2000Raw {
+    lab1_l: f64,
+    lab1_a: f64,
+    lab1_b: f64,
+    lab2_l: f64,
+    lab2_a: f64,
+    lab2_b: f64,
+    delta_e: f64,
+}
+
+#[derive(Copy, Clone, PartialEq, Debug)]
+struct Cie2000 {
+    c1: Lab<D65, f64>,
+    c2: Lab<D65, f64>,
+    delta_e: f64,
+}
+
+impl From<Cie2000Raw> for Cie2000 {
+    fn from(src: Cie2000Raw) -> Cie2000 {
+        Cie2000 {
+            c1: Lab::new(src.lab1_l, src.lab1_a, src.lab1_b),
+            c2: Lab::new(src.lab2_l, src.lab2_a, src.lab2_b),
+            delta_e: src.delta_e,
+        }
+    }
+}
+
+fn load_data() -> Vec<Cie2000> {
+    let file_name = "tests/convert/data_ciede_2000.csv";
+    let mut rdr = csv::Reader::from_path(file_name)
+        .expect("csv file could not be loaded in tests for cie 2000 data");
+    let mut color_data: Vec<Cie2000> = Vec::new();
+    for record in rdr.deserialize() {
+        let r: Cie2000Raw =
+            record.expect("color data could not be decoded in tests for cie 2000 data");
+        color_data.push(r.into())
+    }
+    color_data
+}
+
+fn check_equal_lab(result: f64, expected: f64) {
+    assert_relative_eq!(result, expected, epsilon = 0.0001);
+}
+
+fn check_equal_lch(result: f64, expected: f64) {
+    assert_relative_eq!(result, expected, epsilon = 0.0001);
+}
+
+pub fn run_tests() {
+    let data = load_data();
+
+    for expected in data.iter() {
+        let result_lab = expected.c1.get_color_difference(&expected.c2);
+        check_equal_lab(result_lab, expected.delta_e);
+
+        let lch1: Lch<_, f64> = Lch::from(expected.c1);
+        let lch2: Lch<_, f64> = Lch::from(expected.c2);
+        let result_lch = lch1.get_color_difference(&lch2);
+        check_equal_lch(result_lch, expected.delta_e);
+    }
+}

palette/tests/convert/mod.rs

@@ -1,4 +1,5 @@
 mod data_cie_15_2004;
+mod data_ciede_2000;
 mod data_color_mine;
 mod lab_lch;
 
@@ -7,6 +8,11 @@ pub fn xyz_yxy_conversion() {
     data_cie_15_2004::run_tests();
 }
 
+#[test]
+pub fn color_difference_ciede() {
+    data_ciede_2000::run_tests();
+}
+
 #[test]
 pub fn color_mine_from_xyz() {
     data_color_mine::run_from_xyz_tests();