Skip to content

fix: dithering math correctness + core performance #45

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
g4bri3lDev merged 2 commits into from
Jul 5, 2026
Merged

fix: dithering math correctness + core performance #45

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6765760
fix(dither): correct ordered-dither brightness, palette scaling, and …
davelee98 Jul 5, 2026
bafc435
perf(dither): eliminate per-pixel powf in ordered dither; hoist gamut…
davelee98 Jul 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
158 changes: 137 additions & 21 deletions packages/rust/core/src/algorithms.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -307,23 +307,47 @@ pub fn direct_map(pixels: &[u8], palette: &Palette, canonical_palette: &Palette)

// ── Ordered (Bayer) dithering ─────────────────────────────────────────────────

// 4×4 Bayer matrix, normalized to [-0.5, 0.5]. Indexed as [y % 4][x % 4].
// Values = (bayer_entry / 16.0) - 0.5 for the standard ordered Bayer matrix.
// Written as exact fractions to keep precision (all are multiples of 1/16).
// 4×4 Bayer matrix, zero-mean thresholds in (-0.5, 0.5). Indexed as [y % 4][x % 4].
// Values = (bayer_entry + 0.5) / 16.0 - 0.5, the standard zero-mean normalization —
// the +0.5 centering makes ordered dithering brightness-neutral (mean threshold = 0).
// Written as exact fractions to keep precision (all are odd multiples of 1/32).
const BAYER_4X4: [[f64; 4]; 4] = [
[-0.5000, 0.0000, -0.3750, 0.1250],
[ 0.2500, -0.2500, 0.3750, -0.1250],
[-0.3125, 0.1875, -0.4375, 0.0625],
[ 0.4375, -0.0625, 0.3125, -0.1875],
[-0.46875, 0.03125, -0.34375, 0.15625],
[ 0.28125, -0.21875, 0.40625, -0.09375],
[-0.28125, 0.21875, -0.40625, 0.09375],
[ 0.46875, -0.03125, 0.34375, -0.15625],
];

/// Threshold amplitude (in sRGB-fraction space) for ordered dither on this palette.
///
/// Ordered dithering trades quantization error for a fixed threshold pattern; the
/// threshold should be scaled to the palette's quantization step. For an evenly-spaced
/// grayscale ramp of `n` levels the step is `1/(n-1)`, so a full ±0.5 threshold is
/// `2·(n-1)×` too large (e.g. 14× for 16-level gray) and swamps fine detail with noise.
///
/// Sparse palettes (mono + color) keep the full ±0.5 spread: their transitions are
/// dominated by black↔white/ink decisions where the wide threshold is the tuned,
/// tested behavior (see the issue-#27 property test).
fn ordered_spread(palette: &Palette) -> f64 {
let is_grayscale =
palette.colors.len() >= 3 && palette.colors.iter().all(|&[r, g, b]| r == g && g == b);
if is_grayscale {
1.0 / (palette.colors.len() - 1) as f64
} else {
1.0
}
}

/// Ordered (Bayer 4×4) dither. Pixels are independent — parallelized with rayon.
///
/// The Bayer threshold is added in sRGB-fraction space, not linear. Linear-space
/// thresholding produces ~3× the perceptual spread in shadows compared to highlights
/// (the sRGB gamma is convex, so a fixed linear ±0.5 step is huge near 0 and tiny near 1).
/// sRGB-space thresholding gives uniform perceptual dot density across the tonal range,
/// matching how error diffusion already accumulates error. See GitHub issue #27.
///
/// The threshold amplitude is scaled by `ordered_spread` to the palette's quantization
/// step so dense grayscale ramps are not swamped by full-range dither noise.
pub fn ordered_dither(pixels: &[u8], width: usize, palette: &Palette) -> Vec<u8> {
ordered_dither_impl(pixels, width, palette, None)
}
Expand All @@ -344,6 +368,21 @@ fn ordered_dither_impl(
canonical_palette: Option<&Palette>,
) -> Vec<u8> {
let (_palette_linear, palette_lab) = build_palette_lab(palette);
let spread = ordered_spread(palette);

// Per-threshold gamma LUT: lut[t][v] = srgb_fraction_to_linear((v/255 + threshold_t).clamp).
// Only 16 thresholds × 256 byte values exist, so this removes all three per-pixel powf
// calls. Built from the exact same inputs as the scalar path, so results are bit-identical.
let lut: Vec<[f64; 256]> = (0..16)
.map(|t| {
let threshold = BAYER_4X4[t / 4][t % 4] * spread;
let mut row = [0.0_f64; 256];
for (v, slot) in row.iter_mut().enumerate() {
*slot = srgb_fraction_to_linear((v as f64 / 255.0 + threshold).clamp(0.0, 1.0));
}
row
})
.collect();

pixels
.par_chunks(3)
Expand All @@ -357,18 +396,12 @@ fn ordered_dither_impl(

let x = i % width;
let y = i / width;

let threshold = BAYER_4X4[y % 4][x % 4];

// Add threshold in sRGB-fraction space, then convert to linear for OKLab match.
let r_srgb = (rgb[0] as f64 / 255.0 + threshold).clamp(0.0, 1.0);
let g_srgb = (rgb[1] as f64 / 255.0 + threshold).clamp(0.0, 1.0);
let b_srgb = (rgb[2] as f64 / 255.0 + threshold).clamp(0.0, 1.0);
let t = (y % 4) * 4 + (x % 4);

let lab = rgb_to_oklab(
srgb_fraction_to_linear(r_srgb),
srgb_fraction_to_linear(g_srgb),
srgb_fraction_to_linear(b_srgb),
lut[t][rgb[0] as usize],
lut[t][rgb[1] as usize],
lut[t][rgb[2] as usize],
);
match_pixel_oklab(lab, &palette_lab, WAB) as u8
})
Expand Down Expand Up @@ -464,7 +497,11 @@ mod tests {
///
/// We measure transitions (adjacent columns with differing palette indices) per band
/// in the mid-tone region and assert the max/min ratio stays modest. Empirically the
/// linear-space implementation produces ratio ≈ 13.5; sRGB-space ≈ 4.3.
/// linear-space implementation produces ratio ≈ 13.5; sRGB-space ≈ 2.2.
///
/// The mono decision midpoint sits at OKLab L=0.5 ≈ sRGB 188, so the top two bands
/// (sRGB ≳ 192) are already in the highlight roll-off where dither activity legitimately
/// tapers; we exclude them along with the pure-black first band.
#[test]
fn ordered_dither_activity_is_perceptually_uniform() {
const W: usize = 256;
Expand Down Expand Up @@ -492,9 +529,10 @@ mod tests {
}
}

// Mid-tone bands (skip the first and last — they're near pure black/white where
// the threshold is clamped and dither activity legitimately falls off).
let mid = &transitions[1..BANDS - 1];
// Mid-tone bands: skip the pure-black first band and the top two near-white bands
// (past the mono midpoint ≈ sRGB 188) where the threshold clamps and dither
// activity legitimately falls off.
let mid = &transitions[1..BANDS - 2];
let max = *mid.iter().max().unwrap();
let min = *mid.iter().min().unwrap();
assert!(min > 0, "every mid-tone band should have at least one transition: {transitions:?}");
Expand All @@ -507,6 +545,84 @@ mod tests {
);
}

#[test]
fn bayer_matrix_is_zero_mean() {
// A biased Bayer matrix systematically darkens (or brightens) every ordered-
// dithered image. The standard (v + 0.5)/16 - 0.5 normalization sums to zero.
let sum: f64 = BAYER_4X4.iter().flatten().sum();
assert!(sum.abs() < 1e-12, "Bayer thresholds must sum to zero, got {sum}");
}

#[test]
fn ordered_dither_is_brightness_neutral_at_midpoint() {
// Find the mono decision midpoint: the smallest gray value whose threshold-free
// OKLab match flips from black (0) to white (1).
let pal = ColorScheme::Mono.palette();
let (_, pal_lab) = build_palette_lab(pal);
let midpoint = (0u16..=255)
.find(|&v| {
let lin = srgb_channel_to_linear(v as u8);
let lab = rgb_to_oklab(lin, lin, lin);
match_pixel_oklab(lab, &pal_lab, WAB) == 1
})
.expect("mono palette must have a black→white transition") as u8;

// Dither a solid midpoint-gray field. With a zero-mean threshold the white
// fraction should sit near 0.5; the old -1/32 bias produced ≈ 7/16.
const N: usize = 32;
let pixels = solid_image(midpoint, midpoint, midpoint, N, N);
let out = ordered_dither(&pixels, N, pal);
let white = out.iter().filter(|&&v| v == 1).count() as f64 / out.len() as f64;
assert!(
(white - 0.5).abs() <= 1.0 / 16.0,
"ordered dither at the mono midpoint should be ~50% white, got {white}"
);
}

#[test]
fn ordered_dither_grayscale16_tracks_ramp_and_uses_many_levels() {
// On a horizontal sRGB ramp the full ±0.5 threshold swamps the 17/255 GS16 step
// (14× too large). After scaling by ordered_spread, block-averaged output should
// track the input ramp and exercise most of the 16 levels.
const W: usize = 256;
const H: usize = 4;
let pal = ColorScheme::Grayscale16.palette();

let mut pixels = Vec::with_capacity(W * H * 3);
for _ in 0..H {
for x in 0..W {
let v = x as u8;
pixels.extend_from_slice(&[v, v, v]);
}
}
let out = ordered_dither(&pixels, W, pal);

// (a) 4-wide block average of palette level tracks the input level.
let levels: Vec<f64> = pal.colors.iter().map(|&[r, _, _]| r as f64).collect();
let mut max_err = 0.0_f64;
for bx in 0..W / 4 {
let mut sum = 0.0;
for y in 0..H {
for dx in 0..4 {
let x = bx * 4 + dx;
sum += levels[out[y * W + x] as usize];
}
}
let avg = sum / (H * 4) as f64;
let input = (bx * 4 + 1) as f64; // center-ish input value
max_err = max_err.max((avg - input).abs());
}
assert!(max_err < 24.0, "block-averaged output should track the ramp, max_err={max_err}");

// (b) most of the 16 levels are used.
let mut seen = [false; 16];
for &idx in &out {
seen[idx as usize] = true;
}
let used = seen.iter().filter(|&&s| s).count();
assert!(used >= 12, "GS16 ordered dither should use ≥12 levels, used {used}");
}

#[test]
fn serpentine_first_row_matches_raster() {
// Row 0 is always scanned left-to-right regardless of serpentine flag
Expand Down
17 changes: 12 additions & 5 deletions packages/rust/core/src/color_space.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,14 +15,21 @@ pub fn srgb_fraction_to_linear(value: f64) -> f64 {
}
}

/// Linear [0.0–1.0] → sRGB [0–255]. Inverse of `srgb_channel_to_linear`.
pub fn linear_channel_to_srgb(linear: f64) -> u8 {
let srgb = if linear <= 0.0031308 {
/// Linear [0.0–1.0] → continuous sRGB fraction [0.0–1.0]. Inverse of
/// `srgb_fraction_to_linear`; the continuous counterpart of `linear_channel_to_srgb`,
/// useful when a gamma-encoded value is needed without quantizing to a byte (e.g. applying
/// a tone curve about the perceptual mid-gray).
pub fn linear_fraction_to_srgb(linear: f64) -> f64 {
if linear <= 0.0031308 {
linear * 12.92
} else {
1.055 * linear.powf(1.0 / 2.4) - 0.055
};
(srgb * 255.0).round() as u8
}
}

/// Linear [0.0–1.0] → sRGB [0–255]. Inverse of `srgb_channel_to_linear`.
pub fn linear_channel_to_srgb(linear: f64) -> u8 {
(linear_fraction_to_srgb(linear) * 255.0).round() as u8
}


Expand Down
85 changes: 43 additions & 42 deletions packages/rust/core/src/color_space_lab.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,17 +1,12 @@
//! OKLab color space and color matching for dithering.

// sRGB -> XYZ matrix (D65 illuminant, BruceLinbloom)
const M_RGB_XYZ: [[f64; 3]; 3] = [
[0.4124564, 0.3575761, 0.1804375],
[0.2126729, 0.7151522, 0.0721750],
[0.0193339, 0.1191920, 0.9503041],
];

// XYZ -> LMS (M1, Ottosson)
const M1: [[f64; 3]; 3] = [
[0.8189330101, 0.3618667424, -0.1288597137],
[0.0329845436, 0.9293118715, 0.0361456387],
[0.0482003018, 0.2643662691, 0.6338517070],
// Linear sRGB -> LMS (Ottosson's published combined matrix). This folds the
// sRGB→XYZ and XYZ→LMS steps into one 3×3, halving the per-pixel matrix work in the
// hot path and avoiding the drift of composing separately-sourced matrices.
const M_RGB_LMS: [[f64; 3]; 3] = [
[0.4122214708, 0.5363325363, 0.0514459929],
[0.2119034982, 0.6806995451, 0.1073969566],
[0.0883024619, 0.2817188376, 0.6299787005],
];

// cbrt(LMS) -> OKLab (M2, Ottosson)
Expand All @@ -36,15 +31,11 @@ pub struct OkLab {
pub b: f64,
}

/// Linear RGB → OKLab. Pipeline: RGB → XYZ → LMS → cbrt → OKLab.
/// Linear RGB → OKLab. Pipeline: RGB → LMS → cbrt → OKLab.
pub fn rgb_to_oklab(r: f64, g: f64, b: f64) -> OkLab {
let x = M_RGB_XYZ[0][0] * r + M_RGB_XYZ[0][1] * g + M_RGB_XYZ[0][2] * b;
let y = M_RGB_XYZ[1][0] * r + M_RGB_XYZ[1][1] * g + M_RGB_XYZ[1][2] * b;
let z = M_RGB_XYZ[2][0] * r + M_RGB_XYZ[2][1] * g + M_RGB_XYZ[2][2] * b;

let l = M1[0][0] * x + M1[0][1] * y + M1[0][2] * z;
let m = M1[1][0] * x + M1[1][1] * y + M1[1][2] * z;
let s = M1[2][0] * x + M1[2][1] * y + M1[2][2] * z;
let l = M_RGB_LMS[0][0] * r + M_RGB_LMS[0][1] * g + M_RGB_LMS[0][2] * b;
let m = M_RGB_LMS[1][0] * r + M_RGB_LMS[1][1] * g + M_RGB_LMS[1][2] * b;
let s = M_RGB_LMS[2][0] * r + M_RGB_LMS[2][1] * g + M_RGB_LMS[2][2] * b;

let l_ = l.cbrt();
let m_ = m.cbrt();
Expand Down Expand Up @@ -76,18 +67,11 @@ const M2_INV: [[f64; 3]; 3] = [
[1.0, -0.0894841775, -1.2914855480],
];

// M1_inv: LMS -> XYZ (Ottosson)
const M1_INV: [[f64; 3]; 3] = [
[ 1.2270138511035211, -0.5577999806518222, 0.2812561489664678],
[-0.0405801784232806, 1.1122568696168302, -0.0716766786656012],
[-0.0763812845057069, -0.4214819784180127, 1.5861632204407947],
];

// XYZ -> linear sRGB (IEC 61966-2-1 D65)
const M_XYZ_RGB: [[f64; 3]; 3] = [
[ 3.2404542, -1.5371385, -0.4985314],
[-0.9692660, 1.8760108, 0.0415560],
[ 0.0556434, -0.2040259, 1.0572252],
// LMS -> linear sRGB (Ottosson's published combined inverse of M_RGB_LMS).
const M_LMS_RGB: [[f64; 3]; 3] = [
[ 4.0767416621, -3.3077115913, 0.2309699292],
[-1.2684380046, 2.6097574011, -0.3413193965],
[-0.0041960863, -0.7034186147, 1.7076147010],
];

/// OKLab → linear RGB. Output is clamped to [0, 1]. Inverse of `rgb_to_oklab`.
Expand All @@ -100,13 +84,9 @@ pub fn oklab_to_rgb(lab: OkLab) -> [f64; 3] {
let m = m_ * m_ * m_;
let s = s_ * s_ * s_;

let x = M1_INV[0][0] * l + M1_INV[0][1] * m + M1_INV[0][2] * s;
let y = M1_INV[1][0] * l + M1_INV[1][1] * m + M1_INV[1][2] * s;
let z = M1_INV[2][0] * l + M1_INV[2][1] * m + M1_INV[2][2] * s;

let r = M_XYZ_RGB[0][0] * x + M_XYZ_RGB[0][1] * y + M_XYZ_RGB[0][2] * z;
let g = M_XYZ_RGB[1][0] * x + M_XYZ_RGB[1][1] * y + M_XYZ_RGB[1][2] * z;
let b = M_XYZ_RGB[2][0] * x + M_XYZ_RGB[2][1] * y + M_XYZ_RGB[2][2] * z;
let r = M_LMS_RGB[0][0] * l + M_LMS_RGB[0][1] * m + M_LMS_RGB[0][2] * s;
let g = M_LMS_RGB[1][0] * l + M_LMS_RGB[1][1] * m + M_LMS_RGB[1][2] * s;
let b = M_LMS_RGB[2][0] * l + M_LMS_RGB[2][1] * m + M_LMS_RGB[2][2] * s;

[r.clamp(0.0, 1.0), g.clamp(0.0, 1.0), b.clamp(0.0, 1.0)]
}
Expand Down Expand Up @@ -151,9 +131,30 @@ mod tests {
#[test]
fn white_l_is_one() {
let lab = rgb_to_oklab(1.0, 1.0, 1.0);
assert_relative_eq!(lab.l, 1.0, epsilon = 1e-4);
assert_relative_eq!(lab.a, 0.0, epsilon = 1e-4);
assert_relative_eq!(lab.b, 0.0, epsilon = 1e-4);
// With Ottosson's combined matrix, white maps to L=1, a=b=0 to full precision.
assert_relative_eq!(lab.l, 1.0, epsilon = 1e-6);
assert_relative_eq!(lab.a, 0.0, epsilon = 1e-6);
assert_relative_eq!(lab.b, 0.0, epsilon = 1e-6);
}

#[test]
fn oklab_roundtrip_is_identity() {
// rgb_to_oklab and oklab_to_rgb must be true inverses — catches any mis-transcribed
// coefficient in the combined forward/inverse matrices (e.g. a bad row breaks white).
for &rgb in &[
[1.0, 1.0, 1.0_f64],
[0.0, 0.0, 0.0],
[0.5, 0.5, 0.5],
[0.8, 0.2, 0.1],
[0.1, 0.6, 0.9],
[0.2126, 0.7152, 0.0722],
] {
let lab = rgb_to_oklab(rgb[0], rgb[1], rgb[2]);
let back = oklab_to_rgb(lab);
for c in 0..3 {
assert_relative_eq!(back[c], rgb[c], epsilon = 1e-6);
}
}
}

#[test]
Expand Down
Loading
Loading