Pixelate: CMYK Facet 13.7x faster, Lab halftone 5.0x faster (#2588)

* Reuse Facet quadrants and halftone dots without changing pixels

* Record full-image Pixelate timings and exact-output checks

---------

Co-authored-by: Brandon Thomas <echelon@gmail.com>

crates/algo/examples/bench_pixelate.rs

@@ -0,0 +1,89 @@
+//! Full filter timing, including tile reads, filtering, writeback and pruning.
+//! RAYON_NUM_THREADS=16 cargo run --release -p photocraft-algo --example bench_pixelate -- facet 6000 4000 u8 3 output.bin
+//! Optional final argument after output.bin: rgb (default), gray, cmyk or lab.
+//! Run the same example at the base commit and the PR commit, then compare output.bin with cmp.
+use photocraft_algo::{FilterParams, apply_in};
+use photocraft_color::{ColorMode, PixelFormat, SampleType};
+use photocraft_geom::Rect;
+use photocraft_raster::Surface;
+use std::io::Write;
+use std::time::Instant;
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let args: Vec<_> = std::env::args().skip(1).collect();
+    let case = args.first().map(String::as_str).unwrap_or("facet");
+    let parse = |i: usize, fallback: i32| -> Result<i32, Box<dyn std::error::Error>> {
+        Ok(match args.get(i) {
+            Some(s) => s.parse()?,
+            None => fallback,
+        })
+    };
+    let (w, h, repeats) = (parse(1, 6000)?, parse(2, 4000)?, parse(4, 3)?);
+    if !(1..=10000).contains(&w) || !(1..=10000).contains(&h) || !(1..=9).contains(&repeats) {
+        return Err("dimensions must be 1..10000, repeats 1..9".into());
+    }
+    let depth = match args.get(3).map(String::as_str).unwrap_or("u8") {
+        "u8" => SampleType::U8,
+        "u16" => SampleType::U16,
+        "f32" => SampleType::F32,
+        _ => return Err("depth must be u8, u16 or f32".into()),
+    };
+    let params = match case {
+        "facet" => FilterParams::Facet,
+        "halftone" => FilterParams::ColorHalftone { max_radius: 8.0, angles: [108.0, 162.0, 90.0, 45.0] },
+        "halftone-small" => FilterParams::ColorHalftone { max_radius: 2.0, angles: [108.0, 162.0, 90.0, 45.0] },
+        "halftone-large" => FilterParams::ColorHalftone { max_radius: 64.0, angles: [108.0, 162.0, 90.0, 45.0] },
+        _ => return Err("case must be facet, halftone, halftone-small or halftone-large".into()),
+    };
+    let mode = match args.get(6).map(String::as_str).unwrap_or("rgb") {
+        "rgb" => ColorMode::Rgb,
+        "gray" => ColorMode::Grayscale,
+        "cmyk" => ColorMode::Cmyk,
+        "lab" => ColorMode::Lab,
+        _ => return Err("mode must be rgb, gray, cmyk or lab".into()),
+    };
+    let format = PixelFormat::new(mode, depth, true);
+    let n = format.channels();
+    let bounds = Rect::new(-17, -29, w - 17, h - 29);
+    let mut source = Surface::new(format);
+    for tile in bounds.tiles() {
+        let r = tile.rect().intersect(&bounds);
+        let mut pixels = Vec::with_capacity(r.width() as usize * r.height() as usize * n);
+        for y in r.y0..r.y1 {
+            for x in r.x0..r.x1 {
+                let seed = (x as u32).wrapping_mul(0x9e37_79b9) ^ (y as u32).wrapping_mul(0x85eb_ca6b);
+                for c in 0..n {
+                    let v = seed.rotate_left(c as u32 * 7).wrapping_add(c as u32 * 71) % 65536;
+                    pixels.push(if c + 1 == n && seed.is_multiple_of(5) { 0.0 } else { v as f32 / 65535.0 });
+                }
+            }
+        }
+        source.write_region(r, &pixels);
+    }
+    let warm = Rect::new(bounds.x0, bounds.y0, bounds.x0 + w.min(128), bounds.y0 + h.min(128));
+    drop(apply_in(&source, &params, warm, bounds, None, bounds));
+    let mut times = Vec::new();
+    for run in 0..repeats {
+        let start = Instant::now();
+        let result = apply_in(&source, &params, bounds, bounds, None, bounds);
+        let seconds = start.elapsed().as_secs_f64();
+        println!("{case} {w}x{h} {depth:?} {mode:?} run={} seconds={seconds:.6}", run + 1);
+        times.push(seconds);
+        if run + 1 == repeats
+            && let Some(path) = args.get(5)
+        {
+            let mut file = std::io::BufWriter::new(std::fs::File::create(path)?);
+            for tile in bounds.tiles() {
+                for sample in result.read_region(tile.rect().intersect(&bounds)) {
+                    file.write_all(&sample.to_bits().to_le_bytes())?;
+                }
+            }
+            file.flush()?;
+        }
+    }
+    times.sort_by(f64::total_cmp);
+    if let Some(median) = times.get(times.len() / 2) {
+        println!("MEDIAN {case} {w}x{h} {depth:?} {mode:?} seconds={median:.6}");
+    }
+    Ok(())
+}

crates/algo/src/pixelate.rs

@@ -5,6 +5,13 @@
 
 use photocraft_geom::Rect;
 
+mod facet;
+#[cfg(test)]
+mod reference;
+mod screen;
+#[cfg(test)]
+mod tests;
+
 use crate::fxutil::{MAXC, cell_point, luma, native, ncol, rgba, set_rgba, subtractive, via_rgb, xy};
 use crate::image::Image;
 use crate::noise::hash01;
@@ -49,6 +56,14 @@ pub(crate) fn color_halftone(src: &Image, out: Rect, ctx: &Ctx, max_radius: f32,
     };
     let w = out.width() as usize;
     let chans = if rgb { 3 } else { cc };
+    // Keep cache metadata on the stack; only the optional, bounded tables allocate.
+    let mut screens: [Option<screen::Screen>; MAXC] = std::array::from_fn(|c| {
+        if c >= chans {
+            return None;
+        }
+        let &(s, co) = rot.get(c)?;
+        screen::Screen::new(out, (ox, oy), (s, co), cell)
+    });
     let mut vals = [0.0f32; MAXC];
     for (i, px) in res.chunks_exact_mut(n).enumerate() {
         let (x, y) = (out.x0 + (i % w) as i32, out.y0 + (i / w) as i32);
@@ -68,14 +83,21 @@ pub(crate) fn color_halftone(src: &Image, out: Rect, ctx: &Ctx, max_radius: f32,
                     }
                     // Back to document space to read the cell's value.
                     let (dx, dy) = (mu * co - mv * s + ox, mu * s + mv * co + oy);
-                    let raw = value(dx, dy, c).clamp(0.0, 1.0);
-                    let ink = if sub && !rgb { raw } else { 1.0 - raw };
-                    // Dot area tracks ink: πr² = ink·cell² until dots touch, then grow to cover the corners.
-                    let touch = std::f32::consts::FRAC_PI_4;
-                    let r = if ink <= touch {
-                        cell * (ink / std::f32::consts::PI).sqrt()
-                    } else {
-                        cell * (0.5 + (ink - touch) / (1.0 - touch) * (std::f32::consts::FRAC_1_SQRT_2 - 0.5))
+                    // The five source taps and ink-to-radius mapping depend only on
+                    // the screen cell, not on the output pixel inside that cell.
+                    let radius = || {
+                        let raw = value(dx, dy, c).clamp(0.0, 1.0);
+                        let ink = if sub && !rgb { raw } else { 1.0 - raw };
+                        let touch = std::f32::consts::FRAC_PI_4;
+                        if ink <= touch {
+                            cell * (ink / std::f32::consts::PI).sqrt()
+                        } else {
+                            cell * (0.5 + (ink - touch) / (1.0 - touch) * (std::f32::consts::FRAC_1_SQRT_2 - 0.5))
+                        }
+                    };
+                    let r = match screens.get_mut(c).and_then(Option::as_mut) {
+                        Some(screen) => screen.radius(cu + di as f32, cv + dj as f32, radius),
+                        None => radius(),
                     };
                     // Tiny dots cannot cover more than their own area.
                     cov = cov.max((r - d + 0.5).clamp(0.0, 1.0).min(std::f32::consts::PI * r * r));
@@ -175,6 +197,10 @@ pub(crate) fn crystallize(src: &Image, out: Rect, ctx: &Ctx, cell_size: f32, see
 /// least-varied of its four overlapping quadrants, clumping similar colours
 /// into flat facets while keeping edges.
 pub(crate) fn facet(src: &Image, out: Rect, ctx: &Ctx) -> Vec<f32> {
+    facet::run(src, out, ctx).unwrap_or_else(|| facet_direct(src, out, ctx))
+}
+
+fn facet_direct(src: &Image, out: Rect, ctx: &Ctx) -> Vec<f32> {
     let n = src.ch;
     let r = 2i32;
     let mut res = src.crop(out);

crates/algo/src/pixelate/facet.rs

@@ -0,0 +1,94 @@
+//! Reuse each 3x3 quadrant at its four output pixels. Keep the original
+//! nine-sample addition order: reassociating the sums can change variance ties.
+use crate::fxutil::{MAXC, luma};
+use crate::{Ctx, image::Image};
+use photocraft_geom::Rect;
+
+#[derive(Clone, Copy, Default)]
+struct Quadrant {
+    variance: f32,
+    mean: [f32; MAXC],
+}
+
+pub(super) fn run(src: &Image, out: Rect, ctx: &Ctx) -> Option<Vec<f32>> {
+    let n = src.ch;
+    if n == 0 || n > MAXC || out.is_empty() {
+        return None;
+    }
+    let win = Rect::new(out.x0.checked_sub(2)?, out.y0.checked_sub(2)?, out.x1.checked_add(2)?, out.y1.checked_add(2)?);
+    if win.intersect(&src.rect) != win {
+        return None;
+    }
+    let (w, h) = (win.width() as usize, win.height() as usize);
+    let len = w.checked_mul(h)?;
+    // A normal 256px tile needs 264 KiB. Bound scratch for explicit tile overrides.
+    if len > 1024 * 1024 {
+        return None;
+    }
+    let stride = (src.rect.width() as usize).checked_mul(n)?;
+    let offset = (win.x0.abs_diff(src.rect.x0) as usize).checked_mul(n)?;
+    let row = |y: i32| -> Option<&[f32]> {
+        let start = (y.abs_diff(src.rect.y0) as usize).checked_mul(stride)?.checked_add(offset)?;
+        src.data.get(start..start.checked_add(w.checked_mul(n)?)?)
+    };
+    let mut lum = Vec::new();
+    lum.try_reserve_exact(len).ok()?;
+    for y in win.y0..win.y1 {
+        lum.extend(row(y)?.chunks_exact(n).map(|px| luma(ctx, px)));
+    }
+    let qw = w.checked_sub(2)?;
+    let mut quadrants: [Vec<Quadrant>; 3] = std::array::from_fn(|_| Vec::new());
+    for line in &mut quadrants {
+        line.try_reserve_exact(qw).ok()?;
+        line.resize(qw, Quadrant::default());
+    }
+    let build_row = |y: usize, line: &mut [Quadrant]| -> Option<()> {
+        let rows = [
+            row(win.y0.checked_add(i32::try_from(y).ok()?)?)?,
+            row(win.y0.checked_add(i32::try_from(y + 1).ok()?)?)?,
+            row(win.y0.checked_add(i32::try_from(y + 2).ok()?)?)?,
+        ];
+        for (x, q) in line.iter_mut().enumerate() {
+            let mut sum = [0.0; MAXC];
+            let (mut s, mut s2) = (0.0, 0.0);
+            for (dy, pixels) in rows.iter().enumerate() {
+                let start = x.checked_mul(n)?;
+                let samples = pixels.get(start..start.checked_add(3 * n)?)?;
+                let lights = lum.get((y + dy) * w + x..(y + dy) * w + x + 3)?;
+                for (pixel, &l) in samples.chunks_exact(n).zip(lights) {
+                    for (total, value) in sum.iter_mut().zip(pixel) {
+                        *total += *value;
+                    }
+                    s += l;
+                    s2 += l * l;
+                }
+            }
+            for v in &mut sum {
+                *v /= 9.0;
+            }
+            *q = Quadrant { variance: s2 / 9.0 - (s / 9.0).powi(2), mean: sum };
+        }
+        Some(())
+    };
+    build_row(0, &mut quadrants[0])?;
+    build_row(1, &mut quadrants[1])?;
+    let width = out.width() as usize;
+    let height = out.height() as usize;
+    let mut result = Vec::new();
+    result.try_reserve_exact(width.checked_mul(height)?.checked_mul(n)?).ok()?;
+    for y in 0..height {
+        build_row(y + 2, &mut quadrants[2])?;
+        for x in 0..width {
+            let mut best = Quadrant { variance: f32::MAX, mean: [0.0; MAXC] };
+            // Preserve quadrant order, including exact ties and non-finite variances.
+            for q in [quadrants[0].get(x)?, quadrants[0].get(x + 2)?, quadrants[2].get(x)?, quadrants[2].get(x + 2)?] {
+                if q.variance < best.variance {
+                    best = *q;
+                }
+            }
+            result.extend(best.mean.iter().take(n));
+        }
+        quadrants.rotate_left(1);
+    }
+    Some(result)
+}

crates/algo/src/pixelate/reference.rs

@@ -0,0 +1,118 @@
+//! Original scalar kernels, kept as independent test oracles.
+use super::*;
+pub(crate) fn halftone(src: &Image, out: Rect, ctx: &Ctx, max_radius: f32, angles: [f32; 4]) -> Vec<f32> {
+    let n = src.ch;
+    let cc = ncol(ctx, n);
+    let rmax = max_radius.max(1.0);
+    // Cell side such that a full-size dot just covers the cell's corners.
+    let cell = rmax * std::f32::consts::SQRT_2;
+    let sub = subtractive(ctx);
+    let rgb = via_rgb(ctx);
+    let mut res = src.crop(out);
+    let b = ctx.bounds;
+    let (ox, oy) = (b.x0 as f32, b.y0 as f32);
+    let rot: Vec<(f32, f32)> = (0..cc.max(3)).map(|c| angles[c.min(3)].to_radians().sin_cos()).collect();
+    let clip = b.intersect(&src.rect);
+    if clip.is_empty() {
+        return res;
+    }
+    // Channel value at a point averaged over five taps (for Lab, of the RGB conversion).
+    let value = |x: f32, y: f32, c: usize| -> f32 {
+        let d = cell * 0.25;
+        let mut acc = 0.0;
+        let mut px = [0.0f32; MAXC];
+        for (dx, dy) in [(0.0, 0.0), (-d, -d), (d, -d), (-d, d), (d, d)] {
+            // Clamp to the bounds so cells straddling the edge read real pixels.
+            let (ix, iy) = (((x + dx).floor() as i32).clamp(clip.x0, clip.x1 - 1), ((y + dy).floor() as i32).clamp(clip.y0, clip.y1 - 1));
+            if rgb {
+                for (k, p) in px.iter_mut().enumerate().take(n) {
+                    *p = src.get(ix, iy, k);
+                }
+                acc += rgba(ctx, &px[..n])[c];
+            } else {
+                acc += src.get(ix, iy, c);
+            }
+        }
+        acc / 5.0
+    };
+    let w = out.width() as usize;
+    let chans = if rgb { 3 } else { cc };
+    let mut vals = [0.0f32; MAXC];
+    for (i, px) in res.chunks_exact_mut(n).enumerate() {
+        let (x, y) = (out.x0 + (i % w) as i32, out.y0 + (i / w) as i32);
+        let (fx, fy) = (x as f32 + 0.5 - ox, y as f32 + 0.5 - oy);
+        for (c, v) in vals.iter_mut().enumerate().take(chans) {
+            let (s, co) = rot[c];
+            // Rotate into screen space.
+            let (u, w2) = (fx * co + fy * s, -fx * s + fy * co);
+            let (cu, cv) = ((u / cell).floor(), (w2 / cell).floor());
+            let mut cov: f32 = 0.0;
+            for dj in -1..=1 {
+                for di in -1..=1 {
+                    let (mu, mv) = ((cu + di as f32 + 0.5) * cell, (cv + dj as f32 + 0.5) * cell);
+                    let d = ((u - mu).powi(2) + (w2 - mv).powi(2)).sqrt();
+                    if d >= rmax + 0.5 {
+                        continue;
+                    }
+                    // Back to document space to read the cell's value.
+                    let (dx, dy) = (mu * co - mv * s + ox, mu * s + mv * co + oy);
+                    let raw = value(dx, dy, c).clamp(0.0, 1.0);
+                    let ink = if sub && !rgb { raw } else { 1.0 - raw };
+                    // Dot area tracks ink: πr² = ink·cell² until dots touch, then grow to cover the corners.
+                    let touch = std::f32::consts::FRAC_PI_4;
+                    let r = if ink <= touch {
+                        cell * (ink / std::f32::consts::PI).sqrt()
+                    } else {
+                        cell * (0.5 + (ink - touch) / (1.0 - touch) * (std::f32::consts::FRAC_1_SQRT_2 - 0.5))
+                    };
+                    // Tiny dots cannot cover more than their own area.
+                    cov = cov.max((r - d + 0.5).clamp(0.0, 1.0).min(std::f32::consts::PI * r * r));
+                }
+            }
+            *v = if sub && !rgb { cov } else { 1.0 - cov };
+        }
+        if rgb {
+            let a = if ctx.alpha { px[n - 1] } else { 1.0 };
+            set_rgba(ctx, px, [vals[0], vals[1], vals[2], a]);
+        } else {
+            px[..cc].copy_from_slice(&vals[..cc]);
+        }
+    }
+    res
+}
+
+pub(crate) fn facet(src: &Image, out: Rect, ctx: &Ctx) -> Vec<f32> {
+    let n = src.ch;
+    let r = 2i32;
+    let mut res = src.crop(out);
+    let mut tmp = [0.0f32; MAXC];
+    for (i, px) in res.chunks_exact_mut(n).enumerate() {
+        let (x, y) = xy(out, i);
+        let mut best = (f32::MAX, [0.0f32; MAXC]);
+        for (qx, qy) in [(-r, -r), (0, -r), (-r, 0), (0, 0)] {
+            let mut sum = [0.0f32; MAXC];
+            let (mut s, mut s2) = (0.0f32, 0.0f32);
+            for yy in y + qy..=y + qy + r {
+                for xx in x + qx..=x + qx + r {
+                    for (c, t) in tmp.iter_mut().enumerate().take(n) {
+                        *t = src.get(xx, yy, c);
+                        sum[c] += *t;
+                    }
+                    let l = luma(ctx, &tmp[..n]);
+                    s += l;
+                    s2 += l * l;
+                }
+            }
+            let cnt = ((r + 1) * (r + 1)) as f32;
+            let var = s2 / cnt - (s / cnt).powi(2);
+            if var < best.0 {
+                for v in sum.iter_mut() {
+                    *v /= cnt;
+                }
+                best = (var, sum);
+            }
+        }
+        px.copy_from_slice(&best.1[..n]);
+    }
+    res
+}

crates/algo/src/pixelate/screen.rs

@@ -0,0 +1,67 @@
+//! Tile-local memoization in rotated screen coordinates. Each cell's radius
+//! is independent of the output pixel, including the five-tap edge clamping.
+use photocraft_geom::Rect;
+
+pub(super) struct Screen {
+    x0: i32,
+    y0: i32,
+    width: usize,
+    height: usize,
+    radii: Vec<Option<f32>>,
+}
+
+impl Screen {
+    pub(super) fn new(out: Rect, origin: (f32, f32), rotation: (f32, f32), cell: f32) -> Option<Self> {
+        if out.is_empty() || !cell.is_finite() || cell < 1.0 {
+            return None;
+        }
+        let (s, co) = rotation;
+        let (mut lo_x, mut lo_y, mut hi_x, mut hi_y) = (f32::INFINITY, f32::INFINITY, f32::NEG_INFINITY, f32::NEG_INFINITY);
+        for y in [out.y0, out.y1.checked_sub(1)?] {
+            for x in [out.x0, out.x1.checked_sub(1)?] {
+                let (fx, fy) = (x as f32 + 0.5 - origin.0, y as f32 + 0.5 - origin.1);
+                let (u, v) = ((fx * co + fy * s) / cell, (-fx * s + fy * co) / cell);
+                if !u.is_finite() || !v.is_finite() || u.abs() > 4_000_000.0 || v.abs() > 4_000_000.0 {
+                    return None;
+                }
+                lo_x = lo_x.min(u);
+                lo_y = lo_y.min(v);
+                hi_x = hi_x.max(u);
+                hi_y = hi_y.max(v);
+            }
+        }
+        // Include the neighbouring dots and a rounding margin. If a key ever
+        // falls outside this box, radius() computes it directly.
+        let (x0, y0) = (lo_x.floor() as i32 - 2, lo_y.floor() as i32 - 2);
+        let width = usize::try_from(hi_x.floor() as i32 + 3 - x0).ok()?;
+        let height = usize::try_from(hi_y.floor() as i32 + 3 - y0).ok()?;
+        let len = width.checked_mul(height)?;
+        // At most 2 MiB per channel, even for an explicit oversized tile.
+        if len > 256 * 1024 {
+            return None;
+        }
+        let mut radii = Vec::new();
+        radii.try_reserve_exact(len).ok()?;
+        radii.resize(len, None);
+        Some(Self { x0, y0, width, height, radii })
+    }
+
+    pub(super) fn radius(&mut self, x: f32, y: f32, compute: impl FnOnce() -> f32) -> f32 {
+        // Keys must represent exact integers: don't alias non-finite or large
+        // coordinates onto cells whose floating-point arithmetic differs.
+        let slot = if x.is_finite() && y.is_finite() && x.abs() <= 4_000_001.0 && y.abs() <= 4_000_001.0 {
+            usize::try_from(x as i32 - self.x0)
+                .ok()
+                .zip(usize::try_from(y as i32 - self.y0).ok())
+                .filter(|&(cx, cy)| cx < self.width && cy < self.height)
+                .and_then(|(cx, cy)| cy.checked_mul(self.width)?.checked_add(cx))
+                .and_then(|i| self.radii.get_mut(i))
+        } else {
+            None
+        };
+        match slot {
+            Some(value) => *value.get_or_insert_with(compute),
+            None => compute(),
+        }
+    }
+}

crates/algo/src/pixelate/tests.rs

@@ -0,0 +1,193 @@
+use super::*;
+use photocraft_color::{ColorMode, PixelFormat, SampleType};
+use photocraft_raster::Surface;
+
+fn fixture(rect: Rect, ctx: &Ctx, pattern: u32) -> Image {
+    let n = PixelFormat::new(ctx.mode, SampleType::F32, ctx.alpha).channels();
+    let mut image = Image::new(rect, n);
+    for (i, pixel) in image.data.chunks_exact_mut(n).enumerate() {
+        let (x, y) = xy(rect, i);
+        let seed = (x as u32).wrapping_mul(0x9e37_79b9) ^ (y as u32).wrapping_mul(0x85eb_ca6b);
+        for (c, sample) in pixel.iter_mut().enumerate() {
+            *sample = match pattern {
+                0 => 0.5,
+                1 => {
+                    if x + y > -4 {
+                        0.9
+                    } else {
+                        0.1
+                    }
+                }
+                2 => (seed.rotate_left(c as u32 * 7) % 65536) as f32 / 65535.0,
+                3 => (seed.rotate_left(c as u32 * 7) % 65536) as f32 / 16384.0 - 1.0,
+                _ => match seed.wrapping_add(c as u32) % 7 {
+                    0 => f32::NAN,
+                    1 => f32::INFINITY,
+                    2 => f32::NEG_INFINITY,
+                    _ => 0.5,
+                },
+            };
+        }
+        if ctx.alpha && seed.is_multiple_of(5) {
+            pixel[n - 1] = 0.0;
+        }
+    }
+    image
+}
+
+fn assert_bits(got: &[f32], want: &[f32]) {
+    assert_eq!(got.len(), want.len());
+    for (i, (a, b)) in got.iter().zip(want).enumerate() {
+        assert_eq!(a.to_bits(), b.to_bits(), "sample {i}: {a} != {b}");
+    }
+}
+
+#[test]
+fn facet_reuse_matches_scalar_in_every_model() {
+    let out = Rect::new(-17, -9, 20, 18);
+    for mode in [
+        ColorMode::Rgb,
+        ColorMode::Grayscale,
+        ColorMode::Cmyk,
+        ColorMode::Lab,
+        ColorMode::Bitmap,
+        ColorMode::Duotone,
+        ColorMode::Indexed,
+        ColorMode::Multichannel,
+    ] {
+        for alpha in [false, true] {
+            let ctx = Ctx { bounds: out, mode, alpha };
+            for pattern in 0..5 {
+                let src = fixture(out.inflate(3), &ctx, pattern);
+                assert_bits(&facet(&src, out, &ctx), &reference::facet(&src, out, &ctx));
+            }
+        }
+    }
+}
+
+#[test]
+fn facet_partial_halo_and_thin_regions_match_scalar() {
+    for out in [Rect::new(-3, 7, -2, 8), Rect::new(-5, 0, 30, 1), Rect::new(9, -8, 10, 29)] {
+        let ctx = Ctx { bounds: out, mode: ColorMode::Rgb, alpha: true };
+        for halo in [0, 1, 2, 3] {
+            let src = fixture(out.inflate(halo), &ctx, 2);
+            assert_bits(&facet(&src, out, &ctx), &reference::facet(&src, out, &ctx));
+        }
+    }
+}
+
+#[test]
+fn halftone_cache_matches_scalar_screen_angles_radii_and_models() {
+    let out = Rect::new(-17, -9, 20, 18);
+    for mode in [
+        ColorMode::Rgb,
+        ColorMode::Grayscale,
+        ColorMode::Cmyk,
+        ColorMode::Lab,
+        ColorMode::Bitmap,
+        ColorMode::Duotone,
+        ColorMode::Indexed,
+        ColorMode::Multichannel,
+    ] {
+        for alpha in [false, true] {
+            let ctx = Ctx { bounds: Rect::new(-23, -15, 29, 31), mode, alpha };
+            for pattern in 0..5 {
+                let src = fixture(ctx.bounds, &ctx, pattern);
+                for radius in [1.0, 2.0, 3.75, 8.0, 64.0] {
+                    for angles in [[108.0, 162.0, 90.0, 45.0], [0.0, -17.5, 359.99, 720.0]] {
+                        assert_bits(&color_halftone(&src, out, &ctx, radius, angles), &reference::halftone(&src, out, &ctx, radius, angles));
+                    }
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn halftone_unusual_coordinates_and_params_keep_direct_semantics() {
+    for out in [Rect::new(16_777_211, -3, 16_777_220, 4), Rect::new(-3, 4, -2, 5)] {
+        let ctx = Ctx { bounds: out, mode: ColorMode::Rgb, alpha: true };
+        let src = fixture(out.inflate(2), &ctx, 2);
+        for radius in [f32::NAN, f32::INFINITY, -10.0, 1.0] {
+            for angles in [[f32::NAN, 0.0, f32::INFINITY, 45.0], [0.0; 4]] {
+                assert_bits(&color_halftone(&src, out, &ctx, radius, angles), &reference::halftone(&src, out, &ctx, radius, angles));
+            }
+        }
+    }
+    assert!(screen::Screen::new(Rect::new(0, 0, 100_000, 100_000), (0.0, 0.0), (0.0, 1.0), 1.0).is_none());
+    assert!(screen::Screen::new(Rect::new(8_000_000, 0, 8_000_010, 10), (0.0, 0.0), (0.0, 1.0), 1.0).is_none());
+    let mut screen = screen::Screen::new(Rect::new(0, 0, 1, 1), (0.0, 0.0), (0.0, 1.0), 1.0).unwrap();
+    assert_eq!(screen.radius(-100.0, 0.0, || 3.0), 3.0);
+    assert_eq!(screen.radius(f32::NAN, 0.0, || 4.0), 4.0);
+    assert_eq!(screen.radius(0.0, 0.0, || 5.0), 5.0);
+    assert_eq!(screen.radius(0.0, 0.0, || panic!("cached dot recomputed")), 5.0);
+}
+
+#[test]
+fn tiled_filters_match_scalar_at_all_depths_with_soft_selection() {
+    use crate::{FilterParams, apply_tiled};
+    let bounds = Rect::new(-17, -9, 26, 24);
+    for mode in [ColorMode::Rgb, ColorMode::Grayscale, ColorMode::Cmyk, ColorMode::Lab] {
+        for depth in [SampleType::U8, SampleType::U16, SampleType::F32] {
+            for alpha in [false, true] {
+                let ctx = Ctx { bounds, mode, alpha };
+                let src = fixture(bounds, &ctx, 2);
+                let mut surface = Surface::new(PixelFormat::new(mode, depth, alpha));
+                surface.write_region(bounds, &src.data);
+                let mut selection = Surface::new(PixelFormat::new(ColorMode::Grayscale, SampleType::U8, false));
+                let coverage: Vec<_> = (0..bounds.width() * bounds.height()).map(|i| (i % 7) as f32 / 6.0).collect();
+                selection.write_region(bounds, &coverage);
+                for params in [FilterParams::Facet, FilterParams::ColorHalftone { max_radius: 3.75, angles: [108.0, 162.0, 90.0, 45.0] }] {
+                    for extent in [None, Some(bounds)] {
+                        assert!(matches!(params.halo(), crate::Halo::Radius(_)));
+                        let halo = match params.halo() {
+                            crate::Halo::Radius(r) => r,
+                            _ => 0,
+                        };
+                        let image = match extent {
+                            Some(e) => Image::read_clamped(&surface, bounds.inflate(halo), e),
+                            None => Image::read(&surface, bounds.inflate(halo)),
+                        };
+                        let mut expected_data = match params {
+                            FilterParams::Facet => reference::facet(&image, bounds, &ctx),
+                            FilterParams::ColorHalftone { max_radius, angles } => reference::halftone(&image, bounds, &ctx, max_radius, angles),
+                            _ => Vec::new(),
+                        };
+                        crate::mix_selection(&mut expected_data, bounds, &selection, &image);
+                        let mut expected = surface.clone();
+                        expected.write_region(bounds, &expected_data);
+                        for tile in [1, 7, 17, 256] {
+                            let got = apply_tiled(&surface, &params, bounds, bounds, Some(&selection), tile, extent);
+                            assert_bits(&got.read_region(bounds), &expected.read_region(bounds));
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+#[test]
+fn extended_channels_and_large_screen_offsets_match_scalar() {
+    let out = Rect::new(16_777_211, -3, 16_777_220, 4);
+    let ctx = Ctx { bounds: Rect::new(0, -10, out.x1 + 10, 10), mode: ColorMode::Multichannel, alpha: true };
+    for n in 5..=MAXC {
+        let mut src = Image::new(out.inflate(3), n);
+        for (i, sample) in src.data.iter_mut().enumerate() {
+            *sample = (i * 73 % 257) as f32 / 256.0;
+        }
+        assert_bits(&facet(&src, out, &ctx), &reference::facet(&src, out, &ctx));
+        assert_bits(&color_halftone(&src, out, &ctx, 1.0, [0.0; 4]), &reference::halftone(&src, out, &ctx, 1.0, [0.0; 4]));
+    }
+}
+
+#[test]
+fn facet_scratch_rejects_large_or_incomplete_storage() {
+    let out = Rect::new(0, 0, 10_000, 10_000);
+    let ctx = Ctx { bounds: out, mode: ColorMode::Rgb, alpha: true };
+    let src = Image { rect: out.inflate(3), ch: 4, data: Vec::new() };
+    assert!(facet::run(&src, out, &ctx).is_none());
+    let out = Rect::new(0, 0, 2, 2);
+    let src = Image { rect: out.inflate(3), ch: 4, data: Vec::new() };
+    assert!(facet::run(&src, out, &ctx).is_none());
+}

docs/pixelate-performance.json

@@ -0,0 +1,278 @@
+[
+  {
+    "mode": "rgb",
+    "case": "facet",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        0.686463,
+        0.684428,
+        0.681198
+      ],
+      "after": [
+        0.311512,
+        0.309932,
+        0.31487
+      ]
+    },
+    "before": 0.684428,
+    "after": 0.311512,
+    "sha256": "d629acb2a4afda38dc5976ac8e552e44e419abb001a4b9a3ec7b42e289c4fac3",
+    "speedup": 2.1971160019517706
+  },
+  {
+    "mode": "rgb",
+    "case": "halftone",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        1.135315,
+        1.13726,
+        1.142235
+      ],
+      "after": [
+        0.597161,
+        0.612252,
+        0.603813
+      ]
+    },
+    "before": 1.13726,
+    "after": 0.603813,
+    "sha256": "d8c056438d7e9c8a954c8913bb7185deaa154618c0eeb9c36a65f87148744b28",
+    "speedup": 1.8834639201209644
+  },
+  {
+    "mode": "rgb",
+    "case": "facet",
+    "width": 7200,
+    "height": 5000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        1.027543,
+        1.029425,
+        1.026628
+      ],
+      "after": [
+        0.454462,
+        0.450456,
+        0.476923
+      ]
+    },
+    "before": 1.027543,
+    "after": 0.454462,
+    "sha256": "cbe05cee9bd9ac15e908dc48f88668b36983a8159c326beed53ee5a84f9c8e85",
+    "speedup": 2.2610097213848466
+  },
+  {
+    "mode": "rgb",
+    "case": "halftone",
+    "width": 7200,
+    "height": 5000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        1.707774,
+        1.719981,
+        1.693232
+      ],
+      "after": [
+        0.897247,
+        0.905924,
+        0.898006
+      ]
+    },
+    "before": 1.707774,
+    "after": 0.898006,
+    "sha256": "bb7810a5ae26c0860ce502b3524e6ed51c75359afb230f4195fdb065f90880c2",
+    "speedup": 1.9017400774605069
+  },
+  {
+    "mode": "rgb",
+    "case": "facet",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u16",
+    "times": {
+      "before": [
+        0.765632,
+        0.755045,
+        0.756924
+      ],
+      "after": [
+        0.368377,
+        0.370144,
+        0.375268
+      ]
+    },
+    "before": 0.756924,
+    "after": 0.370144,
+    "sha256": "8090fc2be1997d904bd8b23dbb85466b318cb25d5ed85b567c506124dd721c57",
+    "speedup": 2.044944670182416
+  },
+  {
+    "mode": "rgb",
+    "case": "halftone",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u16",
+    "times": {
+      "before": [
+        1.203186,
+        1.189726,
+        1.198921
+      ],
+      "after": [
+        0.6647,
+        0.662166,
+        0.665923
+      ]
+    },
+    "before": 1.198921,
+    "after": 0.6647,
+    "sha256": "c176f05bc3e98762add441c38c8cddca4fdbcf720c54473f70e2656614d4a018",
+    "speedup": 1.8037024221453286
+  },
+  {
+    "mode": "rgb",
+    "case": "facet",
+    "width": 6000,
+    "height": 4000,
+    "depth": "f32",
+    "times": {
+      "before": [
+        0.806876,
+        0.807626,
+        0.804216
+      ],
+      "after": [
+        0.429701,
+        0.435696,
+        0.441177
+      ]
+    },
+    "before": 0.806876,
+    "after": 0.435696,
+    "sha256": "049ea951c124c77a3090751be168a07949a81340eee7c622b28603eefb5f2b3d",
+    "speedup": 1.8519242774778744
+  },
+  {
+    "mode": "rgb",
+    "case": "halftone",
+    "width": 6000,
+    "height": 4000,
+    "depth": "f32",
+    "times": {
+      "before": [
+        1.283232,
+        1.272004,
+        1.266286
+      ],
+      "after": [
+        0.730128,
+        0.722809,
+        0.72555
+      ]
+    },
+    "before": 1.272004,
+    "after": 0.72555,
+    "sha256": "d91c9f6123ea5ae9842e4761ca6c49a27371cc6101813ada6a2806e92ef61e9b",
+    "speedup": 1.753158293708221
+  },
+  {
+    "mode": "rgb",
+    "case": "halftone-small",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        1.417763,
+        1.418794,
+        1.407005
+      ],
+      "after": [
+        0.720082,
+        0.716397,
+        0.71179
+      ]
+    },
+    "before": 1.417763,
+    "after": 0.716397,
+    "sha256": "d92f681ce4285bebda335bc3ed2288e8f1f468216f208fc5eb57654691b04ea0",
+    "speedup": 1.979018616772544
+  },
+  {
+    "mode": "rgb",
+    "case": "halftone-large",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        1.138943,
+        1.138331,
+        1.140677
+      ],
+      "after": [
+        0.583799,
+        0.587533,
+        0.592007
+      ]
+    },
+    "before": 1.138943,
+    "after": 0.587533,
+    "sha256": "b4b7835d3eff7548e2f71a31c4f87e728bf1754d2d680269bedd7a428091916e",
+    "speedup": 1.9385174960385205
+  },
+  {
+    "mode": "cmyk",
+    "case": "facet",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        7.888404,
+        7.841185,
+        7.897717
+      ],
+      "after": [
+        0.574157,
+        0.569647,
+        0.576679
+      ]
+    },
+    "before": 7.888404,
+    "after": 0.574157,
+    "sha256": "055e8618c5554839c6c3566a8373eca86afaaf13ac2010454828352cdf7ca425",
+    "speedup": 13.739106202658855
+  },
+  {
+    "mode": "lab",
+    "case": "halftone",
+    "width": 6000,
+    "height": 4000,
+    "depth": "u8",
+    "times": {
+      "before": [
+        3.576011,
+        3.590308,
+        3.588703
+      ],
+      "after": [
+        0.719498,
+        0.71831,
+        0.714674
+      ]
+    },
+    "before": 3.588703,
+    "after": 0.71831,
+    "sha256": "f52df18bebfbed6f7a8d7b2071c0c5221629f64579ef58571893d7b50301f2b3",
+    "speedup": 4.996036530188916
+  }
+]

docs/pixelate-performance.md

@@ -0,0 +1,55 @@
+# Facet and Color Halftone performance
+
+2026-10-10. Two changes remove repeated work from the existing Pixelate filters:
+
+- **Facet:** calculate luma once per source pixel and calculate each 3x3 quadrant once. A three-row ring shares that quadrant with its four output pixels. The nine-sample addition order, variance calculation, quadrant order and strict tie comparison stay the same.
+- **Color Halftone:** memoize the five-tap ink calculation and dot radius in the rotated screen grid, separately for each channel. Pixels still use the same candidate dots, distances and anti-alias coverage. Lab conversion, CMYK ink direction and alpha handling stay the same.
+
+Facet uses about 291 KiB of additional scratch for a 256x256 output tile. The luma plane is capped at 1,048,576 samples. Halftone tables are capped at 2 MiB per channel. Scratch allocations are fallible; an incomplete halo, oversized table or unsuitable screen coordinates uses the original calculation. The default tiled filter path stays parallel on native targets and single-threaded on wasm.
+
+## Full filter timings
+
+EC2 m7i.4xlarge, 16 vCPUs (Intel Xeon Platinum 8488C), Ubuntu 24.04, rustc 1.99.0, default release profile, `RAYON_NUM_THREADS=16`. Base: `7722172585a01cbdb93c06f0f5ff2634fcb17999`. The changed algorithms are in `1bc3a9f1`. Before and after were built in separate target directories; no other builds or benchmarks ran during measurement.
+
+Each entry is the median of three alternating pairs: before/after, after/before, before/after. The timer covers `apply_in`, including tile reads, filtering, assembly, depth conversion, writeback and pruning. Fixture construction, the small warmup and output serialization are outside the timer. The fixture is a dense deterministic native-color image with varying alpha, negative document origin and transparent pixels. These are CPU filter application times; window presentation and undo bookkeeping are outside this measurement.
+
+| Filter | Size | Model + alpha | Depth | Before (ms) | After (ms) | Speedup |
+|---|---:|---|---|---:|---:|---:|
+| Facet | 24 MP | RGB | u8 | 684.4 | 311.5 | 2.20x |
+| Halftone r8 | 24 MP | RGB | u8 | 1137.3 | 603.8 | 1.88x |
+| Facet | 36 MP | RGB | u8 | 1027.5 | 454.5 | 2.26x |
+| Halftone r8 | 36 MP | RGB | u8 | 1707.8 | 898.0 | 1.90x |
+| Facet | 24 MP | RGB | u16 | 756.9 | 370.1 | 2.04x |
+| Halftone r8 | 24 MP | RGB | u16 | 1198.9 | 664.7 | 1.80x |
+| Facet | 24 MP | RGB | f32 | 806.9 | 435.7 | 1.85x |
+| Halftone r8 | 24 MP | RGB | f32 | 1272.0 | 725.6 | 1.75x |
+| Halftone r2 | 24 MP | RGB | u8 | 1417.8 | 716.4 | 1.98x |
+| Halftone r64 | 24 MP | RGB | u8 | 1138.9 | 587.5 | 1.94x |
+| Facet | 24 MP | CMYK | u8 | 7888.4 | 574.2 | 13.74x |
+| Halftone r8 | 24 MP | LAB | u8 | 3588.7 | 718.3 | 5.00x |
+
+Every pair was compared with `cmp` over the complete output, serialized as little-endian f32 sample bits in tile order. All 36 comparisons matched exactly, including the 16-bit and float cases. [Raw three-run timings and output SHA-256 hashes](pixelate-performance.json) are committed alongside this table.
+
+## Reproduce
+
+Build the same example against the base and the changed tree. Keep their Cargo target directories separate:
+
+```sh
+git worktree add --detach ../photocraft-pixelate-base 7722172585a01cbdb93c06f0f5ff2634fcb17999
+cp crates/algo/examples/bench_pixelate.rs ../photocraft-pixelate-base/crates/algo/examples/
+(cd ../photocraft-pixelate-base && CARGO_TARGET_DIR=target cargo build --release -p photocraft-algo --example bench_pixelate)
+CARGO_TARGET_DIR=target cargo build --release -p photocraft-algo --example bench_pixelate
+
+RAYON_NUM_THREADS=16 ../photocraft-pixelate-base/target/release/examples/bench_pixelate facet 6000 4000 u8 1 before.bin rgb
+RAYON_NUM_THREADS=16 target/release/examples/bench_pixelate facet 6000 4000 u8 1 after.bin rgb
+cmp before.bin after.bin
+sha256sum before.bin after.bin
+```
+
+Repeat each pair three times, reversing the order for the middle pair. Use `7200 5000` for 36 MP; `u16` or `f32` for the other depths; `halftone`, `halftone-small` or `halftone-large` for the screen cases; and `cmyk` or `lab` for the other color models. The example prints the timed operation separately from file writing.
+
+## Correctness and repository checks
+
+The independent scalar oracles retain the original kernels. Tests compare every sample bit across all eight color models, alpha on/off, constant colors and edges, random samples, HDR/negative values, non-finite samples, thin regions, partial halos, extreme screen offsets, extended channel counts and scratch rejection. Tiled surface tests add 8/16/32-bit storage, soft selections, tile sizes 1/7/17/256 and both transparent and repeated edges.
+
+Passed: 401 algorithm tests; clippy with warnings denied; formatting; dependency layering; the full `cargo xtask wasm` set; `cargo xtask perf --quick` on llvmpipe; and scorecard generation. No command, UI or format contracts changed.

docs/roadmap.md

@@ -204,6 +204,16 @@ Bounded halo tiles reduce the large-mask workspace allocation bound from 192/288
 to about 115/163 MB, including output and 16 workers' scratch, excluding the source.
 See [method and results](octagonal-mask-performance.md).
 
+2026-10-10: Pixelate's Facet shares exact 3x3 quadrant statistics and Color Halftone
+shares dot radii within each rotated screen grid. Full tiled CPU application at 24 MP
+on a 16-vCPU m7i.4xlarge improves Facet from 0.684 s to 0.312 s in RGB8
+(2.20x), and from 7.888 s to 0.574 s in CMYK8 (13.74x).
+Color Halftone improves from 1.137 s to 0.604 s in RGB8 (1.88x)
+and from 3.589 s to 0.718 s in Lab8 (5.00x).
+Three alternating paired release runs also cover 36 MP RGB, 16/32-bit storage,
+and screen radii 2/8/64. All 36 complete output comparisons match exactly.
+See [method, raw timings and checks](pixelate-performance.md).
+
 ### Where we're going (priority order)
 
 1. **Ship the fixes:** cut 0.2.1 once the current batch lands, so users validate them.