Surface Blur: 4x to 20x faster on 16-bit and 32-bit layers (ranked window sums)

The 8-bit sliding histogram (#631) can't serve 16-bit or float layers: a 16-bit threshold spans
thousands of levels and floats have none. The weight `1 - |v - v0| / t` is linear on each side of
the centre value, so the weighted window sum only needs the count, Σv and Σv² of the samples within
the threshold below and above it. A Fenwick tree over the window's distinct values holds those and
slides along the row, so a pixel costs O(radius · log n) instead of O(radius²), for any depth,
with the direct sum's own weight test deciding which samples count. Used from radius 10, where it
overtakes the direct sum; non-finite samples keep the direct sum.

crates/algo/src/blur.rs

@@ -483,9 +483,140 @@ pub(crate) fn surface(src: &Image, out: Rect, ctx: &Ctx, radius: f32, threshold:
     {
         return res;
     }
+    // 16-bit and float (and 8-bit where the histogram declined): the ranked sums overtake the
+    // direct sum near radius 8 on a 256 px tile and are about 10x faster at radius 100.
+    if r >= RANKED_MIN_RADIUS
+        && let Some(res) = surface_ranked(src, out, ctx, r, t)
+    {
+        return res;
+    }
     surface_direct(src, out, ctx, r, t)
 }
 
+/// Radius from which [`surface_ranked`] beats [`surface_direct`] (measured by `surface_timing`).
+const RANKED_MIN_RADIUS: i32 = 10;
+
+/// Fenwick tree over value ranks, each node holding (count, Σv, Σv²) of its range.
+struct RankSums(Vec<[f64; 3]>);
+
+impl RankSums {
+    fn add(&mut self, rank: usize, sign: f64, v: f64) {
+        let mut i = rank + 1;
+        while let Some(node) = self.0.get_mut(i) {
+            node[0] += sign;
+            node[1] += sign * v;
+            node[2] += sign * v * v;
+            i += i & i.wrapping_neg();
+        }
+    }
+
+    /// Sums over the ranks `0..end`.
+    fn prefix(&self, end: usize) -> [f64; 3] {
+        let mut acc = [0.0; 3];
+        let mut i = end.min(self.0.len().saturating_sub(1));
+        while i > 0 {
+            if let Some(node) = self.0.get(i) {
+                acc[0] += node[0];
+                acc[1] += node[1];
+                acc[2] += node[2];
+            }
+            i -= i & i.wrapping_neg();
+        }
+        acc
+    }
+
+    /// Sums over the ranks `start..end`.
+    fn range(&self, start: usize, end: usize) -> [f64; 3] {
+        let (a, b) = (self.prefix(start), self.prefix(end));
+        [b[0] - a[0], b[1] - a[1], b[2] - a[2]]
+    }
+}
+
+/// [`surface`] at any bit depth. The weight `1 - |v - v0| / t` is linear in `v` on each side of
+/// the centre value, so over the samples within the threshold below it `Σw = ((t - v0)·n + Σv) / t`
+/// and `Σv·w = ((t - v0)·Σv + Σv²) / t`, and symmetrically above it. A Fenwick tree over the
+/// window's distinct values holds the count, Σv and Σv² and slides along the row (one column out
+/// and one in per step), so a pixel costs O(radius · log n) instead of O(radius^2). `None` when a
+/// sample isn't finite; the direct sum handles those.
+fn surface_ranked(src: &Image, out: Rect, ctx: &Ctx, r: i32, t: f32) -> Option<Vec<f32>> {
+    let n = src.ch;
+    let cc = if ctx.alpha { n - 1 } else { n };
+    let win = out.inflate(r);
+    let (ww, wh) = (win.width() as usize, win.height() as usize);
+    // Every colour sample of the window (outside the source reads as 0, as `Image::get` does),
+    // and whether its pixel counts (the direct sum skips transparent ones).
+    let mut vals = Vec::with_capacity(ww * wh * cc);
+    let mut opaque = Vec::with_capacity(ww * wh);
+    for y in win.y0..win.y1 {
+        for x in win.x0..win.x1 {
+            opaque.push(!(ctx.alpha && src.get(x, y, n - 1) <= 0.0));
+            for c in 0..cc {
+                let v = src.get(x, y, c);
+                if !v.is_finite() {
+                    return None;
+                }
+                // `+ 0.0` turns -0.0 into 0.0, so equal values share a rank.
+                vals.push(v + 0.0);
+            }
+        }
+    }
+    let mut levels = vals.clone();
+    levels.sort_unstable_by(f32::total_cmp);
+    levels.dedup();
+    let rank: Vec<usize> = vals.iter().map(|v| levels.partition_point(|l| l < v)).collect();
+    let mut tree = RankSums(vec![[0.0; 3]; levels.len() + 1]);
+    let (ow, oh) = (out.width() as usize, out.height() as usize);
+    if ow == 0 || oh == 0 {
+        return Some(Vec::new());
+    }
+    let side = (2 * r + 1) as usize;
+    let t64 = f64::from(t);
+    let mut res = vec![0.0f32; ow * oh * n];
+    for oy in 0..oh {
+        for c in 0..cc {
+            let column = |tree: &mut RankSums, wx: usize, sign: f64| {
+                for wy in oy..oy + side {
+                    let p = wy * ww + wx;
+                    if opaque.get(p) == Some(&true)
+                        && let (Some(&k), Some(&v)) = (rank.get(p * cc + c), vals.get(p * cc + c))
+                    {
+                        tree.add(k, sign, f64::from(v));
+                    }
+                }
+            };
+            for wx in 0..side {
+                column(&mut tree, wx, 1.0);
+            }
+            for ox in 0..ow {
+                let v0 = src.px(out.x0 + ox as i32, out.y0 + oy as i32)[c];
+                // Ranks strictly within the threshold below v0 (and v0 itself), then above it.
+                let lo = levels.partition_point(|&l| 1.0 - (l - v0).abs() / t <= 0.0 && l < v0);
+                let mid = levels.partition_point(|&l| l <= v0);
+                let hi = levels.partition_point(|&l| l <= v0 || 1.0 - (l - v0).abs() / t > 0.0);
+                let (below, above) = (tree.range(lo, mid), tree.range(mid, hi));
+                let v0d = f64::from(v0);
+                let wsum = ((t64 - v0d) * below[0] + below[1] + (t64 + v0d) * above[0] - above[1]) / t64;
+                let acc = ((t64 - v0d) * below[1] + below[2] + (t64 + v0d) * above[1] - above[2]) / t64;
+                res[(oy * ow + ox) * n + c] = if below[0] + above[0] > 0.0 && wsum > 0.0 { (acc / wsum) as f32 } else { v0 };
+                if ox + 1 < ow {
+                    column(&mut tree, ox, -1.0);
+                    column(&mut tree, ox + side, 1.0);
+                }
+            }
+            // Empty the tree for the next row and channel.
+            for wx in ow - 1..ow - 1 + side {
+                column(&mut tree, wx, -1.0);
+            }
+        }
+        if ctx.alpha {
+            for ox in 0..ow {
+                res[(oy * ow + ox) * n + n - 1] = src.px(out.x0 + ox as i32, out.y0 + oy as i32)[n - 1];
+            }
+        }
+    }
+    Some(res)
+}
+
 /// [`surface`] by summing the whole (2r + 1)^2 window of every pixel.
 fn surface_direct(src: &Image, out: Rect, ctx: &Ctx, r: i32, t: f32) -> Vec<f32> {
     let n = src.ch;
@@ -702,4 +833,118 @@ mod tests {
         img.data[i] = 0.5 / 255.0;
         assert!(surface_8bit(&img, out, &ctx, 4, 0.1, 26).is_none());
     }
+
+    /// The direct sum with f64 accumulators: the reference both fast paths are held to.
+    fn surface_reference(src: &Image, out: Rect, ctx: &Ctx, r: i32, t: f32) -> Vec<f32> {
+        let n = src.ch;
+        let cc = if ctx.alpha { n - 1 } else { n };
+        let mut res = Vec::new();
+        for y in out.y0..out.y1 {
+            for x in out.x0..out.x1 {
+                let p = src.px(x, y);
+                for &v0 in p.iter().take(cc) {
+                    let c = res.len() % n;
+                    let (mut acc, mut wsum) = (0.0f64, 0.0f64);
+                    for yy in y - r..=y + r {
+                        for xx in x - r..=x + r {
+                            if ctx.alpha && src.get(xx, yy, n - 1) <= 0.0 {
+                                continue;
+                            }
+                            let v = src.get(xx, yy, c);
+                            let w = (1.0 - (v - v0).abs() / t).max(0.0);
+                            acc += f64::from(v) * f64::from(w);
+                            wsum += f64::from(w);
+                        }
+                    }
+                    res.push(if wsum > 0.0 { (acc / wsum) as f32 } else { v0 });
+                }
+                if ctx.alpha {
+                    res.push(p[n - 1]);
+                }
+            }
+        }
+        res
+    }
+
+    #[test]
+    fn surface_blur_ranked_sums_match_the_weighted_window_at_any_depth() {
+        // 16-bit levels and arbitrary floats (including negatives and values above 1), with some
+        // fully transparent pixels, which the sum skips.
+        for (quant, channels, alpha) in [(Some(65535.0f32), 4, true), (None, 4, true), (None, 1, false)] {
+            let src_rect = Rect::new(-20, -20, 48, 40);
+            let mut img = Image::new(src_rect, channels);
+            let mut s = 0x1b87_3593u32;
+            for px in img.data.chunks_mut(channels) {
+                for v in px.iter_mut() {
+                    s ^= s << 13;
+                    s ^= s >> 17;
+                    s ^= s << 5;
+                    let raw = (s % 10_000) as f32 / 10_000.0;
+                    *v = quant.map_or(raw * 1.5 - 0.25, |q| (raw * q).round() / q);
+                }
+                if alpha && s.is_multiple_of(7) {
+                    px[channels - 1] = 0.0;
+                }
+            }
+            let ctx = Ctx { bounds: src_rect, mode: crate::ColorMode::Rgb, alpha };
+            // The output reaches the source edge, where samples read as 0.
+            let out = Rect::new(-6, 0, 30, 22);
+            for (r, threshold) in [(1, 2.0), (4, 15.0), (12, 60.0), (16, 255.0)] {
+                let t = (threshold / 255.0) * 2.5;
+                let ranked = surface_ranked(&img, out, &ctx, r, t).expect("finite input");
+                let reference = surface_reference(&img, out, &ctx, r, t);
+                assert_eq!(ranked.len(), reference.len());
+                let err = reference.iter().zip(&ranked).map(|(a, b)| (a - b).abs()).fold(0.0f32, f32::max);
+                assert!(err < 1e-6, "{quant:?} r {r} threshold {threshold}: max error {err}");
+            }
+        }
+    }
+
+    #[test]
+    fn surface_blur_ranked_declines_non_finite_samples_and_handles_empty_output() {
+        let src_rect = Rect::new(0, 0, 8, 8);
+        let mut img = Image::new(src_rect, 1);
+        img.data.fill(0.5);
+        let ctx = Ctx { bounds: src_rect, mode: crate::ColorMode::Grayscale, alpha: false };
+        assert_eq!(surface_ranked(&img, Rect::new(2, 2, 2, 6), &ctx, 2, 0.1), Some(Vec::new()));
+        img.data[9] = f32::NAN;
+        assert!(surface_ranked(&img, Rect::new(0, 0, 8, 8), &ctx, 2, 0.1).is_none());
+        // The public entry still answers, through the direct sum.
+        assert_eq!(surface(&img, Rect::new(0, 0, 8, 8), &ctx, 12.0, 15.0).len(), 64);
+    }
+
+    /// Release-mode timings of the direct sum and the ranked path on one 256 px tile of
+    /// 16-bit-like and float noise (`cargo test -p photocraft-algo --release -- --ignored --nocapture surface_timing`).
+    #[test]
+    #[ignore]
+    fn surface_timing() {
+        for (label, quant) in [("16-bit", Some(65535.0f32)), ("float", None)] {
+            for r in [2, 4, 8, 16, 32, 64, 100] {
+                let out = Rect::new(0, 0, 256, 256);
+                let src_rect = out.inflate(r);
+                let mut img = Image::new(src_rect, 4);
+                let mut s = 0x2545_f491u32;
+                let len = img.data.len();
+                for (i, v) in img.data.iter_mut().enumerate() {
+                    s ^= s << 13;
+                    s ^= s >> 17;
+                    s ^= s << 5;
+                    // Smooth gradient plus noise, so the threshold keeps a realistic share.
+                    let x = (i / 4) as f32 / len as f32 * 4.0;
+                    let raw = (0.5 + 0.4 * (x * 6.0).sin() + (s % 1000) as f32 / 20000.0).clamp(0.0, 1.0);
+                    *v = quant.map_or(raw, |q| (raw * q).round() / q);
+                }
+                let ctx = Ctx { bounds: src_rect, mode: crate::ColorMode::Rgb, alpha: true };
+                let t = (15.0 / 255.0) * 2.5;
+                let a = std::time::Instant::now();
+                let direct = surface_direct(&img, out, &ctx, r, t);
+                let td = a.elapsed();
+                let a = std::time::Instant::now();
+                let ranked = surface_ranked(&img, out, &ctx, r, t).expect("finite");
+                let tr = a.elapsed();
+                let err = direct.iter().zip(&ranked).map(|(a, b)| (a - b).abs()).fold(0.0f32, f32::max);
+                println!("{label} r {r:3}: direct {td:>10.2?}  ranked {tr:>10.2?}  max err {err:.2e}");
+            }
+        }
+    }
 }