Merge PR #1976: GUM-3201: Reduce guided-filter coefficient allocations

Cargo.lock

@@ -3403,6 +3403,7 @@ dependencies = [
 name = "lightcraft-pipeline"
 version = "0.11.0"
 dependencies = [
+ "dhat",
  "image",
  "lightcraft-color",
  "lightcraft-develop",

contributors/people.toml

@@ -51,6 +51,11 @@ sync_github_name = true
 display_name = "Netter"
 sync_github_name = false
 
+[people.ternarybits]
+real_name = "Ted Mao"
+display_name = "ternarybits"
+sync_github_name = false
+
 [people.nledford]
 real_name = "Nathaniel Ledford"
 

crates/pipeline/Cargo.toml

@@ -19,6 +19,7 @@ rayon = { workspace = true }
 serde = { workspace = true }
 
 [dev-dependencies]
+dhat = { workspace = true }
 lightcraft-scenes = { workspace = true }
 proptest = { workspace = true }
 image = { workspace = true }

crates/pipeline/examples/guided_memory.rs

@@ -0,0 +1,84 @@
+//! Isolated fast-guided-filter benchmark (procedural log luminance; no decoded RGB or tone median).
+//! `RAYON_NUM_THREADS=1 cargo run --release -p lightcraft-pipeline --example guided_memory -- new`
+//! Replace `new` with `old` for the materialized-coefficient reference. Optional positional arguments:
+//! width height sigma iterations (defaults: 6000 4000 120 5). Append `--heap` for one dhat-counted run.
+#![forbid(unsafe_code)]
+#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
+
+use lightcraft_pipeline::local::{BASE_EPS, guided, guided_fast, guided_fast_step};
+use lightcraft_raster::blur::gaussian;
+use lightcraft_raster::resample::{Filter, resize};
+use lightcraft_raster::{Plane, par_join, par_rows};
+use std::hint::black_box;
+use std::time::Instant;
+
+#[global_allocator]
+static ALLOC: dhat::Alloc = dhat::Alloc;
+
+/// Previous algorithm; kept here only to make before/after measurements reproducible.
+fn materialized(p: &Plane, sigma: f32) -> Plane {
+    let s = guided_fast_step(sigma);
+    if s <= 1 {
+        return guided(p, sigma, BASE_EPS);
+    }
+    let lo = resize(p, p.width.div_ceil(s).max(1), p.height.div_ceil(s).max(1), Filter::Box);
+    let (a, b) = coefficients(&lo, sigma / s as f32);
+    let (a, b) = (resize(&a, p.width, p.height, Filter::Bilinear), resize(&b, p.width, p.height, Filter::Bilinear));
+    let mut q = Plane::new(p.width, p.height);
+    par_rows(&mut q.data, p.width, |y, row| {
+        for ((v, pv), (av, bv)) in row.iter_mut().zip(p.row(y)).zip(a.row(y).iter().zip(b.row(y))) {
+            *v = av * pv + bv;
+        }
+    });
+    q
+}
+
+fn coefficients(lo: &Plane, sg: f32) -> (Plane, Plane) {
+    let (mean, corr) = par_join(|| gaussian(lo, sg), || gaussian(&lo.map(|v| v * v), sg));
+    let a = corr.zip_map(&mean, |c, m| {
+        let v = (c - m * m).max(0.0);
+        v / (v + BASE_EPS)
+    });
+    let b = mean.zip_map(&a, |m, a| m - a * m);
+    par_join(|| gaussian(&a, sg), || gaussian(&b, sg))
+}
+
+fn main() -> Result<(), Box<dyn std::error::Error>> {
+    let args: Vec<_> = std::env::args().skip(1).collect();
+    let mode = args.first().map(String::as_str).unwrap_or("new");
+    if !matches!(mode, "old" | "new") {
+        return Err("mode must be old or new".into());
+    }
+    let parse = |i: usize, default: usize| -> Result<usize, Box<dyn std::error::Error>> {
+        Ok(args.get(i).filter(|v| !v.starts_with("--")).map(|v| v.parse()).transpose()?.unwrap_or(default))
+    };
+    let (w, h, sigma, n) = (parse(1, 6000)?, parse(2, 4000)?, parse(3, 120)? as f32, parse(4, 5)?.max(1));
+    if w == 0 || h == 0 || w.checked_mul(h).is_none_or(|pixels| pixels > 100_000_000) {
+        return Err("dimensions must be positive and at most 100 MP".into());
+    }
+    let heap = args.iter().any(|a| a == "--heap");
+    let p = Plane::from_fn(w, h, |x, y| ((x * 7919 + y * 104729) % 97) as f32 / 97.0 - 4.0 + if x > w / 2 { 2.0 } else { 0.0 });
+    let filter = || if mode == "old" { materialized(&p, sigma) } else { guided_fast(&p, sigma, BASE_EPS) };
+    // Initialize Rayon and warm code paths before counting. The source is excluded from heap stats.
+    black_box(filter());
+    let profiler = heap.then(|| dhat::Profiler::builder().testing().build());
+    let mut best = f64::MAX;
+    for _ in 0..if heap { 1 } else { n } {
+        let t = Instant::now();
+        let q = black_box(filter());
+        best = best.min(t.elapsed().as_secs_f64() * 1000.0);
+        if heap {
+            let stats = dhat::HeapStats::get();
+            println!(
+                "peak_bytes={} total_bytes={} total_blocks={} retained_output_bytes={}",
+                stats.max_bytes,
+                stats.total_bytes,
+                stats.total_blocks,
+                q.data.len() * size_of::<f32>()
+            );
+        }
+    }
+    drop(profiler);
+    println!("{mode} {w}x{h} sigma={sigma} step={} heap={heap} min_wall_ms={best:.3}", guided_fast_step(sigma));
+    Ok(())
+}

crates/pipeline/src/local.rs

@@ -178,10 +178,55 @@ fn guided_coeffs(p: &Plane, sigma: f32, eps: f32) -> (Plane, Plane) {
     par_join(|| gaussian(&a, sigma), || gaussian(&b, sigma))
 }
 
-/// `q = a·p + b` per pixel.
+/// `q = a·p + b` per pixel. Reduced coefficients are sampled with the same separable
+/// weights and accumulation order as a bilinear resize, without full-size coefficient planes.
 fn guided_apply(p: &Plane, a: &Plane, b: &Plane) -> Plane {
     let mut q = Plane::new(p.width, p.height);
     let w = p.width;
+    if (a.width, a.height) != (p.width, p.height) {
+        use lightcraft_raster::resample::{Filter, weights};
+        let wx = weights(a.width, w, Filter::Bilinear);
+        let wy = weights(a.height, p.height, Filter::Bilinear);
+        // Bilinear upsampling has at most four taps (including zero-weight boundary taps).
+        // Keep only those horizontally sampled rows per band, reusing them for neighbouring
+        // output rows. Scratch scales with width and active bands, not image area.
+        const BAND: usize = 32;
+        for_rows(&mut q.data, w.saturating_mul(BAND), |band, rows| {
+            let mut cache: [(usize, Vec<f32>, Vec<f32>); 4] = std::array::from_fn(|_| (usize::MAX, vec![0.0; w], vec![0.0; w]));
+            let mut bias = vec![0.0; w];
+            for (dy, row) in rows.chunks_mut(w).enumerate() {
+                let y = band * BAND + dy;
+                let (top, ys) = &wy[y];
+                bias.fill(0.0);
+                for (dy, yw) in ys.iter().enumerate() {
+                    let cy = top + dy;
+                    let (cached_y, aa, bb) = &mut cache[cy % 4];
+                    if *cached_y != cy {
+                        let (ar, br) = (a.row(cy), b.row(cy));
+                        for ((av, bv), (left, xs)) in aa.iter_mut().zip(bb.iter_mut()).zip(&wx) {
+                            *av = 0.0;
+                            *bv = 0.0;
+                            for (dx, xw) in xs.iter().enumerate() {
+                                *av += ar[left + dx] * xw;
+                                *bv += br[left + dx] * xw;
+                            }
+                        }
+                        *cached_y = cy;
+                    }
+                    for (av, ax) in row.iter_mut().zip(aa) {
+                        *av += *ax * yw;
+                    }
+                    for (bv, bx) in bias.iter_mut().zip(bb) {
+                        *bv += *bx * yw;
+                    }
+                }
+                for ((v, pv), bv) in row.iter_mut().zip(p.row(y)).zip(&bias) {
+                    *v = *v * pv + bv;
+                }
+            }
+        });
+        return q;
+    }
     for_rows(&mut q.data, w, |y, row| {
         for ((v, pv), (av, bv)) in row.iter_mut().zip(p.row(y)).zip(a.row(y).iter().zip(b.row(y))) {
             *v = av * pv + bv;
@@ -191,7 +236,7 @@ fn guided_apply(p: &Plane, a: &Plane, b: &Plane) -> Plane {
 }
 
 /// Fast guided filter (He & Sun 2015): coefficients computed on a `s×` subsampled plane and
-/// bilinearly upsampled; the output keeps full-resolution edges. Equivalent to [`guided`] for
+/// sampled bilinearly at full resolution; the output keeps full-resolution edges. Equivalent to [`guided`] for
 /// large windows at a fraction of the cost.
 pub fn guided_fast(p: &Plane, sigma: f32, eps: f32) -> Plane {
     use lightcraft_raster::resample::{Filter, resize};
@@ -202,8 +247,7 @@ pub fn guided_fast(p: &Plane, sigma: f32, eps: f32) -> Plane {
     let (lw, lh) = (p.width.div_ceil(s).max(1), p.height.div_ceil(s).max(1));
     let lo = resize(p, lw, lh, Filter::Box);
     let (ma, mb) = guided_coeffs(&lo, sigma / s as f32, eps);
-    let up = |c: &Plane| resize(c, p.width, p.height, Filter::Bilinear);
-    guided_apply(p, &up(&ma), &up(&mb))
+    guided_apply(p, &ma, &mb)
 }
 
 /// Separable median (`2r+1` px along each row, then each column; clamped edges). It keeps the
@@ -733,6 +777,60 @@ mod tests {
 mod fast_tests {
     use super::*;
 
+    /// The previous fast path, including its two materialized coefficient resizes.
+    fn resized_reference(p: &Plane, sigma: f32, eps: f32) -> Plane {
+        use lightcraft_raster::resample::{Filter, resize};
+        let s = guided_fast_step(sigma);
+        if s <= 1 {
+            return guided(p, sigma, eps);
+        }
+        let lo = resize(p, p.width.div_ceil(s).max(1), p.height.div_ceil(s).max(1), Filter::Box);
+        let (a, b) = guided_coeffs(&lo, sigma / s as f32, eps);
+        let (a, b) = (resize(&a, p.width, p.height, Filter::Bilinear), resize(&b, p.width, p.height, Filter::Bilinear));
+        Plane::from_fn(p.width, p.height, |x, y| a.get(x, y) * p.get(x, y) + b.get(x, y))
+    }
+
+    #[test]
+    fn sampled_coefficients_match_materialized_resizes() {
+        for (w, h) in [(0, 0), (1, 1), (1, 19), (23, 1), (2, 3), (17, 13), (64, 48), (161, 119)] {
+            let planes = [
+                Plane::filled(w, h, -3.25),
+                Plane::from_fn(w, h, |x, y| ((x * 7919 + y * 104729) % 97) as f32 / 7.0 - 8.0),
+                Plane::from_fn(w, h, |x, y| {
+                    if x == w / 2 && y == h / 2 {
+                        9.0
+                    } else if x < w / 2 {
+                        -5.0
+                    } else {
+                        2.0
+                    }
+                }),
+            ];
+            // Step 1 fallback, every small step, Draft's cap, and the maximum step (16).
+            for sigma in [0.0, 2.99, 5.99, 6.0, 9.0, 12.0, 24.0, 48.0, 200.0] {
+                for eps in [BASE_EPS, CLARITY_EPS, 0.0001] {
+                    for p in &planes {
+                        let expected = resized_reference(p, sigma, eps);
+                        let actual = guided_fast(p, sigma, eps);
+                        assert_eq!((actual.width, actual.height), (w, h));
+                        assert_eq!(actual.data, expected.data, "{w}x{h}, sigma {sigma}, eps {eps}");
+                    }
+                }
+            }
+        }
+    }
+
+    #[test]
+    fn sampled_tone_base_keeps_the_median_before_the_guided_filter() {
+        let p = Plane::from_fn(73, 51, |x, y| if x % 9 == 0 { 3.0 } else { -4.0 + ((x + y) % 7) as f32 * 0.03 });
+        for sigma in [6.0, 24.0, 120.0] {
+            for detail in [0, 1, 3, BASE_DETAIL_MAX] {
+                let median = median_hv(&p, detail);
+                assert_eq!(tone_base(&p, sigma, detail).data, resized_reference(&median, sigma, BASE_EPS).data, "{sigma}, {detail}");
+            }
+        }
+    }
+
     #[test]
     fn fast_guided_close_to_reference() {
         let p = Plane::from_fn(160, 120, |x, y| ((x as f32 * 0.13).sin() + (y as f32 * 0.07).cos()) * 2.0 + if x > 80 { 1.5 } else { 0.0 });

docs/gpu-pipeline.md

@@ -4,6 +4,10 @@ The CPU pipeline (`crates/pipeline`) is the reference ("oracle"). `lightcraft-gp
 evaluates the same stages with wgpu compute shaders (WGSL) on Metal / Vulkan / DX12. Everything is
 pure Rust (wgpu, naga); the drivers are the system's. No GL backend is compiled in.
 
+CPU fast guided-filter allocation measurements and reproduction commands are documented in
+[guided-filter-memory.md](guided-filter-memory.md). That optimization retains the full-size base
+and does not change the GPU's coefficient buffers.
+
 Warped geometry uses a CPU-computed coverage bit mask for the source boundary.
 The CPU reference evaluates that boundary in double precision; GPU float
 rounding otherwise can turn a blank edge pixel into a photo pixel. The mask

docs/guided-filter-memory.md

@@ -0,0 +1,53 @@
+# CPU fast guided-filter memory
+
+`local::guided_fast` computes blurred coefficients on the same reduced grid as before, then
+samples them while applying `a * p + b`. It keeps the raster resizer's weights, boundary
+normalization and horizontal-then-vertical accumulation order. Odd dimensions therefore use
+their actual reduced/full dimension ratios. Four coefficient row pairs and one bias row are cached
+per active 32-row band; adjacent output rows reuse their horizontal samples.
+
+The two full-resolution coefficient planes and their resize intermediates are gone. The returned
+full-resolution base remains: downstream CPU tone mapping and the spatial-plane cache need it.
+The step-1 filter, Box reduction, Gaussian coefficients, V2's median input, process selection and
+Draft sigma cap are unchanged. The GPU implementation is unchanged; this reduces CPU memory only.
+
+## Reproduce
+
+The dev-only example includes the previous materialized algorithm and uses the workspace's
+existing `dhat` allocator. No allocator or profiling dependency is added to product builds.
+
+```sh
+cargo build --release -p lightcraft-pipeline --example guided_memory
+RAYON_NUM_THREADS=1 target/release/examples/guided_memory old 6000 4000 120 15
+RAYON_NUM_THREADS=1 target/release/examples/guided_memory new 6000 4000 120 15
+RAYON_NUM_THREADS=1 target/release/examples/guided_memory old 6000 4000 120 1 --heap
+RAYON_NUM_THREADS=1 target/release/examples/guided_memory new 6000 4000 120 1 --heap
+```
+
+If `CARGO_TARGET_DIR` is set, use that directory instead of `target` for the executable.
+Arguments are mode, width, height, integer sigma and timing iterations. The example warms up
+before measurement; `--heap` counts one filter call with the source allocated beforehand.
+Time the runs without `--heap`: allocation profiling significantly changes timings. It reports
+the minimum wall time; on a shared machine repeat when other builds/renders have stopped.
+Use `LIGHTCRAFT_PROFILE=1` and `cargo xtask bench` for broader pipeline measurements.
+
+## Observed allocation and time
+
+Measured 2026-10-16 UTC, macOS aarch64, Rust 1.96.0, release build, one Rayon thread. Bytes are
+requested heap bytes counted by dhat, not process RSS. Timings are minimums of 15 unprofiled
+calls on a shared machine with concurrent builds; they are indicative, not a speed guarantee.
+
+| Dimensions / sigma / step | Old peak bytes | New peak bytes | Old total allocated bytes | New total allocated bytes | Old / new wall ms |
+|---|---:|---:|---:|---:|---:|
+| 6000 × 4000 / 120 / 16 | 289,125,000 | 97,780,856 | 313,083,488 | 135,643,856 | 42.135 / 45.131 |
+| 1601 × 1067 / 24 / 8 | 20,822,412 | 7,331,172 | 25,102,604 | 11,562,560 | 3.359 / 3.121 |
+
+The 24 MP result saves about 191 MB of peak allocation (66%). Its **96,000,000-byte output is
+included**; its **96,000,000-byte input is excluded**. Application RGB images, log luminance,
+V2 median input and its construction temporaries, other spatial planes and caches also remain
+outside this isolated measurement. More Rayon threads can hold more band scratch concurrently.
+This is not a bound on total pipeline memory.
+
+Regression tests compare the sampled result exactly against materialized bilinear resizes for
+empty/singleton/one-dimensional, odd/even and multiple-band images; flat fields, signed texture,
+edges and impulses; small/full and reduced filters; and median tone bases.