Box blur: count repeated edge pixels instead of visiting them

The starting window of each horizontal and vertical box pass visited every
position of the radius, including any number outside the image that repeat
the edge pixel, so a radius from a project could make one pass take
arbitrarily long. The pixels inside are summed once and the repeated edge
pixels are added as a count times the edge value.

crates/raster/src/blur.rs

@@ -27,13 +27,21 @@ fn box_h_into(src: &[Px], out: &mut [Px], w: usize, r: usize, repeat: bool) {
                 row[i as usize]
             }
         };
+        // The starting window: the pixels inside the row once each, and the repeated edge pixels
+        // counted rather than visited (a radius far past the row would otherwise loop over it).
         let mut acc = [0.0f32; 4];
-        for i in -(r as isize)..=(r as isize) {
-            let p = get(i);
+        let inside = r.min(w - 1);
+        for p in &row[..=inside] {
             for c in 0..4 {
                 acc[c] += p[c];
             }
         }
+        if repeat {
+            let (left, right) = (r as f32, (r - inside) as f32);
+            for c in 0..4 {
+                acc[c] += row[0][c] * left + row[w - 1][c] * right;
+            }
+        }
         for x in 0..w {
             o[x] = acc.map(|v| v * norm);
             let add = get(x as isize + r as isize + 1);
@@ -65,14 +73,23 @@ fn box_v_into(src: &[Px], out: &mut [Px], w: usize, h: usize, r: usize, repeat:
     out.par_chunks_mut(w * block).enumerate().for_each(|(bi, chunk)| {
         let ya = (bi * block) as isize;
         let mut acc = vec![[0.0f32; 4]; w];
-        for i in ya - r as isize..=ya + r as isize {
-            if let Some(rw) = row(i) {
-                for (a, p) in acc.iter_mut().zip(rw) {
-                    for c in 0..4 {
-                        a[c] += p[c];
-                    }
+        // The starting window: rows inside the image once each, the repeated edge rows counted.
+        let (lo, hi) = (ya - r as isize, ya + r as isize);
+        let add = |acc: &mut [Px], rw: &[Px], k: f32| {
+            for (a, p) in acc.iter_mut().zip(rw) {
+                for c in 0..4 {
+                    a[c] += p[c] * k;
                 }
             }
+        };
+        for i in lo.max(0)..=hi.min(h as isize - 1) {
+            add(&mut acc, &src[i as usize * w..(i as usize + 1) * w], 1.0);
+        }
+        if repeat {
+            let above = (-lo).max(0) as f32;
+            let below = (hi - (h as isize - 1)).max(0) as f32;
+            add(&mut acc, &src[0..w], above);
+            add(&mut acc, &src[(h - 1) * w..h * w], below);
         }
         for (yy, o) in chunk.chunks_mut(w).enumerate() {
             let y = ya + yy as isize;
@@ -264,6 +281,47 @@ pub fn radial_blur(img: &Image, center: (f64, f64), amount: f64, zoom: bool) ->
 mod tests {
     use super::*;
 
+    #[test]
+    fn box_windows_count_repeated_edges_and_take_any_radius() {
+        // Against a window summed pixel by pixel, for radii inside and far past the image.
+        let (w, h) = (7usize, 5usize);
+        let src: Vec<Px> = (0..w * h).map(|i| [i as f32, (i % 3) as f32, 1.0, 0.5]).collect();
+        let at = |x: isize, y: isize, repeat: bool| -> Px {
+            if !repeat && (x < 0 || y < 0 || x >= w as isize || y >= h as isize) {
+                return [0.0; 4];
+            }
+            src[y.clamp(0, h as isize - 1) as usize * w + x.clamp(0, w as isize - 1) as usize]
+        };
+        for repeat in [true, false] {
+            for r in [1usize, 3, 6, 9, 40] {
+                let (mut hx, mut vx) = (vec![[0.0; 4]; w * h], vec![[0.0; 4]; w * h]);
+                box_h_into(&src, &mut hx, w, r, repeat);
+                box_v_into(&src, &mut vx, w, h, r, repeat);
+                for y in 0..h as isize {
+                    for x in 0..w as isize {
+                        let (mut eh, mut ev) = ([0.0f32; 4], [0.0f32; 4]);
+                        for d in -(r as isize)..=r as isize {
+                            for c in 0..4 {
+                                eh[c] += at(x + d, y, repeat)[c];
+                                ev[c] += at(x, y + d, repeat)[c];
+                            }
+                        }
+                        let n = (2 * r + 1) as f32;
+                        let i = y as usize * w + x as usize;
+                        for c in 0..4 {
+                            assert!((hx[i][c] - eh[c] / n).abs() < 1e-3, "h r={r} repeat={repeat}");
+                            assert!((vx[i][c] - ev[c] / n).abs() < 1e-3, "v r={r} repeat={repeat}");
+                        }
+                    }
+                }
+            }
+        }
+        // A radius of 2^40 finishes at once.
+        let mut out = vec![[0.0; 4]; w * h];
+        box_h_into(&src, &mut out, w, 1 << 40, true);
+        box_v_into(&src, &mut out.clone(), w, h, 1 << 40, true);
+    }
+
     #[test]
     fn gaussian_preserves_mass() {
         let mut img = Image::new(64, 64);