From 4d240cf6163a2827df8332a5a455ea7e9c495dea Mon Sep 17 00:00:00 2001 From: Matthew James Sivak Date: Tue, 13 Oct 2026 14:10:21 +0200 Subject: [PATCH] LR-EDIT-EFFECTS-DEHAZE: raster::morph square min / max filters, Image::transpose Split out of #1286 (Dehaze V3) so it can be reviewed and land on its own. `min_filter` / `max_filter` are separable van Herk / Gil-Werman filters: three comparisons per pixel per pass whatever the radius, edges clamped, NaNs skipped like f32::min, a short buffer or empty plane returned unchanged. `Image::transpose` (row-parallel) does the column pass. Co-Authored-By: Claude Opus 5.5 --- crates/raster/src/lib.rs | 29 ++++++ crates/raster/src/morph.rs | 187 +++++++++++++++++++++++++++++++++++++ 2 files changed, 216 insertions(+) create mode 100644 crates/raster/src/morph.rs diff --git a/crates/raster/src/lib.rs b/crates/raster/src/lib.rs index 41ae75af8..08b30cfa7 100644 --- a/crates/raster/src/lib.rs +++ b/crates/raster/src/lib.rs @@ -4,12 +4,14 @@ //! `Plane = Image`). //! - Resampling ([`resample`]): separable Lanczos-3 / Mitchell / box, rayon-parallel. //! - Blur ([`blur`]): separable Gaussian via repeated box filters (O(1) per pixel in the radius). +//! - Morphology ([`morph`]): square min / max filters (van Herk / Gil-Werman, O(1) per pixel in the radius). //! - [`Histogram`] of display-encoded RGB + luminance. #![forbid(unsafe_code)] #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] pub mod blur; pub mod histogram; +pub mod morph; pub mod resample; pub use histogram::{HdrBins, Histogram}; @@ -98,6 +100,21 @@ impl Image { pub fn is_empty(&self) -> bool { self.data.is_empty() } + /// The image mirrored about its main diagonal (`width` and `height` swap; row-parallel). A + /// buffer shorter than its size reads as `T::default()` where it runs out. + pub fn transpose(&self) -> Self + where + T: Send + Sync, + { + let (w, h) = (self.width, self.height); + let mut t = Image::::new(h, w); + par_rows(&mut t.data, h, |x, row| { + for (y, v) in row.iter_mut().enumerate() { + *v = self.data.get(y * w + x).copied().unwrap_or_default(); + } + }); + t + } /// [`Self::map`] into the same buffer (no second image in memory). pub fn map_in_place(&mut self, f: impl Fn(T) -> T + Sync + Send) where @@ -258,6 +275,18 @@ mod tests { use super::*; use lightcraft_geom::Orientation; + #[test] + fn transpose_swaps_the_axes() { + let p = Plane::from_fn(5, 3, |x, y| (y * 10 + x) as f32); + let t = p.transpose(); + assert_eq!((t.width, t.height), (3, 5)); + assert!((0..3).all(|y| (0..5).all(|x| t.get(y, x) == p.get(x, y)))); + assert_eq!(t.transpose(), p); + // a buffer shorter than its size reads as zero where it runs out, never a panic + let short = Plane { width: 3, height: 2, data: vec![1.0; 4] }; + assert_eq!(short.transpose().data, [1.0, 1.0, 1.0, 0.0, 1.0, 0.0]); + } + #[test] fn crop_and_rotate() { let img = Image::::from_fn(4, 3, |x, y| (y * 10 + x) as u32); diff --git a/crates/raster/src/morph.rs b/crates/raster/src/morph.rs new file mode 100644 index 000000000..e428f91c2 --- /dev/null +++ b/crates/raster/src/morph.rs @@ -0,0 +1,187 @@ +//! Grey-scale morphology: square min (erosion) and max (dilation) filters. +//! +//! Separable van Herk / Gil-Werman filters: each 1-D pass splits the (edge-extended) line into +//! blocks of the window's width and keeps a running min from the left and from the right within +//! each block; any window then spans at most two blocks, so its min is the min of one suffix and +//! one prefix value. That is three comparisons per pixel per pass, whatever the radius. +//! +//! Edges are clamped (the edge pixel repeats), which for a min or max is the same as clipping the +//! window to the image. NaNs are skipped like [`f32::min`] does: a window yields NaN only when every +//! value in it is NaN. + +use crate::Plane; + +/// Square min filter (erosion) of radius `r`: each output pixel is the minimum over the +/// `(2r+1)²` window around it, clipped to the image. `r == 0` returns a copy. +pub fn min_filter(p: &Plane, r: usize) -> Plane { + filter(p, r, f32::min) +} + +/// Square max filter (dilation) of radius `r` (see [`min_filter`]). +pub fn max_filter(p: &Plane, r: usize) -> Plane { + filter(p, r, f32::max) +} + +fn filter(p: &Plane, r: usize, pick: fn(f32, f32) -> f32) -> Plane { + let (w, h) = (p.width, p.height); + if r == 0 || p.data.is_empty() || p.data.len() != w.saturating_mul(h) { + return p.clone(); + } + let mut out = p.clone(); + rows(&mut out, r, pick); + let mut t = out.transpose(); + rows(&mut t, r, pick); + t.transpose() +} + +/// [`line`] over every row of `p`, in place (row-parallel, one scratch per worker). +fn rows(p: &mut Plane, r: usize, pick: fn(f32, f32) -> f32) { + let w = p.width; + if w == 0 { + return; + } + #[cfg(feature = "parallel")] + { + use rayon::prelude::*; + p.data.par_chunks_mut(w).for_each_init(Scratch::default, |s, row| line(row, r, pick, s)); + } + #[cfg(not(feature = "parallel"))] + { + let mut s = Scratch::default(); + p.data.chunks_mut(w).for_each(|row| line(row, r, pick, &mut s)); + } +} + +#[derive(Default)] +struct Scratch { + ext: Vec, + pre: Vec, + suf: Vec, +} + +/// One 1-D pass over `row`, in place. +fn line(row: &mut [f32], r: usize, pick: fn(f32, f32) -> f32, s: &mut Scratch) { + let n = row.len(); + if n == 0 { + return; + } + // a radius past the line covers it all; capping keeps the scratch bounded + let r = r.min(n - 1); + if r == 0 { + return; + } + let k = 2 * r + 1; + // the line with `r` copies of each edge pixel on either side + s.ext.clear(); + s.ext.extend(std::iter::repeat_n(row[0], r)); + s.ext.extend_from_slice(row); + s.ext.extend(std::iter::repeat_n(row[n - 1], r)); + let m = s.ext.len(); + s.pre.clear(); + s.pre.resize(m, 0.0); + s.suf.clear(); + s.suf.resize(m, 0.0); + for (i, &v) in s.ext.iter().enumerate() { + s.pre[i] = if i % k == 0 { v } else { pick(s.pre[i - 1], v) }; + } + for i in (0..m).rev() { + let v = s.ext[i]; + s.suf[i] = if i + 1 == m || (i + 1) % k == 0 { v } else { pick(s.suf[i + 1], v) }; + } + // window [x, x + 2r] of `ext` is output x + for (x, o) in row.iter_mut().enumerate() { + *o = pick(s.suf[x], s.pre[x + 2 * r]); + } +} + +#[cfg(test)] +mod tests { + use super::*; + + /// The definition: min (or max) over the clipped window. + fn naive(p: &Plane, r: usize, max: bool) -> Plane { + let (w, h) = (p.width, p.height); + let pick = |a: f32, b: f32| if max { a.max(b) } else { a.min(b) }; + Plane::from_fn(w, h, |x, y| { + let mut m = p.get(x, y); + for yy in y.saturating_sub(r)..=(y + r).min(h - 1) { + for xx in x.saturating_sub(r)..=(x + r).min(w - 1) { + m = pick(m, p.get(xx, yy)); + } + } + m + }) + } + + /// splitmix64-ish deterministic noise. + fn noise(w: usize, h: usize, seed: u64) -> Plane { + Plane::from_fn(w, h, |x, y| { + let mut z = seed.wrapping_add((y * w + x) as u64).wrapping_mul(0x9e37_79b9_7f4a_7c15); + z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9); + z ^= z >> 31; + (z % 10_000) as f32 / 10_000.0 + }) + } + + #[test] + fn matches_the_naive_filter() { + for (w, h) in [(1, 1), (1, 9), (9, 1), (7, 5), (31, 17), (64, 40)] { + for r in [1, 2, 3, 5, 8, 20, 100] { + let p = noise(w, h, (w * 131 + h * 7 + r) as u64); + assert_eq!(min_filter(&p, r), naive(&p, r, false), "min {w}x{h} r={r}"); + assert_eq!(max_filter(&p, r), naive(&p, r, true), "max {w}x{h} r={r}"); + } + } + } + + #[test] + fn radius_zero_and_empty_are_copies() { + let p = noise(6, 4, 1); + assert_eq!(min_filter(&p, 0), p); + let e = Plane::new(0, 0); + assert_eq!(min_filter(&e, 3), e); + let e = Plane::new(5, 0); + assert_eq!(max_filter(&e, 3), e); + } + + #[test] + fn radius_past_the_image_is_the_global_extreme() { + let p = noise(13, 7, 3); + let lo = p.data.iter().copied().fold(f32::INFINITY, f32::min); + let hi = p.data.iter().copied().fold(f32::NEG_INFINITY, f32::max); + assert!(min_filter(&p, 1_000_000).data.iter().all(|&v| v == lo)); + assert!(max_filter(&p, usize::MAX).data.iter().all(|&v| v == hi)); + } + + #[test] + fn a_dark_dot_grows_to_the_window() { + let mut p = Plane::filled(11, 11, 1.0); + p.set(5, 5, 0.0); + let q = min_filter(&p, 2); + for y in 0..11 { + for x in 0..11 { + let inside = (3..=7).contains(&x) && (3..=7).contains(&y); + assert_eq!(q.get(x, y), if inside { 0.0 } else { 1.0 }, "({x}, {y})"); + } + } + } + + #[test] + fn nans_are_skipped() { + let mut p = Plane::filled(5, 5, 0.5); + p.set(2, 2, f32::NAN); + p.set(0, 0, 0.25); + let q = min_filter(&p, 1); + assert_eq!(q.get(2, 2), 0.5); + assert_eq!(q.get(1, 1), 0.25); + let all = Plane::filled(3, 3, f32::NAN); + assert!(min_filter(&all, 1).data.iter().all(|v| v.is_nan())); + } + + #[test] + fn a_short_buffer_is_returned_unchanged() { + // a plane whose data doesn't fill its size must not panic + let p = Plane { width: 4, height: 4, data: vec![1.0; 3] }; + assert_eq!(min_filter(&p, 2), p); + } +}