shifathosn
/
photocraft
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compute exact octagonal mask shifts with bounded workspace (#2154) 

Co-authored-by: Brandon Thomas <echelon@gmail.com>
This commit is contained in:
Ayan 2026-10-10 15:46:06 +05:00 committed by GitHub
parent ae96037672
commit d99562f038
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
5 changed files with 681 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
crates/algo/src/matting.rs
View File

@ -217,36 +217,15 @@ pub fn gaussian_blur(src: &[f32], w: usize, h: usize, sigma: f32) -> Vec<f32> {
out
}
/// Grey dilation (`grow`) or erosion by about `r` pixels: alternating 3×3 square and cross steps
/// (an octagonal structuring element).
/// Grey dilation (`grow`) or erosion with the octagonal element of `r` alternating
/// 3×3 square/cross steps. Masks without NaNs or negative zero use O(w*h*log(r+1))
/// work; those values retain the original comparison order. Invalid sizes return an empty buffer.
pub fn morph(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
let mut cur = src.to_vec();
let mut nxt = cur.clone();
let pick = |a: f32, b: f32| if grow { a.max(b) } else { a.min(b) };
for k in 0..r {
let square = k % 2 == 0;
for y in 0..h {
for x in 0..w {
let mut v = cur[y * w + x];
for dy in -1i32..=1 {
for dx in -1i32..=1 {
if (dx != 0 && dy != 0 && !square) || (dx == 0 && dy == 0) {
continue;
}
let (nx, ny) = (x as i32 + dx, y as i32 + dy);
if nx >= 0 && ny >= 0 && (nx as usize) < w && (ny as usize) < h {
v = pick(v, cur[ny as usize * w + nx as usize]);
}
}
}
nxt[y * w + x] = v;
}
}
std::mem::swap(&mut cur, &mut nxt);
}
cur
morphology::morph(src, w, h, r, grow)
}
mod morphology;
/// Running max (or min) over a window of radius `r` (edge-clipped), van Herk / Gil–Werman:
/// three comparisons per sample whatever the radius.
fn running_extreme(src: &[f32], out: &mut [f32], r: usize, max: bool, g: &mut Vec<f32>, hb: &mut Vec<f32>) {
@ -276,7 +255,7 @@ fn max_square(src: &[f32], w: usize, h: usize, r: usize, max: bool) -> Vec<f32>
let (mut g, mut hb) = (Vec::new(), Vec::new());
let mut tmp = vec![0.0f32; w * h];
for y in 0..h {
running_extreme(&src[y * w..(y + 1) * w], &mut tmp[y * w..(y + 1) * w], r, max, &mut g, &mut hb);
running_extreme(&src[y * w..(y + 1) * w], &mut tmp[y * w..(y + 1) * w], r.min(w.saturating_sub(1)), max, &mut g, &mut hb);
}
let mut out = vec![0.0f32; w * h];
let (mut col, mut res) = (vec![0.0f32; h], vec![0.0f32; h]);
@ -284,7 +263,7 @@ fn max_square(src: &[f32], w: usize, h: usize, r: usize, max: bool) -> Vec<f32>
for y in 0..h {
col[y] = tmp[y * w + x];
}
running_extreme(&col, &mut res, r, max, &mut g, &mut hb);
running_extreme(&col, &mut res, r.min(h.saturating_sub(1)), max, &mut g, &mut hb);
for y in 0..h {
out[y * w + x] = res[y];
}
@ -314,6 +293,13 @@ pub fn edge_width(img: &RgbImage, r: usize) -> Vec<f32> {
/// Refines a soft mask buffer `m` (`w × h`) with `img` as guide (required when `radius > 0`).
pub fn refine_buffer(img: Option<&RgbImage>, m: &[f32], w: usize, h: usize, p: &RefineParams) -> Vec<f32> {
refine_buffer_with(img, m, w, h, p, morph)
}
fn refine_buffer_with<F>(img: Option<&RgbImage>, m: &[f32], w: usize, h: usize, p: &RefineParams, morphology: F) -> Vec<f32>
where
F: Fn(&[f32], usize, usize, usize, bool) -> Vec<f32> + Copy,
{
let mut a = m.to_vec();
if let (true, Some(img)) = (p.uses_image(), img) {
let bin: Vec<bool> = m.iter().map(|v| *v >= 0.5).collect();
@ -362,7 +348,7 @@ pub fn refine_buffer(img: Option<&RgbImage>, m: &[f32], w: usize, h: usize, p: &
}
let d = p.shift_px();
if d != 0 {
a = morph(&a, w, h, d.unsigned_abs() as usize, d > 0);
a = morphology(&a, w, h, d.unsigned_abs() as usize, d > 0);
}
a
}
@ -371,6 +357,20 @@ pub fn refine_buffer(img: Option<&RgbImage>, m: &[f32], w: usize, h: usize, p: &
/// its non-zero pixels; the result is clipped to `canvas`. `sampler` supplies the guide image
/// (only read when `radius > 0`). Returns `None` when nothing remains selected.
pub fn refine_mask(sampler: &dyn Sampler, mask: &(dyn Fn(Rect) -> Vec<f32> + Sync), content: Rect, canvas: Rect, p: &RefineParams) -> Option<Region> {
refine_mask_with(sampler, mask, content, canvas, p, morph)
}
fn refine_mask_with<F>(
sampler: &dyn Sampler,
mask: &(dyn Fn(Rect) -> Vec<f32> + Sync),
content: Rect,
canvas: Rect,
p: &RefineParams,
morphology: F,
) -> Option<Region>
where
F: Fn(&[f32], usize, usize, usize, bool) -> Vec<f32> + Copy + Sync + Send,
{
let content = content.intersect(&canvas);
if content.is_empty() {
return None;
@ -405,7 +405,7 @@ pub fn refine_mask(sampler: &dyn Sampler, mask: &(dyn Fn(Rect) -> Vec<f32> + Syn
return (*t, vec![q(first); cw * ch]);
}
let img = p.uses_image().then(|| sampler.rgb(r));
let out = refine_buffer(img.as_ref(), &m, rw, rh, p);
let out = refine_buffer_with(img.as_ref(), &m, rw, rh, p, morphology);
(*t, core(&out))
};
#[cfg(not(target_arch = "wasm32"))]

185
crates/algo/src/matting/morphology.rs Normal file
View File

@ -0,0 +1,185 @@
//! Exact flat octagonal morphology: a square followed by a Manhattan ball.
pub(super) fn morph(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
let Some(n) = w.checked_mul(h).filter(|&n| n <= src.len()) else { return Vec::new() };
if n == 0 || r == 0 {
return src.to_vec();
}
// The old max/min visitation order matters for NaN payloads and signed zero.
// Keep that path, and small radii where setup cannot pay off.
if r <= 2 || src[..n].iter().any(|v| v.is_nan() || (*v == 0.0 && v.is_sign_negative())) {
return repeated(src, w, h, r, grow);
}
let mut out = if n >= 1048576 && w >= 256 && h >= 128 && r <= 64 { tiled(&src[..n], w, h, r, grow) } else { dense(&src[..n], w, h, r, grow) };
out.extend_from_slice(&src[n..]);
out
}
/// Bound scratch by a 256x128 core plus the radius halo. The returned image is
/// the only full-size allocation. Constant input tiles need no filter scratch.
fn tiled(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
let mut out = vec![0.0; src.len()];
let run = |band: usize, rows: &mut [f32]| {
let y = band * 128;
let height = rows.len() / w;
let (y0, y1) = (y.saturating_sub(r), y.saturating_add(height).saturating_add(r).min(h));
for x in (0..w).step_by(256) {
let width = (w - x).min(256);
let (x0, x1) = (x.saturating_sub(r), x.saturating_add(width).saturating_add(r).min(w));
let first = src[y0 * w + x0];
if (y0..y1).all(|yy| src[yy * w + x0..yy * w + x1].iter().all(|&v| v == first)) {
for row in rows.chunks_mut(w) {
row[x..x + width].fill(first);
}
continue;
}
let tw = x1 - x0;
let mut input = Vec::with_capacity(tw * (y1 - y0));
for yy in y0..y1 {
input.extend_from_slice(&src[yy * w + x0..yy * w + x1]);
}
let filtered = dense(&input, tw, y1 - y0, r, grow);
for (yy, row) in rows.chunks_mut(w).enumerate() {
let base = (y - y0 + yy) * tw + x - x0;
row[x..x + width].copy_from_slice(&filtered[base..base + width]);
}
}
};
#[cfg(not(target_arch = "wasm32"))]
{
use rayon::prelude::*;
out.par_chunks_mut(w * h.min(128)).enumerate().for_each(|(band, rows)| run(band, rows));
}
#[cfg(target_arch = "wasm32")]
out.chunks_mut(w * h.min(128)).enumerate().for_each(|(band, rows)| run(band, rows));
out
}
fn dense(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
// ceil(r/2) unit squares plus floor(r/2) unit crosses. Saturate radii
// at the image extent: additional growth/erosion cannot reach another pixel.
let extent = w.max(h) - 1;
let square = (r / 2 + r % 2).min(extent);
let diamond = (r / 2).min(w.saturating_add(h).saturating_sub(2));
let mut cur = square_extreme(src, w, h, square, grow);
if diamond > 0 {
let mut next = vec![0.0; src.len()];
let mut built = 0usize;
for bit in (0..usize::BITS - diamond.leading_zeros()).rev() {
if built > 0 {
cross(&cur, &mut next, w, h, built, grow);
std::mem::swap(&mut cur, &mut next);
built *= 2;
}
if (diamond >> bit) & 1 != 0 {
cross(&cur, &mut next, w, h, 1, grow);
std::mem::swap(&mut cur, &mut next);
built += 1;
}
}
}
cur
}
fn square_extreme(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
#[cfg(not(target_arch = "wasm32"))]
if src.len() >= 262144 {
return parallel_square(src, w, h, r, grow);
}
super::max_square(src, w, h, r, grow)
}
/// Two full-size work planes, reused across blocked transposes. Filtering each
/// column contiguously avoids the full-photo stride of the serial square pass.
#[cfg(not(target_arch = "wasm32"))]
fn parallel_square(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
use rayon::prelude::*;
let mut rows = vec![0.0; src.len()];
rows.par_chunks_mut(w)
.enumerate()
.for_each_init(|| (Vec::new(), Vec::new()), |(g, hb), (y, out)| super::running_extreme(&src[y * w..(y + 1) * w], out, r.min(w - 1), grow, g, hb));
let mut cols = vec![0.0; src.len()];
// Independent bands of up to 32 columns, transposed in 32x32 cache blocks.
cols.par_chunks_mut(h * w.min(32)).enumerate().for_each(|(band, out)| {
let x0 = band * 32;
for y0 in (0..h).step_by(32) {
let y1 = y0.saturating_add(32).min(h);
for (x, col) in out.chunks_mut(h).enumerate() {
for (y, v) in col[y0..y1].iter_mut().enumerate() {
*v = rows[(y0 + y) * w + x0 + x];
}
}
}
});
cols.par_chunks_mut(h).for_each_init(
|| (Vec::new(), Vec::new(), vec![0.0; h]),
|(g, hb, out), col| {
super::running_extreme(col, out, r.min(h - 1), grow, g, hb);
col.copy_from_slice(out);
},
);
rows.par_chunks_mut(w * h.min(32)).enumerate().for_each(|(band, out)| {
let y0 = band * 32;
let height = out.len() / w;
for x0 in (0..w).step_by(32) {
for x in x0..x0.saturating_add(32).min(w) {
let col = &cols[x * h + y0..x * h + y0 + height];
for (y, &v) in col.iter().enumerate() {
out[y * w + x] = v;
}
}
}
});
rows
}
/// A radius-2d Manhattan ball is the union of five radius-d balls, centered
/// here and d pixels in each axial direction. Clamp the centers at the canvas
/// edge: every point in the clipped ball still belongs to one of those balls.
fn cross(src: &[f32], dst: &mut [f32], w: usize, h: usize, d: usize, grow: bool) {
let pick = |a: f32, b: f32| if grow { a.max(b) } else { a.min(b) };
let row = |y: usize, out: &mut [f32]| {
let a = &src[y * w..(y + 1) * w];
let up = y.saturating_sub(d) * w;
let down = y.saturating_add(d).min(h - 1) * w;
for (x, v) in out.iter_mut().enumerate() {
*v = pick(pick(pick(pick(a[x], a[x.saturating_sub(d)]), a[x.saturating_add(d).min(w - 1)]), src[up + x]), src[down + x]);
}
};
if dst.len() >= 262144 {
crate::photo_util::par_rows(dst, w, 1, row);
} else {
dst.chunks_mut(w).enumerate().for_each(|(y, out)| row(y, out));
}
}
fn repeated(src: &[f32], w: usize, h: usize, r: usize, grow: bool) -> Vec<f32> {
let mut cur = src.to_vec();
let mut nxt = cur.clone();
let pick = |a: f32, b: f32| if grow { a.max(b) } else { a.min(b) };
for k in 0..r {
let square = k % 2 == 0;
for y in 0..h {
for x in 0..w {
let mut v = cur[y * w + x];
for dy in -1i32..=1 {
for dx in -1i32..=1 {
if (dx != 0 && dy != 0 && !square) || (dx == 0 && dy == 0) {
continue;
}
let (nx, ny) = (x as i32 + dx, y as i32 + dy);
if nx >= 0 && ny >= 0 && (nx as usize) < w && (ny as usize) < h {
v = pick(v, cur[ny as usize * w + nx as usize]);
}
}
}
nxt[y * w + x] = v;
}
}
std::mem::swap(&mut cur, &mut nxt);
}
cur
}
#[cfg(test)]
mod tests;

317
crates/algo/src/matting/morphology/tests.rs Normal file
View File

@ -0,0 +1,317 @@
use super::*;
#[test]
fn every_small_binary_mask_matches_repeated_octagons() {
for (w, h) in [(1, 1), (1, 7), (7, 1), (2, 3), (3, 3), (4, 3)] {
for bits in 0..1usize << (w * h) {
let src: Vec<_> = (0..w * h).map(|i| ((bits >> i) & 1) as f32).collect();
for r in [0, 1, 2, 3, 4, 5, 7, 10, 13] {
for grow in [false, true] {
assert_eq!(morph(&src, w, h, r, grow), repeated(&src, w, h, r, grow), "{w}x{h} bits={bits} r={r} grow={grow}");
}
}
}
}
}
#[test]
fn gray_and_hdr_values_match_at_edges_and_every_radius() {
for (w, h) in [(2, 29), (29, 2), (37, 29), (1, 51), (51, 1)] {
for depth in [8, 16, 32] {
let src: Vec<_> = (0..w * h)
.map(|i| {
let v = ((i * 101 + i / w * 37) % 251) as f32 / 250.0;
match depth {
8 => (v * 255.0).round() / 255.0,
16 => (v * 65535.0).round() / 65535.0,
_ => v * 4.0 - 1.0,
}
})
.collect();
for r in 0..=65 {
for grow in [false, true] {
assert_eq!(morph(&src, w, h, r, grow), repeated(&src, w, h, r, grow), "{w}x{h} depth={depth} r={r} grow={grow}");
}
}
}
}
}
#[test]
fn impulses_match_the_closed_octagonal_footprint() {
let (w, h) = (97usize, 83usize);
for (sx, sy) in [(0, 0), (0, 41), (48, 41), (96, 82), (5, 79)] {
for r in 1..=65 {
for grow in [true, false] {
let mut src = vec![if grow { 0.0 } else { 1.0 }; w * h];
src[sy * w + sx] = if grow { 1.0 } else { 0.0 };
let out = morph(&src, w, h, r, grow);
let (a, b) = (r / 2 + r % 2, r / 2);
for y in 0..h {
for x in 0..w {
let affected = x.abs_diff(sx).saturating_sub(a) + y.abs_diff(sy).saturating_sub(a) <= b;
let expected = if affected == grow { 1.0 } else { 0.0 };
assert_eq!(out[y * w + x], expected, "seed {sx},{sy} r={r} grow={grow} at {x},{y}");
}
}
}
}
}
}
#[test]
fn exceptional_floats_keep_the_old_visitation_order() {
let (w, h) = (13, 9);
for special in [f32::from_bits(0x7fc00001), f32::NEG_INFINITY, f32::INFINITY, -0.0] {
let src: Vec<_> = (0..w * h).map(|i| if i % 7 == 0 { special } else { (i % 11) as f32 / 10.0 }).collect();
for r in [1, 3, 4, 17, 64] {
for grow in [false, true] {
let a = morph(&src, w, h, r, grow);
let b = repeated(&src, w, h, r, grow);
assert_eq!(a.iter().map(|v| v.to_bits()).collect::<Vec<_>>(), b.iter().map(|v| v.to_bits()).collect::<Vec<_>>());
}
}
}
let src = vec![f32::from_bits(0x7fc00002); w * h];
assert_eq!(morph(&src, w, h, 17, true).iter().map(|v| v.to_bits()).collect::<Vec<_>>(), src.iter().map(|v| v.to_bits()).collect::<Vec<_>>());
}
#[test]
fn invalid_sizes_and_large_radii_are_bounded() {
assert!(morph(&[], usize::MAX, 2, 3, true).is_empty());
assert!(morph(&[0.0], 2, 2, 3, true).is_empty());
assert_eq!(morph(&[], 0, 5, usize::MAX, true), Vec::<f32>::new());
assert_eq!(morph(&[0.0, 0.5, 1.0, 0.1], 2, 2, usize::MAX, true), vec![1.0; 4]);
assert_eq!(morph(&[0.0, 0.5, 1.0, 0.1, 42.0], 2, 2, 5, true), vec![1.0, 1.0, 1.0, 1.0, 42.0]);
}
#[test]
fn saturated_radius_does_not_pad_a_thin_axis_to_the_long_axis() {
let (w, h) = (2, 30001);
let src: Vec<_> = (0..w * h).map(|i| (i % 251) as f32).collect();
assert_eq!(morph(&src, w, h, usize::MAX, true), vec![250.0; w * h]);
assert_eq!(morph(&src, w, h, usize::MAX, false), vec![0.0; w * h]);
}
#[cfg(not(target_arch = "wasm32"))]
#[test]
fn worker_counts_are_bit_identical() {
let (w, h) = (769, 351);
let src: Vec<_> = (0..w * h).map(|i| ((i * 101 + i / w * 37) % 251) as f32 / 250.0).collect();
let mut reference = None;
for workers in [1, 2, 4] {
let pool = rayon::ThreadPoolBuilder::new().num_threads(workers).build().unwrap();
let result = pool.install(|| morph(&src, w, h, 33, false));
if let Some(a) = &reference { assert_eq!(a, &result) } else { reference = Some(result) }
}
}
#[cfg(not(target_arch = "wasm32"))]
#[test]
fn parallel_square_matches_serial_for_partial_blocks_and_thin_images() {
for (w, h) in [(769, 351), (513, 512), (3, 90001), (90001, 3)] {
let src: Vec<_> = (0..w * h).map(|i| ((i * 101 + i / w * 37) % 251) as f32 / 250.0).collect();
for r in [0, 1, 3, 16, 32] {
for grow in [true, false] {
assert_eq!(parallel_square(&src, w, h, r, grow), crate::matting::max_square(&src, w, h, r, grow), "{w}x{h} r={r} grow={grow}");
}
}
}
}
#[test]
fn shifted_tiles_match_a_whole_mask_at_canvas_edges() {
use crate::{
matting::{RefineParams, refine_buffer, refine_mask, region_reader},
segment::{ImageSampler, RgbImage},
selection::Region,
};
use photocraft_geom::Rect;
let (w, h) = (611usize, 317usize);
let img = RgbImage::from_fn(w, h, |_, _| [0.5; 3]);
let sampler = ImageSampler { img: &img, origin: (0, 0) };
let mask: Vec<_> = (0..w * h).map(|i| if (i % w < 63 && i / w < 191) || (i % w > 280 && i / w > 171) { 255 } else { 0 }).collect();
let rect = Rect::new(0, 0, w as i32, h as i32);
let region = Region { bbox: rect, mask };
let input: Vec<_> = region.mask.iter().map(|v| *v as f32 / 255.0).collect();
for shift in [-100.0, -51.0, 51.0, 100.0] {
// Constant guide keeps this check independent of image detail.
let p = RefineParams { radius: 32.0, shift_edge: shift, ..Default::default() };
let whole = refine_buffer(Some(&img), &input, w, h, &p);
let tiled = refine_mask(&sampler, &region_reader(&region), rect, rect, &p).unwrap();
for y in 0..h {
for x in 0..w {
assert!((tiled.at(x as i32, y as i32) - whole[y * w + x]).abs() <= 1.5 / 255.0, "shift={shift} at {x},{y}");
}
}
}
}
#[test]
fn bounded_tiles_match_dense_windows_on_partial_cores_and_fractional_values() {
let (w, h) = (1031usize, 1021usize);
for kind in 0..2 {
let src: Vec<_> = (0..w * h)
.map(|i| {
let (x, y) = (i % w, i / w);
if kind == 0 {
x as f32 / 1000.0 + y as f32 / 3000.0 - 0.3
} else if (x.abs_diff(255) + y.abs_diff(127)) < 35 || x == 768 && y == 896 {
0.75
} else {
0.0
}
})
.collect();
for r in [3, 16, 33, 64] {
for grow in [false, true] {
assert_eq!(morph(&src, w, h, r, grow), dense(&src, w, h, r, grow), "kind={kind} r={r} grow={grow}");
}
}
}
for value in [0.0, 0.5, 1.0, f32::INFINITY, f32::NEG_INFINITY] {
let mut src = vec![value; w * h];
src.push(42.0);
for grow in [true, false] {
assert_eq!(morph(&src, w, h, 64, grow), src);
}
}
}
#[cfg(not(target_arch = "wasm32"))]
#[test]
fn bounded_tiles_match_across_worker_counts() {
let (w, h) = (1031usize, 1021usize);
let src: Vec<_> = (0..w * h).map(|i| (i % w) as f32 / 1000.0 + (i / w) as f32 / 3000.0).collect();
let reference = dense(&src, w, h, 64, false);
for workers in [1, 2, 4] {
let pool = rayon::ThreadPoolBuilder::new().num_threads(workers).build().unwrap();
assert_eq!(pool.install(|| morph(&src, w, h, 64, false)), reference);
}
}
#[test]
#[ignore = "release comparison on tiles, 24 and 36 MP masks"]
#[allow(clippy::assertions_on_constants)]
fn octagonal_morphology_release_comparison() {
use std::{hint::black_box, time::Instant};
assert!(!cfg!(debug_assertions), "run with --release");
for (w, h) in [(384usize, 384usize), (6000, 4000), (7200, 5000)] {
let src: Vec<_> = (0..w * h)
.map(|i| {
let (x, y) = (i % w, i / w);
let dx = x as f64 - w as f64 * 0.51;
let dy = y as f64 - h as f64 * 0.48;
let d = (dx * dx + dy * dy).sqrt();
(0.5 + (w.min(h) as f64 * 0.27 - d) / 7.0).clamp(0.0, 1.0) as f32
})
.collect();
for r in [4, 16, 32, 64] {
for grow in [true, false] {
let reference = repeated(&src, w, h, r, grow);
assert_eq!(morph(&src, w, h, r, grow), reference);
let (mut old_ms, mut new_ms) = (Vec::new(), Vec::new());
for run in 0..3 {
let mut measure = |fast: bool| {
let start = Instant::now();
let out = if fast { morph(black_box(&src), w, h, r, grow) } else { repeated(black_box(&src), w, h, r, grow) };
let elapsed = start.elapsed().as_secs_f64() * 1000.0;
assert_eq!(out, reference);
if fast { new_ms.push(elapsed) } else { old_ms.push(elapsed) }
};
measure(run % 2 == 0);
measure(run % 2 != 0);
}
old_ms.sort_by(f64::total_cmp);
new_ms.sort_by(f64::total_cmp);
println!("OCTAGON {w}x{h} r={r} grow={grow} old_ms={old_ms:?} new_ms={new_ms:?} speedup={:.2}", old_ms[1] / new_ms[1]);
}
}
}
}
#[test]
#[ignore = "release comparison of the full tiled Select and Mask CPU path"]
#[allow(clippy::assertions_on_constants)]
fn tiled_refinement_release_comparison() {
use crate::{
matting::{RefineParams, refine_mask_with, region_reader},
segment::{ImageSampler, RgbImage, trim_region},
selection::Region,
};
use photocraft_geom::Rect;
use std::{hint::black_box, time::Instant};
assert!(!cfg!(debug_assertions), "run with --release");
for (w, h) in [(6000usize, 4000usize), (7200, 5000)] {
let image = RgbImage::from_fn(w, h, |x, y| {
let t = x as f32 / w as f32;
let u = y as f32 / h as f32;
[0.2 + 0.5 * t, 0.15 + 0.55 * u, 0.3 + 0.15 * ((x / 17 + y / 23) % 3) as f32]
});
let sampler = ImageSampler { img: &image, origin: (0, 0) };
let canvas = Rect::new(0, 0, w as i32, h as i32);
let mask: Vec<_> = (0..w * h)
.map(|i| {
let dx = (i % w) as f64 - w as f64 * 0.51;
let dy = (i / w) as f64 - h as f64 * 0.48;
let d = (dx * dx + dy * dy).sqrt();
((0.5 + (w.min(h) as f64 * 0.27 - d) / 7.0).clamp(0.0, 1.0) * 255.0).round() as u8
})
.collect();
let region = trim_region(Region { bbox: canvas, mask }).unwrap();
let reader = region_reader(&region);
let p = RefineParams { radius: 64.0, shift_edge: 100.0, ..Default::default() };
let old = refine_mask_with(&sampler, &reader, region.bbox, canvas, &p, repeated).unwrap();
let new = refine_mask_with(&sampler, &reader, region.bbox, canvas, &p, morph).unwrap();
assert_eq!(new.bbox, old.bbox);
assert_eq!(new.mask, old.mask);
let (mut old_ms, mut new_ms) = (Vec::new(), Vec::new());
for run in 0..3 {
let mut measure = |fast: bool| {
let start = Instant::now();
let out = if fast {
refine_mask_with(&sampler, black_box(&reader), region.bbox, canvas, &p, morph)
} else {
refine_mask_with(&sampler, black_box(&reader), region.bbox, canvas, &p, repeated)
}
.unwrap();
let elapsed = start.elapsed().as_secs_f64() * 1000.0;
assert_eq!(out.bbox, old.bbox);
assert_eq!(out.mask, old.mask);
if fast { new_ms.push(elapsed) } else { old_ms.push(elapsed) }
};
measure(run % 2 == 0);
measure(run % 2 != 0);
}
old_ms.sort_by(f64::total_cmp);
new_ms.sort_by(f64::total_cmp);
println!("REFINE {w}x{h} radius=64 shift=100 old_ms={old_ms:?} new_ms={new_ms:?} speedup={:.2}", old_ms[1] / new_ms[1]);
}
}
#[test]
#[ignore = "release comparison on a dense fractional 24 MP matte"]
#[allow(clippy::assertions_on_constants)]
fn dense_matte_release_comparison() {
use std::{hint::black_box, time::Instant};
assert!(!cfg!(debug_assertions), "run with --release");
let (w, h, r) = (6000usize, 4000usize, 64usize);
let src: Vec<_> = (0..w * h).map(|i| 0.6 * (i % w) as f32 / w as f32 + 0.4 * (i / w) as f32 / h as f32).collect();
let reference = repeated(&src, w, h, r, true);
assert_eq!(morph(&src, w, h, r, true), reference);
let (mut old_ms, mut new_ms) = (Vec::new(), Vec::new());
for run in 0..3 {
let mut measure = |fast: bool| {
let start = Instant::now();
let out = if fast { morph(black_box(&src), w, h, r, true) } else { repeated(black_box(&src), w, h, r, true) };
let elapsed = start.elapsed().as_secs_f64() * 1000.0;
assert_eq!(out, reference);
if fast { new_ms.push(elapsed) } else { old_ms.push(elapsed) }
};
measure(run % 2 == 0);
measure(run % 2 != 0);
}
old_ms.sort_by(f64::total_cmp);
new_ms.sort_by(f64::total_cmp);
println!("DENSE {w}x{h} r={r} grow=true old_ms={old_ms:?} new_ms={new_ms:?} speedup={:.2}", old_ms[1] / new_ms[1]);
}

135
docs/octagonal-mask-performance.md Normal file
View File

@ -0,0 +1,135 @@
# Exact octagonal mask shifts
Select and Mask applies Shift Edge as flat grayscale dilation or erosion. Previously,
`matting::morph` ran one full 3x3 pass for each radius step, alternating a square and a
cross. A 64 px shift needed 64 image passes. The filter supports fractional coverage.
## Decomposition
For radius `r`, let `a = ceil(r / 2)` and `b = floor(r / 2)`. The original footprint is
exactly a square of radius `a` followed by a Manhattan ball of radius `b`:
```text
max(abs(dx) - a, 0) + max(abs(dy) - a, 0) <= b
```
Equivalently, `max(abs(dx), abs(dy)) <= r` and `abs(dx) + abs(dy) <= r + a`.
Reordering the unit squares and crosses preserves the footprint on a rectangular
canvas: a path to any included pixel can stay between its endpoints, inside the canvas.
No outside pixel participates in the old clipped passes or in the new result.
The square uses the existing van Herk / Gil-Werman prefix/suffix window extrema, with
constant comparisons per sample. Native dense buffers of at least 262,144 pixels run
the horizontal windows in parallel, transpose in 32x32 cache blocks, filter contiguous columns, then transpose
back. The two full-size float work planes are reused. Small inputs and wasm keep the
serial square pass. Each axis caps the window radius at that axis's extent.
The Manhattan ball is built from the radius's binary digits. A radius-`2d` ball is the
union of five radius-`d` balls, centered at the output and at its four axial offsets of
`d` pixels. For any point in the larger ball, choose its larger absolute coordinate.
Moving the center `d` pixels along that axis leaves Manhattan distance at most `d`.
Clamping that center to the canvas retains this property and adds no out-of-footprint
pixels. An odd radius adds one unit cross after doubling.
For masks without NaNs or negative zero, this changes the work from **O(width * height * radius)** to
**O(width * height * log(radius + 1))**. Large native cross passes also run in parallel.
Small radii 0-2 retain the original path. NaNs and negative zero retain the original
comparison order, including NaN payloads and signed-zero behavior. Other values,
including infinities and HDR mask values, use exact max/min comparisons with no
quantization or floating-point sums.
## Measurements
Recorded 2026-10-10 on an AWS c7i.4xlarge, Intel Xeon Platinum 8488C, 16 vCPUs,
32 GiB RAM, Rust 1.99.0, standard workspace release profile. Native Rayon uses 16
workers. The 384x384 case is serial in both paths; it represents the size of a 256 px
core with a 64 px halo. Large masks use 256x128 output cores with a radius halo,
scheduled across the native workers. Constant halo windows fill their cores directly.
Input is a circular grayscale mask with a 7 px fractional edge. Input generation is
outside timing. Each case warms both paths, then measures three runs each in alternating
old/new order. The table gives the median. Every run compares every output pixel exactly
against the original repeated-pass implementation.
| Mask | Radius | Operation | Before p50 | After p50 | Speedup |
|---|---:|---|---:|---:|---:|
| 384 x 384 | 4 px | Dilate | 5.559 ms | 3.301 ms | 1.68x |
| 384 x 384 | 4 px | Erode | 4.765 ms | 2.941 ms | 1.62x |
| 384 x 384 | 16 px | Dilate | 21.254 ms | 4.666 ms | 4.55x |
| 384 x 384 | 16 px | Erode | 18.393 ms | 4.080 ms | 4.51x |
| 384 x 384 | 32 px | Dilate | 42.498 ms | 5.469 ms | 7.77x |
| 384 x 384 | 32 px | Erode | 37.368 ms | 4.874 ms | 7.67x |
| 384 x 384 | 64 px | Dilate | 84.287 ms | 6.601 ms | 12.77x |
| 384 x 384 | 64 px | Erode | 72.744 ms | 5.603 ms | 12.98x |
| 6000 x 4000 | 4 px | Dilate | 986.909 ms | 41.678 ms | 23.68x |
| 6000 x 4000 | 4 px | Erode | 847.584 ms | 40.829 ms | 20.76x |
| 6000 x 4000 | 16 px | Dilate | 3630.533 ms | 50.315 ms | 72.16x |
| 6000 x 4000 | 16 px | Erode | 3069.906 ms | 47.808 ms | 64.21x |
| 6000 x 4000 | 32 px | Dilate | 7153.449 ms | 60.465 ms | 118.31x |
| 6000 x 4000 | 32 px | Erode | 6024.692 ms | 57.223 ms | 105.28x |
| 6000 x 4000 | 64 px | Dilate | 14414.481 ms | 99.749 ms | 144.51x |
| 6000 x 4000 | 64 px | Erode | 11942.596 ms | 83.131 ms | 143.66x |
| 7200 x 5000 | 4 px | Dilate | 1483.109 ms | 59.551 ms | 24.91x |
| 7200 x 5000 | 4 px | Erode | 1263.256 ms | 57.847 ms | 21.84x |
| 7200 x 5000 | 16 px | Dilate | 5463.782 ms | 67.911 ms | 80.46x |
| 7200 x 5000 | 16 px | Erode | 4593.086 ms | 65.145 ms | 70.51x |
| 7200 x 5000 | 32 px | Dilate | 10766.328 ms | 81.414 ms | 132.24x |
| 7200 x 5000 | 32 px | Erode | 9048.115 ms | 78.240 ms | 115.65x |
| 7200 x 5000 | 64 px | Dilate | 21348.501 ms | 124.007 ms | 172.16x |
| 7200 x 5000 | 64 px | Erode | 17968.390 ms | 106.289 ms | 169.05x |
## Full CPU refinement
The second comparison runs the actual tiled `refine_mask` pipeline, swapping only the
morphology callback between the previous loop and the new implementation. It includes
guide sampling, guided filtering, boundary distances, quantization, tile assembly and
output trimming. Both paths use the same synthetic RGB guide and a trimmed circular
selection, radius 64 and Shift Edge +100%. Each result's bounding box and every mask
byte must match. Warmup and three alternating paired runs use the same method as above.
Document sizes are 24 and 36 MP; the selection is cropped by the normal workflow.
| Document | Before p50 | After p50 | Speedup |
|---|---:|---:|---:|
| 6000 x 4000 | 3019.176 ms | 761.687 ms | 3.96x |
| 7200 x 5000 | 3705.934 ms | 920.891 ms | 4.02x |
A third kernel comparison uses a dense fractional gradient, `0.6*x/w + 0.4*y/h`,
so no tile can take the constant-window shortcut:
| Document | Before p50 | After p50 | Speedup |
|---|---:|---:|---:|
| 6000 x 4000 | 14207.853 ms | 394.622 ms | 36.00x |
Kernel tables include filter allocations. The CPU refinement table includes the rest
of the refinement function. Document history and UI rendering are outside that timing.
## Bounded workspace
For image-shaped mask buffers without NaNs or negative zero, at least 1,048,576 pixels, width at least 256,
height at least 128 and radii 3-64, the returned f32 output is the only full-image allocation. Each active worker uses at most a 384x256 halo window, with three f32
planes and less than 6 KiB of line scratch. At 16 workers the derived allocation bound,
including output, is about **115 MB at 24 MP** or **163 MB at 36 MP**, versus the old
**192 MB / 288 MB** two-plane workspace. This is about **40% / 43% less**. These are
allocation bounds, not measured process RSS; the source and caller-retained copies are
separate. For N pixels and P workers, the bound is `4*N + P*(12*384*256 + 6144)`
bytes. Constant windows need no filter scratch. Smaller calls and larger radii use
the dense decomposition, with two f32 work planes plus line scratch.
## Correctness and reproduction
Coverage includes all binary masks on small rectangular canvases up to 4x3; every radius
0-65 on fractional 8-bit, 16-bit and HDR float values; impulse footprints at corners,
edges and interior positions; thin images; partial transpose blocks; worker-count
identity; infinities, NaN payloads and signed zero. Shifted tiled masks are compared
with whole-buffer refinement across canvas and tile edges. Checked dimensions reject
short input and multiplication overflow. Radii larger than the canvas saturate to
reachable extents on the decomposition path, including very thin canvases.
```sh
cargo test --release -p photocraft-algo octagonal_morphology_release_comparison -- --ignored --nocapture --test-threads=1
cargo test --release -p photocraft-algo tiled_refinement_release_comparison -- --ignored --nocapture --test-threads=1
cargo test --release -p photocraft-algo dense_matte_release_comparison -- --ignored --nocapture --test-threads=1
```
The benchmark's old path is also the production fallback for exceptional values and
small radii. It uses the previous loop without a changed footprint or tie order.

13
docs/roadmap.md
View File

@ -167,6 +167,19 @@ undo/redo. Unmapped Camera Raw fields/versions remain opaque; Photoshop
acceptance of generated exports and pixel parity are still unverified. Corpus floors above are
unchanged.
2026-10-10: Select and Mask's Shift Edge now computes the exact octagonal footprint
with a square window extreme plus a dyadically built Manhattan ball. On a c7i.4xlarge
(16 vCPUs), full tiled CPU refinement at radius 64 / Shift Edge +100% takes
3.019 s -> 0.762 s at 24 MP (3.96x),
3.706 s -> 0.921 s at 36 MP (4.02x).
Three alternating paired release runs include guide sampling, guided filtering,
boundary distances, quantization and assembly. Every output mask matches exactly.
The dense fractional 24 MP radius-64 dilation kernel improves 36.00x;
circular soft-mask kernels improve 20.76-172.16x at 24-36 MP, radii 4-64.
Bounded halo tiles reduce the large-mask workspace allocation bound from 192/288 MB
to about 115/163 MB, including output and 16 workers' scratch, excluding the source.
See [method and results](octagonal-mask-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.