From f1fd2a1edd582a01c28bf96710fbcaa554fc18a3 Mon Sep 17 00:00:00 2001 From: Lucien Cartier-Tilet Date: Thu, 20 Aug 2026 23:21:45 +0200 Subject: [PATCH] feat: add color module with sRGB to Lab conversion (scalar + SIMD) - src/color.rs: Rgb/Lab types, planar SoA storage, wide-based SIMD RGB to Lab conversion with scalar remainder tail - src/palette.rs: NORD_PALETTE constant - Cargo.toml: add image/ndarray/rayon/wide deps, tuned release profile --- Cargo.toml | 8 + rust-toolchain.toml | 2 +- src/color.rs | 491 ++++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 3 + src/palette.rs | 20 ++ 5 files changed, 523 insertions(+), 1 deletion(-) create mode 100644 src/color.rs create mode 100644 src/palette.rs diff --git a/Cargo.toml b/Cargo.toml index 81678ae..4f871f9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,3 +15,11 @@ path = "src/main.rs" name = "nordify" [dependencies] +image = "0.25.10" +ndarray = { version = "0.17.2", features = ["rayon"] } +rayon = "1.12.0" +wide = "1.6.1" + +[profile.release] +lto = true +codegen-units = 1 diff --git a/rust-toolchain.toml b/rust-toolchain.toml index 9ce1c37..ea9b96b 100644 --- a/rust-toolchain.toml +++ b/rust-toolchain.toml @@ -1,4 +1,4 @@ [toolchain] -channel = "stable" +channel = "nightly" components = ["clippy", "rust-docs", "rust-src", "rust-analyzer", "rustfmt"] profile = "default" \ No newline at end of file diff --git a/src/color.rs b/src/color.rs new file mode 100644 index 0000000..95e5d3b --- /dev/null +++ b/src/color.rs @@ -0,0 +1,491 @@ +use std::sync::LazyLock; +use wide::f64x4; + +pub trait Color {} + +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub struct Rgb { + pub r: u8, + pub g: u8, + pub b: u8, +} + +impl Color for Rgb {} + +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Lab { + pub l: f64, + pub a: f64, + pub b: f64, +} + +impl Color for Lab {} + +/// Constants for the sRGB transfer function. `srgb_to_linear` uses these +/// to convert a normalized (`0.0..=1.0`) sRGB channel value to linear-light +/// intensity: +/// +/// ```txt +/// srgb_channel > SRGB_LINEAR_THRESHOLD +/// ? ((srgb_channel + SRGB_ALPHA) / (1.0 + SRGB_ALPHA)) ^ SRGB_GAMMA +/// : srgb_channel / SRGB_LINEAR_DIVISOR +/// ``` +/// +/// The sRGB standard (IEC 61966-2-1) fixes all four values. They are not +/// tunable. A pure power curve has an infinite slope at zero. This causes +/// numerical instability for very dark values. The sRGB standard uses a +/// plain linear segment near black instead. `SRGB_LINEAR_THRESHOLD` and +/// `SRGB_LINEAR_DIVISOR` make the linear segment match the power-law +/// segment at the threshold. Both value and slope match. This leaves no +/// visible seam between the two pieces. +const SRGB_LINEAR_THRESHOLD: f64 = 0.04045; +const SRGB_ALPHA: f64 = 0.055; +const SRGB_GAMMA: f64 = 2.4; +const SRGB_LINEAR_DIVISOR: f64 = 12.92; + +/// The linear sRGB to CIE 1931 XYZ conversion matrix (D65 white point). +/// +/// Each row multiplies a `(linear_r, linear_g, linear_b)` triple to produce +/// one XYZ tristimulus value. Row 0 produces X. Row 1 produces Y (luminance). +/// Row 2 produces Z. See `linear_rgb_to_lab` and `linear_rgb_to_lab_simd`, +/// which compute `XYZ_MATRIX[i][0] * linear_r + XYZ_MATRIX[i][1] * linear_g + +/// XYZ_MATRIX[i][2] * linear_b` for each row `i`. This is the middle step +/// of RGB to Lab. It must use linear-light RGB (after `SRGB_TO_LINEAR`). +/// It must not use gamma-encoded values. The transform is only physically +/// meaningful between two linear color spaces. +/// +/// The values come from the sRGB standard (IEC 61966-2-1). The standard +/// defines three primary chromaticities and the D65 white point. The matrix +/// solves for mixing the primaries at full intensity to reproduce D65 white. +/// This is a standard, publicly fixed matrix. It matches the `nordify` +/// Python reference implementation. +const XYZ_MATRIX: [[f64; 3]; 3] = [ + [0.4124564, 0.3575761, 0.1804375], + [0.2126729, 0.7151522, 0.0721750], + [0.0193339, 0.1191920, 0.9503041], +]; + +/// The CIE standard illuminant D65 reference white point in CIE XYZ +/// coordinates. `XYZ_MATRIX` produces XYZ values scaled so the sRGB +/// primaries mix to D65 white. Dividing by `WHITE_POINT` (in +/// `linear_rgb_to_lab` and `linear_rgb_to_lab_simd`, right before +/// `xyz_compress`) re-normalizes so white maps to `(1.0, 1.0, 1.0)`. +/// This makes the Lab convention "white has L* = 100" work correctly. +/// `Y` is exactly `1.0` by definition. Y represents luminance, +/// normalized so white has luminance 1. `X` and `Z` are not `1.0` +/// because D65 white does not have equal energy in those directions. +const WHITE_POINT: [f64; 3] = [0.95047, 1.00000, 1.08883]; + +/// Constants for the CIE Lab forward compression function `f(t)`. +/// `xyz_compress` and `xyz_compress_simd` use these to convert a +/// white-point-normalized XYZ component into the perceptually-uniform +/// value that the L*/a*/b* formulas use: +/// +/// ```txt +/// t > LAB_CUBE_THRESHOLD ? cbrt(t) : LAB_LINEAR_SLOPE * t + LAB_LINEAR_INTERCEPT +/// ``` +/// +/// This piecewise split avoids the numerical instability of a cube root +/// infinite slope near zero. CIE fixes all three values so the linear +/// segment meets the cube-root curve smoothly at `LAB_CUBE_THRESHOLD`. +/// `LAB_LINEAR_INTERCEPT` (`16/116`) has a second role. The final +/// `L = 116 * f(Y) - 16` formula subtracts this value back out. A +/// pure-black input (`f(0) = 16/116`) maps to `L = 0`. +const LAB_CUBE_THRESHOLD: f64 = 0.008856; +const LAB_LINEAR_SLOPE: f64 = 7.787; +const LAB_LINEAR_INTERCEPT: f64 = 16.0 / 116.0; + +/// Converts a normalized (`0.0..=1.0`) gamma-corrected sRGB channel +/// value to its linear-light equivalent. This reverses the sRGB transfer +/// function. +fn srgb_to_linear(srgb_channel: f64) -> f64 { + if srgb_channel > SRGB_LINEAR_THRESHOLD { + ((srgb_channel + SRGB_ALPHA) / (1.0 + SRGB_ALPHA)).powf(SRGB_GAMMA) + } else { + srgb_channel / SRGB_LINEAR_DIVISOR + } +} + +/// Precomputed `srgb_to_linear` result for every possible `u8` channel +/// value (`0..=255`). This table builds once from the scalar formula so +/// it cannot drift out of sync. An sRGB channel takes 256 distinct values. +/// This table turns "linearize a channel" into an O(1) array read +/// (`SRGB_TO_LINEAR[value as usize]`) instead of a `powf` call. The SIMD +/// chunk path (`linearize_chunk_simd`) uses it with 4 lookups packed into +/// an `f64x4`. Scalar code (`From for Lab`, `convert_from` remainder +/// tail) also uses it. This table avoids `wide`'s `powf_simd`, whose +/// precision is unspecified and can vary by platform. This table gives +/// deterministic results everywhere. This matters for tolerance-based tests. +static SRGB_TO_LINEAR: LazyLock<[f64; 256]> = LazyLock::new(|| { + let mut table = [0.0f64; 256]; + for (channel_value, linear) in table.iter_mut().enumerate() { + *linear = srgb_to_linear(channel_value as f64 / 255.0); + } + table +}); + +/// Applies the CIE Lab forward compression function `f(t)` to a single +/// white-point-normalized XYZ tristimulus component. Use this to derive +/// L*/a*/b* from X/Y/Z. +fn xyz_compress(xyz_component: f64) -> f64 { + if xyz_component > LAB_CUBE_THRESHOLD { + xyz_component.cbrt() + } else { + LAB_LINEAR_SLOPE * xyz_component + LAB_LINEAR_INTERCEPT + } +} + +fn xyz_compress_simd(xyz_component: f64x4) -> f64x4 { + let cube_root = xyz_component.cbrt(); + let linear = + f64x4::splat(LAB_LINEAR_SLOPE) * xyz_component + f64x4::splat(LAB_LINEAR_INTERCEPT); + xyz_component + .simd_gt(f64x4::splat(LAB_CUBE_THRESHOLD)) + .select(cube_root, linear) +} + +/// Computes Lab from linear-light RGB (after sRGB decoding). The scalar +/// `From for Lab` and `convert_from` scalar remainder tail share this +/// function. +fn linear_rgb_to_lab(linear_r: f64, linear_g: f64, linear_b: f64) -> Lab { + let x = xyz_compress( + (XYZ_MATRIX[0][0] * linear_r + XYZ_MATRIX[0][1] * linear_g + XYZ_MATRIX[0][2] * linear_b) + / WHITE_POINT[0], + ); + let y = xyz_compress( + (XYZ_MATRIX[1][0] * linear_r + XYZ_MATRIX[1][1] * linear_g + XYZ_MATRIX[1][2] * linear_b) + / WHITE_POINT[1], + ); + let z = xyz_compress( + (XYZ_MATRIX[2][0] * linear_r + XYZ_MATRIX[2][1] * linear_g + XYZ_MATRIX[2][2] * linear_b) + / WHITE_POINT[2], + ); + + Lab { + l: 116.0 * y - 16.0, + a: 500.0 * (x - y), + b: 200.0 * (y - z), + } +} + +/// Converts 4 linear-light RGB triples at once. This is the vectorized +/// counterpart to `linear_rgb_to_lab`. It returns `(l, a, b)` lane vectors. +fn linear_rgb_to_lab_simd( + linear_r: f64x4, + linear_g: f64x4, + linear_b: f64x4, +) -> (f64x4, f64x4, f64x4) { + let x = xyz_compress_simd( + (f64x4::splat(XYZ_MATRIX[0][0]) * linear_r + + f64x4::splat(XYZ_MATRIX[0][1]) * linear_g + + f64x4::splat(XYZ_MATRIX[0][2]) * linear_b) + / f64x4::splat(WHITE_POINT[0]), + ); + let y = xyz_compress_simd( + (f64x4::splat(XYZ_MATRIX[1][0]) * linear_r + + f64x4::splat(XYZ_MATRIX[1][1]) * linear_g + + f64x4::splat(XYZ_MATRIX[1][2]) * linear_b) + / f64x4::splat(WHITE_POINT[1]), + ); + let z = xyz_compress_simd( + (f64x4::splat(XYZ_MATRIX[2][0]) * linear_r + + f64x4::splat(XYZ_MATRIX[2][1]) * linear_g + + f64x4::splat(XYZ_MATRIX[2][2]) * linear_b) + / f64x4::splat(WHITE_POINT[2]), + ); + + let l = f64x4::splat(116.0) * y - f64x4::splat(16.0); + let a = f64x4::splat(500.0) * (x - y); + let b = f64x4::splat(200.0) * (y - z); + + (l, a, b) +} + +/// Linearizes 4 raw (`0..=255`) sRGB channel values at once. This uses 4 +/// independent `SRGB_TO_LINEAR` lookups packed into one vector. Table +/// lookups avoid `powf_simd` documented precision non-determinism. +fn linearize_chunk_simd(raw_channel_chunk: &[f64]) -> f64x4 { + f64x4::new([ + SRGB_TO_LINEAR[raw_channel_chunk[0] as usize], + SRGB_TO_LINEAR[raw_channel_chunk[1] as usize], + SRGB_TO_LINEAR[raw_channel_chunk[2] as usize], + SRGB_TO_LINEAR[raw_channel_chunk[3] as usize], + ]) +} + +impl From for Lab { + /// Converts an 8-bit sRGB color to CIE Lab. The conversion uses + /// linear-light RGB and XYZ with a D65 white point. This gives + /// perceptually accurate color distance. + fn from(value: Rgb) -> Self { + let linear_r = SRGB_TO_LINEAR[value.r as usize]; + let linear_g = SRGB_TO_LINEAR[value.g as usize]; + let linear_b = SRGB_TO_LINEAR[value.b as usize]; + linear_rgb_to_lab(linear_r, linear_g, linear_b) + } +} + +pub trait PlanarColor { + fn len(&self) -> usize; + fn is_empty(&self) -> bool; + fn with_capacity(n: usize) -> Self; + fn clear(&mut self); + fn push(&mut self, color: T); +} + +/// Planar (struct-of-arrays) RGB storage. One contiguous `Vec` holds +/// each channel. Channel values are raw (`0..=255`, matching the `u8` +/// inputs). They are not normalized to `0..1`. This keeps `SRGB_TO_LINEAR` +/// indexing (`value as usize`) exact. Reconstructing an integer index from +/// a normalized float would lose precision. Distance and argmin comparisons +/// work the same either way. Scaling every coordinate by the same constant +/// preserves relative distances. +pub struct RgbPlanar { + r: Vec, + g: Vec, + b: Vec, +} + +impl RgbPlanar { + pub fn new(r: Vec, g: Vec, b: Vec) -> Self { + assert_eq!(r.len(), g.len()); + assert_eq!(r.len(), b.len()); + Self { r, g, b } + } + + pub fn r(&self) -> &[f64] { + &self.r + } + pub fn g(&self) -> &[f64] { + &self.g + } + pub fn b(&self) -> &[f64] { + &self.b + } +} + +impl PlanarColor for RgbPlanar { + fn len(&self) -> usize { + self.r.len() + } + + fn is_empty(&self) -> bool { + self.r.is_empty() + } + + fn with_capacity(n: usize) -> Self { + Self { + r: Vec::with_capacity(n), + g: Vec::with_capacity(n), + b: Vec::with_capacity(n), + } + } + + fn clear(&mut self) { + self.r.clear(); + self.g.clear(); + self.b.clear(); + } + + fn push(&mut self, color: Rgb) { + self.r.push(f64::from(color.r)); + self.g.push(f64::from(color.g)); + self.b.push(f64::from(color.b)); + } +} + +pub struct LabPlanar { + l: Vec, + a: Vec, + b: Vec, +} + +impl LabPlanar { + pub fn new(l: Vec, a: Vec, b: Vec) -> Self { + assert_eq!(l.len(), a.len()); + assert_eq!(l.len(), b.len()); + Self { l, a, b } + } + + pub fn l(&self) -> &[f64] { + &self.l + } + pub fn a(&self) -> &[f64] { + &self.a + } + pub fn b(&self) -> &[f64] { + &self.b + } + + /// Repopulates `self` from `rgb` per-color RGB values. This converts + /// to Lab via linear-light RGB and XYZ. It reuses the existing + /// allocation in `self`. It clears and refills in place instead of + /// allocating a new `LabPlanar`. This processes 4 colors at a time + /// via `wide::f64x4`. A scalar tail handles `rgb.len() % 4` leftover + /// colors. + pub fn convert_from(&mut self, rgb: &RgbPlanar) { + self.clear(); + let len = rgb.len(); + let chunks = len / 4; + + for chunk_index in 0..chunks { + let i = chunk_index * 4; + let linear_r = linearize_chunk_simd(&rgb.r()[i..i + 4]); + let linear_g = linearize_chunk_simd(&rgb.g()[i..i + 4]); + let linear_b = linearize_chunk_simd(&rgb.b()[i..i + 4]); + + let (l, a, b) = linear_rgb_to_lab_simd(linear_r, linear_g, linear_b); + let (l, a, b) = (l.to_array(), a.to_array(), b.to_array()); + + for lane in 0..4 { + self.push(Lab { + l: l[lane], + a: a[lane], + b: b[lane], + }); + } + } + + for i in (chunks * 4)..len { + let linear_r = SRGB_TO_LINEAR[rgb.r()[i] as usize]; + let linear_g = SRGB_TO_LINEAR[rgb.g()[i] as usize]; + let linear_b = SRGB_TO_LINEAR[rgb.b()[i] as usize]; + self.push(linear_rgb_to_lab(linear_r, linear_g, linear_b)); + } + } +} + +impl PlanarColor for LabPlanar { + fn len(&self) -> usize { + self.l.len() + } + + fn is_empty(&self) -> bool { + self.l.is_empty() + } + + fn with_capacity(n: usize) -> Self { + Self { + l: Vec::with_capacity(n), + a: Vec::with_capacity(n), + b: Vec::with_capacity(n), + } + } + + fn clear(&mut self) { + self.l.clear(); + self.a.clear(); + self.b.clear(); + } + + fn push(&mut self, color: Lab) { + self.l.push(color.l); + self.a.push(color.a); + self.b.push(color.b); + } +} + +impl From<&RgbPlanar> for LabPlanar { + fn from(value: &RgbPlanar) -> Self { + let mut out = Self::with_capacity(value.len()); + out.convert_from(value); + out + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::palette::NORD_PALETTE; + + const TOLERANCE: f64 = 1e-6; + + fn assert_lab_close(actual: Lab, expected: Lab) { + assert!( + (actual.l - expected.l).abs() < TOLERANCE, + "l: {} vs {}", + actual.l, + expected.l + ); + assert!( + (actual.a - expected.a).abs() < TOLERANCE, + "a: {} vs {}", + actual.a, + expected.a + ); + assert!( + (actual.b - expected.b).abs() < TOLERANCE, + "b: {} vs {}", + actual.b, + expected.b + ); + } + + fn edge_case_colors() -> Vec { + let mut colors: Vec = NORD_PALETTE.to_vec(); + colors.push(Rgb { r: 0, g: 0, b: 0 }); + colors.push(Rgb { + r: 255, + g: 255, + b: 255, + }); + colors.push(Rgb { r: 255, g: 0, b: 0 }); + colors + } + + fn to_rgb_planar(colors: &[Rgb]) -> RgbPlanar { + let mut planar = RgbPlanar::with_capacity(colors.len()); + for &color in colors { + planar.push(color); + } + planar + } + + fn assert_lab_planar_matches_scalar(colors: &[Rgb], lab_planar: &LabPlanar) { + assert_eq!(lab_planar.len(), colors.len()); + for (index, &color) in colors.iter().enumerate() { + let expected = Lab::from(color); + let actual = Lab { + l: lab_planar.l()[index], + a: lab_planar.a()[index], + b: lab_planar.b()[index], + }; + assert_lab_close(actual, expected); + } + } + + #[test] + fn lab_planar_from_rgb_planar_matches_scalar() { + let colors = edge_case_colors(); + let rgb_planar = to_rgb_planar(&colors); + let lab_planar = LabPlanar::from(&rgb_planar); + + assert_lab_planar_matches_scalar(&colors, &lab_planar); + } + + #[test] + fn lab_planar_convert_from_handles_non_multiple_of_lane_width() { + for len in [5, 7] { + let colors: Vec = edge_case_colors().into_iter().take(len).collect(); + let rgb_planar = to_rgb_planar(&colors); + let lab_planar = LabPlanar::from(&rgb_planar); + + assert_lab_planar_matches_scalar(&colors, &lab_planar); + } + } + + #[test] + fn lab_planar_convert_from_reuses_allocation() { + let long_colors = edge_case_colors(); + let short_colors: Vec = long_colors.iter().take(3).copied().collect(); + + let long_rgb = to_rgb_planar(&long_colors); + let short_rgb = to_rgb_planar(&short_colors); + + let mut lab_planar = LabPlanar::with_capacity(long_colors.len()); + lab_planar.convert_from(&long_rgb); + lab_planar.convert_from(&short_rgb); + + assert_lab_planar_matches_scalar(&short_colors, &lab_planar); + } +} diff --git a/src/lib.rs b/src/lib.rs index 6eec06a..82d5d02 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,6 @@ +mod color; +mod palette; + pub fn add(a: u32, b: u32) -> u32 { a + b } diff --git a/src/palette.rs b/src/palette.rs new file mode 100644 index 0000000..15d77b8 --- /dev/null +++ b/src/palette.rs @@ -0,0 +1,20 @@ +use crate::color::Rgb; + +pub const NORD_PALETTE: [Rgb; 16] = [ + Rgb {r: 46, g: 52, b: 64,}, // nord0 + Rgb {r: 59, g: 66, b: 82,}, // nord1 + Rgb {r: 67, g: 76, b: 94,}, // nord2 + Rgb {r: 76, g: 86, b: 106,}, // nord3 + Rgb {r: 216, g: 222, b: 233,}, // nord4 + Rgb {r: 229, g: 233, b: 240,}, // nord5 + Rgb {r: 236, g: 239, b: 244,}, // nord6 + Rgb {r: 143, g: 188, b: 187,}, // nord7 + Rgb {r: 136, g: 192, b: 208,}, // nord8 + Rgb {r: 129, g: 161, b: 193,}, // nord9 + Rgb {r: 94, g: 129, b: 172,}, // nord10 + Rgb {r: 191, g: 97, b: 106,}, // nord11 + Rgb {r: 208, g: 135, b: 112,}, // nord12 + Rgb {r: 235, g: 203, b: 139,}, // nord13 + Rgb {r: 163, g: 190, b: 140,}, // nord14 + Rgb {r: 180, g: 142, b: 173,}, // nord15 +];