Testing out color distance formulas in Rust

2024-10-22 09:16:57 -07:00
parent 5d5ea914c7
commit ba1337147b
5 changed files with 193 additions and 14 deletions
@@ -241,8 +241,11 @@ version = "0.1.0"
 dependencies = [
 "anyhow",
 "clap",
+ "empfindung",
 "image",
+ "lab",
 "ordered-float",
+ "rand",
 ]

 [[package]]
@@ -303,6 +306,16 @@ version = "1.13.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0"

+[[package]]
+name = "empfindung"
+version = "0.2.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8af095ec32b4a626b56a2e188515e4e5480d82b02e5d8a830fdb61934226b0df"
+dependencies = [
+ "lab",
+ "rgb",
+]
+
 [[package]]
 name = "equivalent"
 version = "1.0.1"
@@ -486,6 +499,12 @@ version = "0.3.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "f5d4a7da358eff58addd2877a45865158f0d78c911d43a5784ceb7bbf52833b0"

+[[package]]
+name = "lab"
+version = "0.11.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bf36173d4167ed999940f804952e6b08197cae5ad5d572eb4db150ce8ad5d58f"
+
 [[package]]
 name = "lebe"
 version = "0.5.2"
@@ -3,11 +3,21 @@ name = "color-cell"
 version = "0.1.0"
 edition = "2021"

+[profile.dev.package.zune-jpeg]
+opt-level = 3
+
+[profile.dev.package.lab]
+opt-level = 3
+
+[profile.dev.package.empfindung]
+opt-level = 3
+
 [dependencies]
 anyhow = "1.0.86"
 clap = { version = "4.5.15", features = ["derive"] }
+empfindung = "0.2.6"
 image = "0.25.2"
+lab = "0.11.0"
 ordered-float = "4.2.2"
+rand = "0.8.5"

-[profile.dev.package.zune-jpeg]
-opt-level = 3
@@ -0,0 +1,87 @@
+use empfindung::cie00;
+use lab::Lab;
+
+use crate::indexed::redmean_sq;
+
+
+fn ciede2000(a: [u8;3], b: [u8;3]) -> f32 {
+    cie00::diff(Lab::from_rgb(&a), Lab::from_rgb(&b)) // / 100.0
+}
+
+fn ciede2000_sq(a: [u8;3], b: [u8;3]) -> f32 {
+    let diff = ciede2000(a, b);
+    diff * diff
+}
+
+fn redmean(a: [u8;3], b: [u8;3]) -> f32 {
+    redmean_sq(a, b).sqrt()
+}
+
+#[cfg(test)]
+mod test {
+    use core::f32;
+    use std::f32::INFINITY;
+
+    use rand::prelude::*;
+
+    use super::ciede2000;
+
+    fn cost(params: [f32; 3]) -> f32 {
+        let decode = |x: u32| -> [u8; 3] {
+            let r = (x >> 8) as u8;
+            let g = ((x >> 4) & 15) as u8;
+            let b = (x & 15) as u8;
+            [r | r<<4, g | g<<4, b | b<<4]
+        };
+        let mut result = 0.0;
+        // TODO: precompute
+        for color1 in 0..4096 {
+            if color1 % 73 != 0 {
+                continue;
+            }
+            let rgb1 = decode(color1);
+            for color2 in (color1 + 1)..4096 {
+                if color2 % 73 != 0 {
+                    continue;
+                }
+                let rgb2 = decode(color2);
+                let true_diff = ciede2000(rgb1, rgb2);
+                let mut approx_diff = 0.0;
+                for i in 0..3 {
+                    let delta = rgb1[i] as f32 - rgb2[i] as f32;
+                    approx_diff += params[i] * delta * delta;
+                }
+                approx_diff = approx_diff.sqrt();
+                let err = true_diff - approx_diff;
+                result += err * err;
+            }
+        }
+        result
+    }
+
+
+    #[test]
+    fn tinker() {
+        let mut params = [1.0; 3];
+        let mut best_params = [1.0; 3];
+        let mut best_cost = f32::INFINITY;
+        let mut scale = 1.0;
+        for i in 0..5000 {
+            if i > 4990 {
+                dbg!((i, best_params, best_cost));
+            }
+            for c in 0..3 {
+                params[c] += scale * thread_rng().gen_range(-0.01..0.01f32);
+            }
+            let new_cost = cost(params);
+            if new_cost < best_cost {
+                best_cost = new_cost;
+                best_params = params;
+            } else {
+                params = best_params;
+                scale *= 0.99999;
+            }
+        }
+        dbg!((best_params, best_cost));
+    }
+}
@@ -1,14 +1,20 @@
-use image::RgbImage;
-use ordered_float::NotNan;
+use std::collections::BinaryHeap;

-// pub fn euclid_sq(x: [u8; 3], y: [u8; 3]) -> f32 {
-//     let mut result = 0.0;
-//     for i in 0..3 {
-//         let delta = x[i] as f32 - y[i] as f32;
-//         result += delta * delta;
-//     }
-//     result
-// }
+use empfindung::cie00;
+use image::RgbImage;
+use lab::Lab;
+use ordered_float::NotNan;
+use rand::Rng;
+
+pub fn euclid_sq(x: [u8; 3], y: [u8; 3]) -> f32 {
+    let mut result = 0.0;
+    let weights = [1.641, 5.764, 2.088];
+    for (i, weight) in weights.iter().enumerate() {
+        let delta = x[i] as f32 - y[i] as f32;
+        result += weight * delta * delta;
+    }
+    result
+}

 pub fn redmean_sq(x: [u8; 3], y: [u8; 3]) -> f32 {
    let mut ds: [f32; 3] = [0.0; 3];
@@ -22,6 +28,27 @@ pub fn redmean_sq(x: [u8; 3], y: [u8; 3]) -> f32 {
    result / 584971.0
 }

+fn rgb_to_ycocg(x: [u8; 3]) -> [f32; 3] {
+    let [r, g, b] = x.map(|x| x as f32);
+    let co = r - b;
+    let tmp = b + co/2.0;
+    let cg = g - tmp;
+    let y = tmp + cg/2.0;
+    [y, co, cg]
+}
+
+pub fn ycocg_sq(x: [u8; 3], y: [u8; 3]) -> f32 {
+    let x = rgb_to_ycocg(x);
+    let y = rgb_to_ycocg(y);
+    let mut result = 0.0;
+    let weights = [15.03, 1.982, 7.702];
+    for (i, weight) in weights.iter().enumerate() {
+        let delta = (x[i] as f32 - y[i] as f32)/255.0;
+        result += weight * delta * delta;
+    }
+    result
+}
+
 type Rgb = [u8; 3];

 const PICO8_PALETTE: [Rgb; 32] = [
@@ -97,19 +124,54 @@ impl PicoImage {
        // Brute-force nearest-color conversion
        // TODO begin using this everywhere
        let mut result = Self::new(img.width(), img.height());
+        let mut discrepancies: Vec<(f32, (Rgb, Rgb), f32)> = vec![];
        for y in 0..img.height() {
            let offset = (y * img.width()) as usize;
            for x in 0..img.width() {
                let original = img.get_pixel(x, y).0;
                let best = (0..16)
                    .min_by_key(|&i| {
-                        let dist_sq = redmean_sq(original, PicoColor(i).to_rgb());
-                        NotNan::new(dist_sq).unwrap()
+                        let candidate = PicoColor(i).to_rgb();
+                        // let diff = redmean_sq(original, candidate);
+                        // let diff = euclid_sq(original, candidate);
+                        let diff = ycocg_sq(original, candidate);
+                        // let diff =
+                        //     cie00::diff(Lab::from_rgb(&original), Lab::from_rgb(&candidate));
+                        NotNan::new(diff).unwrap()
                    })
                    .unwrap();
+
+                // Log places where our color formula disagrees with CIEDE2000
+                let best_rgb = PicoColor(best).to_rgb();
+                let ciede2000 = cie00::diff(Lab::from_rgb(&original), Lab::from_rgb(&best_rgb));
+                let diff = 24.36 * ycocg_sq(original, best_rgb).sqrt();
+                let discrepancy = (
+                    ciede2000,
+                    (original, best_rgb),
+                    diff
+                );
+                if discrepancies.len() < 10 {
+                    discrepancies.push(discrepancy);
+                } else {
+                    let idx = rand::thread_rng().gen_range(0..discrepancies.len());
+                    if discrepancies[idx].0 < discrepancy.0 {
+                        discrepancies[idx] = discrepancy;
+                    }
+                }
+
+                // if x % 43 == 22 && y % 43 == 22 {
+                //     let best = PicoColor(best).to_rgb();
+                //     let ciede2000 = cie00::diff(Lab::from_rgb(&original), Lab::from_rgb(&best));
+                //     let redmean = redmean_sq(original, best);
+                //     dbg!((original, best));
+                //     dbg!((redmean, ciede2000 * ciede2000 / 100_00.0));
+                // }
                result.data[offset + x as usize] = best;
            }
        }
+        for d in discrepancies {
+            dbg!(d);
+        }
        result
    }

@@ -1,3 +1,4 @@
+mod color_dist;
 mod encode;
 mod indexed;