From faa039b4c36dd7bd9818bf17dc32d3d4d2fb144d Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Wed, 15 Jul 2026 20:44:23 +0900 Subject: [PATCH] perf(vision): accelerate exact distance transform --- Cargo.toml | 1 + README.md | 30 +- bench/opencv_vision_comparison/README.md | 3 +- bench/opencv_vision_comparison/performance.py | 32 + crates/spatialrust-py/Cargo.toml | 4 + crates/spatialrust-py/spatialrust.pyi | 14 +- crates/spatialrust-py/src/lib.rs | 106 +++- crates/spatialrust-py/tests/test_bindings.py | 7 + crates/spatialrust-vision/Cargo.toml | 3 +- crates/spatialrust-vision/benches/dense.rs | 17 +- crates/spatialrust-vision/src/dense.rs | 581 ++++++++++++++++-- docs/ROADMAP.md | 2 + docs/site/algorithms.html | 2 +- notes/2026-07-15_exact_edt_acceleration.md | 39 ++ 14 files changed, 765 insertions(+), 76 deletions(-) create mode 100644 notes/2026-07-15_exact_edt_acceleration.md diff --git a/Cargo.toml b/Cargo.toml index b61a7ab..78481b6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -91,6 +91,7 @@ thiserror = "2" wgpu = "24" criterion = { version = "0.5", features = ["html_reports"] } proptest = "1" +rayon = "1.12" image = { version = "0.24.9", default-features = false } kamadak-exif = "0.6.1" exr = { version = "=1.72.0", default-features = false } diff --git a/README.md b/README.md index 53f0a89..c5e4db3 100644 --- a/README.md +++ b/README.md @@ -150,17 +150,18 @@ ratio; these are machine-specific measurements, not universal guarantees. | Workload | VGA | 1080p | 4K | | --- | ---: | ---: | ---: | -| AI CHW preprocess, allocate | **SpatialRust 4.52×** | **SpatialRust 7.19×** | **SpatialRust 10.51×** | -| AI CHW preprocess, reuse vs OpenCV allocate | **SpatialRust 8.54×** | **SpatialRust 11.46×** | **SpatialRust 17.13×** | -| Bilinear resize, allocate | OpenCV 28.3× | OpenCV 61.0× | OpenCV 62.8× | -| Bilinear resize, reuse | OpenCV 28.6× | OpenCV 110.8× | OpenCV 124.0× | -| RGB to gray, allocate | OpenCV 11.9× | OpenCV 8.8× | OpenCV 12.8× | -| RGB to gray, reuse | OpenCV 6.4× | OpenCV 12.5× | OpenCV 8.6× | -| Gaussian blur 5×5 | OpenCV 125.1× | OpenCV 121.6× | OpenCV 177.8× | -| Sobel X 3×3 | OpenCV 36.6× | OpenCV 35.6× | OpenCV 38.6× | -| Morphology open 5×5 | OpenCV 909.0× | OpenCV 1,027.2× | OpenCV 1,082.8× | -| Canny | OpenCV 13.0× | OpenCV 14.5× | OpenCV 15.3× | -| Exact Euclidean distance transform | OpenCV 10.60× | OpenCV 12.63× | OpenCV 12.35× | +| AI CHW preprocess, allocate | **SpatialRust 4.48×** | **SpatialRust 9.27×** | **SpatialRust 9.14×** | +| AI CHW preprocess, reuse vs OpenCV allocate | **SpatialRust 8.16×** | **SpatialRust 14.56×** | **SpatialRust 15.78×** | +| Bilinear resize, allocate | OpenCV 26.46× | OpenCV 64.02× | OpenCV 93.61× | +| Bilinear resize, reuse | OpenCV 27.36× | OpenCV 145.81× | OpenCV 113.87× | +| RGB to gray, allocate | OpenCV 11.97× | OpenCV 5.98× | OpenCV 12.98× | +| RGB to gray, reuse | OpenCV 6.01× | OpenCV 13.81× | OpenCV 2.09× | +| Gaussian blur 5×5 | OpenCV 139.02× | OpenCV 107.91× | OpenCV 167.28× | +| Sobel X 3×3 | OpenCV 14.38× | OpenCV 20.31× | OpenCV 23.30× | +| Morphology open 5×5 | OpenCV 805.78× | OpenCV 578.52× | OpenCV 774.40× | +| Canny | OpenCV 10.66× | OpenCV 12.54× | OpenCV 12.65× | +| Exact Euclidean distance transform, allocate | OpenCV 2.37× | OpenCV 1.99× | OpenCV 1.63× | +| Exact Euclidean distance transform, reuse | OpenCV 1.61× | OpenCV 1.20× | OpenCV 1.07× | The current CPU result is deliberately mixed: SpatialRust's fused typed CHW path wins, while OpenCV's tuned general-purpose image kernels lead the present @@ -169,6 +170,13 @@ samples are produced by the [performance harness](bench/opencv_vision_comparison the dated [Epic 111 receipt](notes/2026-07-15_epic111_opencv_comparison_v2.md) records the exact environment and methodology. +The EDT fast path is exact on the canonical masks and reduced the native 4K +allocation benchmark from 451.63 ms to about 75 ms. With caller-owned output +and [`DistanceTransformWorkspace`](https://rsasaki0109.github.io/SpatialRust/spatialrust_vision/struct.DistanceTransformWorkspace.html), +the native Criterion median is about 43 ms. The Python API comparison above +still gives OpenCV a narrow 1.07× 4K reuse lead; see the +[acceleration receipt](notes/2026-07-15_exact_edt_acceleration.md). + #### Vision accuracy The same deterministic RGB inputs passed all VGA, 1080p, and 4K gates: diff --git a/bench/opencv_vision_comparison/README.md b/bench/opencv_vision_comparison/README.md index e7f1e77..96d5b96 100644 --- a/bench/opencv_vision_comparison/README.md +++ b/bench/opencv_vision_comparison/README.md @@ -24,7 +24,8 @@ runtime dependency. The shared report contract and workload registry are in [`../opencv_comparison`](../opencv_comparison/README.md). The performance suite measures allocate/reuse bilinear resize, RGB-to-gray, -AI CHW preprocessing, Gaussian blur, Sobel X, morphology open, and Canny at +AI CHW preprocessing, Gaussian blur, Sobel X, morphology open, Canny, and exact +Euclidean distance transform (allocate and caller-owned-output/workspace reuse) at VGA, 1080p, and 4K. OpenCV and SpatialRust calls are seeded and interleaved; short calls are batched to reduce timer noise. Reports preserve raw samples, mean/median/p95, standard deviation, coefficient of variation, median absolute diff --git a/bench/opencv_vision_comparison/performance.py b/bench/opencv_vision_comparison/performance.py index cb146c1..0f058c0 100644 --- a/bench/opencv_vision_comparison/performance.py +++ b/bench/opencv_vision_comparison/performance.py @@ -338,6 +338,33 @@ def main() -> None: seed=113, min_sample_time_ms=MIN_SAMPLE_TIME_MS, ) + distance_cv_out = np.empty((height, width), dtype=np.float32) + distance_sr_out = np.empty((height, width), dtype=np.float32) + distance_workspace = sr.DistanceTransformWorkspace() + cv2.distanceTransform( + distance_mask, + cv2.DIST_L2, + cv2.DIST_MASK_PRECISE, + dst=distance_cv_out, + ) + sr.distance_transform_edt( + distance_mask, out=distance_sr_out, workspace=distance_workspace + ) + _, _, cv_distance_reuse, sr_distance_reuse = timed_pair( + lambda: cv2.distanceTransform( + distance_mask, + cv2.DIST_L2, + cv2.DIST_MASK_PRECISE, + dst=distance_cv_out, + ), + lambda: sr.distance_transform_edt( + distance_mask, out=distance_sr_out, workspace=distance_workspace + ), + warmup=args.warmup, + repeats=repeats, + seed=114, + min_sample_time_ms=MIN_SAMPLE_TIME_MS, + ) rows = ( ("resize_bilinear", "opencv", "allocate", cv_resize_alloc), @@ -361,6 +388,8 @@ def main() -> None: ("canny", "spatialrust", "allocate", sr_canny), ("distance_transform_edt", "opencv", "allocate", cv_distance), ("distance_transform_edt", "spatialrust", "allocate", sr_distance), + ("distance_transform_edt", "opencv", "reuse", cv_distance_reuse), + ("distance_transform_edt", "spatialrust", "reuse", sr_distance_reuse), ) measurements.extend( measurement(workload, implementation, mode, width, height, timing) @@ -380,6 +409,9 @@ def main() -> None: "morphology_open": speed_comparison(cv_morphology, sr_morphology), "canny": speed_comparison(cv_canny, sr_canny), "distance_transform_edt": speed_comparison(cv_distance, sr_distance), + "distance_transform_edt_reuse": speed_comparison( + cv_distance_reuse, sr_distance_reuse + ), } environment_receipt = environment( diff --git a/crates/spatialrust-py/Cargo.toml b/crates/spatialrust-py/Cargo.toml index 674ac56..4a7c0ae 100644 --- a/crates/spatialrust-py/Cargo.toml +++ b/crates/spatialrust-py/Cargo.toml @@ -53,3 +53,7 @@ spatialrust = { path = "../spatialrust", features = [ # Keep this crate out of the main Rust workspace so `cargo test --workspace` # in CI does not need a Python toolchain. [workspace] + +[profile.release] +lto = "thin" +codegen-units = 1 diff --git a/crates/spatialrust-py/spatialrust.pyi b/crates/spatialrust-py/spatialrust.pyi index 2db757d..a154bf0 100644 --- a/crates/spatialrust-py/spatialrust.pyi +++ b/crates/spatialrust-py/spatialrust.pyi @@ -285,8 +285,20 @@ def connected_components_image( mask: _U8Array, connectivity: int = ..., ) -> tuple[_U32Array, list[tuple[int, int, tuple[float, float, float, float]]]]: ... + +@final +class DistanceTransformWorkspace: + """Reusable host scratch storage for exact unit-spacing distance transforms.""" + + def __init__(self) -> None: ... + @property + def capacity(self) -> int: ... + def distance_transform_edt( - mask: _U8Array, spacing: tuple[float, float] = ... + mask: _U8Array, + spacing: tuple[float, float] = ..., + out: Optional[_F32Array] = ..., + workspace: Optional[DistanceTransformWorkspace] = ..., ) -> _F32Array: ... def find_mask_contours( mask: _U8Array, epsilon: float = ... diff --git a/crates/spatialrust-py/src/lib.rs b/crates/spatialrust-py/src/lib.rs index a278392..80f2228 100644 --- a/crates/spatialrust-py/src/lib.rs +++ b/crates/spatialrust-py/src/lib.rs @@ -75,6 +75,7 @@ use spatialrust::vision::{ connected_components as label_components, decode_rle as decode_mask_runs, detect_and_describe_orb as detect_and_describe_orb_op, detect_fast as detect_fast_op, detect_harris as detect_harris_op, detect_shi_tomasi as detect_shi_tomasi_op, + distance_transform_edt_u8_into as distance_transform_edt_u8_into_op, distance_transform_edt_with_spacing as distance_transform_edt_op, encode_rle as encode_mask_runs, equalize_histogram as equalize_histogram_op, estimate_homography_ransac as estimate_homography_ransac_op, @@ -92,12 +93,12 @@ use spatialrust::vision::{ solve_pnp as solve_pnp_op, stereo_block_match as stereo_block_match_op, stitch_panorama_pair as stitch_panorama_pair_op, threshold as threshold_op, AbsolutePose, AdaptiveThresholdMethod, BinaryMask, BorderMode, BoundingBox2, CameraMatrix3, CannyOptions, - ConfidenceMap, Connectivity, CornerSelectionOptions, DescriptorBuffer, FastOptions, - HarrisOptions, Interpolation, Kernel2D, Keypoint2, MaskRle, MatchOptions, MorphologyOperation, - MorphologyShape, ObjectImageCorrespondence, OrbOptions, OrbScoreType, PanoramaOptions, - PerspectiveTransform, PointCorrespondence2, PointMap, RgbdOdometryOptions, RleOrder, - RobustEstimationOptions, ShiTomasiOptions, SoftNmsMethod, StereoBmOptions, StructuringElement, - ThresholdType, + ConfidenceMap, Connectivity, CornerSelectionOptions, DescriptorBuffer, + DistanceTransformWorkspace, FastOptions, HarrisOptions, Interpolation, Kernel2D, Keypoint2, + MaskRle, MatchOptions, MorphologyOperation, MorphologyShape, ObjectImageCorrespondence, + OrbOptions, OrbScoreType, PanoramaOptions, PerspectiveTransform, PointCorrespondence2, + PointMap, RgbdOdometryOptions, RleOrder, RobustEstimationOptions, ShiTomasiOptions, + SoftNmsMethod, StereoBmOptions, StructuringElement, ThresholdType, }; use spatialrust::vision::{dense_flow_block_match as dense_flow_native, DenseFlowOptions}; use spatialrust::voxelize::{ @@ -3122,18 +3123,106 @@ fn connected_components_image<'py>( Ok((labels.into_pyarray_bound(py), stats)) } +/// Reusable host scratch storage for exact unit-spacing distance transforms. +#[pyclass(name = "DistanceTransformWorkspace")] +struct PyDistanceTransformWorkspace { + inner: DistanceTransformWorkspace, +} + +#[pymethods] +impl PyDistanceTransformWorkspace { + #[new] + fn new() -> Self { + Self { inner: DistanceTransformWorkspace::new() } + } + + #[getter] + fn capacity(&self) -> usize { + self.inner.capacity() + } +} + /// Computes the exact Euclidean distance to the nearest zero-valued mask pixel. #[pyfunction] -#[pyo3(signature = (mask, spacing=(1.0, 1.0)))] +#[pyo3(signature = (mask, spacing=(1.0, 1.0), out=None, workspace=None))] fn distance_transform_edt<'py>( py: Python<'py>, mask: PyReadonlyArray2<'_, u8>, spacing: (f32, f32), + out: Option>>, + workspace: Option>, ) -> PyResult>> { + if workspace.is_some() && spacing != (1.0, 1.0) { + return Err(PyValueError::new_err("workspace reuse currently requires spacing=(1.0, 1.0)")); + } + let mask_view = mask.as_array(); + let (height, width) = (mask_view.shape()[0], mask_view.shape()[1]); + if let Some(mut workspace) = workspace { + let packed; + let input = match mask_view.as_slice() { + Some(slice) => slice, + None => { + packed = mask_view.iter().copied().collect::>(); + packed.as_slice() + } + }; + if let Some(out) = out { + { + let mut out_rw = out.readwrite(); + let mut out_view = out_rw.as_array_mut(); + if out_view.shape() != [height, width] { + return Err(PyValueError::new_err(format!( + "out shape must be ({height}, {width}), found {:?}", + out_view.shape() + ))); + } + let Some(out_slice) = out_view.as_slice_mut() else { + return Err(PyValueError::new_err( + "out must be a contiguous float32 array of shape (H, W)", + )); + }; + distance_transform_edt_u8_into_op( + input, + width, + height, + out_slice, + &mut workspace.inner, + ) + .map_err(to_py_err)?; + } + return Ok(out); + } + let mut output = vec![0.0_f32; width * height]; + distance_transform_edt_u8_into_op(input, width, height, &mut output, &mut workspace.inner) + .map_err(to_py_err)?; + let array = Array2::from_shape_vec((height, width), output).map_err(to_py_err)?; + return Ok(array.into_pyarray_bound(py)); + } + let image = gray_u8_image_from_numpy(mask)?; - let (width, height) = (image.width(), image.height()); let binary = image.into_vec().into_iter().map(|value| u8::from(value != 0)).collect(); let mask = BinaryMask::try_new(width, height, binary).map_err(to_py_err)?; + if let Some(out) = out { + { + let mut out_rw = out.readwrite(); + let mut out_view = out_rw.as_array_mut(); + if out_view.shape() != [height, width] { + return Err(PyValueError::new_err(format!( + "out shape must be ({height}, {width}), found {:?}", + out_view.shape() + ))); + } + let Some(out_slice) = out_view.as_slice_mut() else { + return Err(PyValueError::new_err( + "out must be a contiguous float32 array of shape (H, W)", + )); + }; + let distances = + distance_transform_edt_op(&mask, spacing.0, spacing.1).map_err(to_py_err)?; + out_slice.copy_from_slice(distances.as_slice()); + } + return Ok(out); + } let distances = distance_transform_edt_op(&mask, spacing.0, spacing.1).map_err(to_py_err)?; let array = Array2::from_shape_vec((distances.height(), distances.width()), distances.into_vec()) @@ -3238,6 +3327,7 @@ fn spatialrust_module(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/crates/spatialrust-py/tests/test_bindings.py b/crates/spatialrust-py/tests/test_bindings.py index ddd1118..63b1bc0 100644 --- a/crates/spatialrust-py/tests/test_bindings.py +++ b/crates/spatialrust-py/tests/test_bindings.py @@ -218,6 +218,13 @@ def test_exact_euclidean_distance_transform(): assert anisotropic[0, 1] == pytest.approx(2.0) assert anisotropic[1, 0] == pytest.approx(3.0) + out = np.empty_like(actual) + workspace = sr.DistanceTransformWorkspace() + reused = sr.distance_transform_edt(mask, out=out, workspace=workspace) + assert reused is out + np.testing.assert_array_equal(reused, actual) + assert workspace.capacity >= mask.size + def test_point_map_to_point_cloud_filters_invalid_and_low_confidence(): points = np.array( diff --git a/crates/spatialrust-vision/Cargo.toml b/crates/spatialrust-vision/Cargo.toml index ec3f58f..bdcf7b2 100644 --- a/crates/spatialrust-vision/Cargo.toml +++ b/crates/spatialrust-vision/Cargo.toml @@ -22,7 +22,7 @@ geometry = ["dep:spatialrust-camera"] odometry = ["geometry", "feature2d"] photography = ["warp", "geometry"] detection = [] -dense = ["detection"] +dense = ["detection", "dep:rayon"] spatial = ["dense", "dep:spatialrust-core", "dep:spatialrust-camera"] ai-adapters = ["preprocess", "dense", "detection", "dep:spatialrust-tensor", "dep:bytemuck", "spatialrust-tensor/image"] video = ["dense", "detection"] @@ -37,6 +37,7 @@ spatialrust-camera = { workspace = true, optional = true } spatialrust-tensor = { workspace = true, optional = true } bytemuck = { workspace = true, optional = true } thiserror.workspace = true +rayon = { workspace = true, optional = true } [dev-dependencies] criterion.workspace = true diff --git a/crates/spatialrust-vision/benches/dense.rs b/crates/spatialrust-vision/benches/dense.rs index 6d61752..6af16a7 100644 --- a/crates/spatialrust-vision/benches/dense.rs +++ b/crates/spatialrust-vision/benches/dense.rs @@ -1,5 +1,7 @@ use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput}; -use spatialrust_vision::{distance_transform_edt, BinaryMask}; +use spatialrust_vision::{ + distance_transform_edt, distance_transform_edt_into, BinaryMask, DistanceTransformWorkspace, +}; fn benchmark_exact_distance_transform(c: &mut Criterion) { let mut group = c.benchmark_group("distance_transform_edt"); @@ -13,6 +15,19 @@ fn benchmark_exact_distance_transform(c: &mut Criterion) { group.bench_function(BenchmarkId::from_parameter(name), |b| { b.iter(|| black_box(distance_transform_edt(black_box(&mask)).unwrap())); }); + let mut output = vec![0.0_f32; width * height]; + let mut workspace = DistanceTransformWorkspace::new(); + distance_transform_edt_into(&mask, &mut output, &mut workspace).unwrap(); + group.bench_function(BenchmarkId::new("reuse", name), |b| { + b.iter(|| { + distance_transform_edt_into( + black_box(&mask), + black_box(&mut output), + black_box(&mut workspace), + ) + .unwrap() + }); + }); } group.finish(); } diff --git a/crates/spatialrust-vision/src/dense.rs b/crates/spatialrust-vision/src/dense.rs index 298b247..9552c9c 100644 --- a/crates/spatialrust-vision/src/dense.rs +++ b/crates/spatialrust-vision/src/dense.rs @@ -2,6 +2,7 @@ use std::collections::{BTreeMap, VecDeque}; +use rayon::prelude::*; use spatialrust_image::{ColorSpace, Image, ImageMetadata, ImageView}; use crate::{BoundingBox2, VisionError, VisionResult}; @@ -105,6 +106,102 @@ pub fn distance_transform_edt(mask: &BinaryMask) -> VisionResult> distance_transform_edt_with_spacing(mask, 1.0, 1.0) } +/// Reusable host scratch storage for the unit-spacing exact distance transform. +/// +/// The workspace makes allocation reuse explicit and never shares mutable state +/// between calls. Grow it once for the largest expected image, then pass it to +/// [`distance_transform_edt_into`] on the same worker thread. +#[derive(Debug, Default)] +pub struct DistanceTransformWorkspace { + horizontal: Vec, + transposed: Vec, + column_major: Vec, +} + +impl DistanceTransformWorkspace { + /// Creates an empty workspace that grows on first use. + #[must_use] + pub const fn new() -> Self { + Self { horizontal: Vec::new(), transposed: Vec::new(), column_major: Vec::new() } + } + + /// Returns the largest pixel count currently reserved by every scratch plane. + #[must_use] + pub fn capacity(&self) -> usize { + self.horizontal.capacity().min(self.transposed.capacity()).min(self.column_major.capacity()) + } + + fn prepare(&mut self, len: usize) { + self.horizontal.resize(len, u16::MAX); + self.transposed.resize(len, 0); + self.column_major.resize(len, 0.0); + } +} + +/// Computes the exact unit-spacing Euclidean distance transform into a reusable slice. +/// +/// `output` must contain exactly `mask.width() * mask.height()` elements. Scratch +/// allocations are owned by `workspace`; the caller retains ownership of the +/// output allocation. +pub fn distance_transform_edt_into( + mask: &BinaryMask, + output: &mut [f32], + workspace: &mut DistanceTransformWorkspace, +) -> VisionResult<()> { + distance_transform_edt_u8_into( + mask.image().as_slice(), + mask.width(), + mask.height(), + output, + workspace, + ) +} + +/// Computes a unit-spacing exact distance transform from a packed `u8` mask. +/// +/// Zero is background and every non-zero value is foreground, so common `0/255` +/// image masks need no normalization copy. `input` and `output` must both have +/// exactly `width * height` elements. +pub fn distance_transform_edt_u8_into( + input: &[u8], + width: usize, + height: usize, + output: &mut [f32], + workspace: &mut DistanceTransformWorkspace, +) -> VisionResult<()> { + let len = width.checked_mul(height).ok_or_else(|| { + VisionError::InvalidDimensions("distance-transform size overflows".into()) + })?; + if input.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "distance-transform input needs {len} elements, found {}", + input.len() + ))); + } + if output.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "distance-transform output needs {len} elements, found {}", + output.len() + ))); + } + if len == 0 { + return Ok(()); + } + if !input.contains(&0) { + return Err(VisionError::InvalidParameter( + "distance transform requires at least one background pixel".to_owned(), + )); + } + if width > u16::MAX as usize || height > u16::MAX as usize { + return Err(VisionError::InvalidDimensions( + "reusable unit distance transform supports dimensions up to 65535".to_owned(), + )); + } + workspace.prepare(len); + distance_transform_unit_into(input, width, height, workspace, output); + Ok(()) +} + /// Computes the exact Euclidean distance transform with physical pixel spacing. /// /// `spacing_x` and `spacing_y` are the positive finite distances between @@ -145,44 +242,38 @@ pub fn distance_transform_edt_with_spacing( )); } - let mut horizontal = vec![f64::INFINITY; len]; - let mut source = vec![f64::INFINITY; width.max(height)]; - let mut transformed = vec![0.0_f64; width.max(height)]; - let mut sites = vec![0_usize; width.max(height)]; - let mut boundaries = vec![0.0_f64; width.max(height).saturating_add(1)]; - - for y in 0..height { - for (x, value) in source[..width].iter_mut().enumerate() { - *value = if mask.contains(x, y) { f64::INFINITY } else { 0.0 }; - } - squared_distance_transform_1d( - &source[..width], - f64::from(spacing_x).powi(2), - &mut transformed[..width], - &mut sites[..width], - &mut boundaries[..=width], - ); - horizontal[y * width..(y + 1) * width].copy_from_slice(&transformed[..width]); - } - - let mut squared = vec![0.0_f64; len]; - for x in 0..width { - for y in 0..height { - source[y] = horizontal[y * width + x]; - } - squared_distance_transform_1d( - &source[..height], - f64::from(spacing_y).powi(2), - &mut transformed[..height], - &mut sites[..height], - &mut boundaries[..=height], - ); - for y in 0..height { - squared[y * width + x] = transformed[y]; - } + if spacing_x == 1.0 + && spacing_y == 1.0 + && width <= u16::MAX as usize + && height <= u16::MAX as usize + { + let distances = distance_transform_unit(mask.image().as_slice(), width, height); + return Ok(Image::try_new_with_metadata( + width, + height, + distances, + ImageMetadata { color_space: ColorSpace::Gray, ..Default::default() }, + )?); } - let distances = squared.into_iter().map(|value| value.sqrt() as f32).collect(); + let mut horizontal = vec![f64::INFINITY; len]; + horizontal_binary_distance_rows( + mask.image().as_slice(), + width, + f64::from(spacing_x).powi(2), + &mut horizontal, + ); + + // Transpose in cache-sized tiles so the expensive second pass reads and + // writes each column linearly instead of gathering it with a full-row + // stride. The envelope arithmetic remains f64; only final distances are f32. + let mut transposed = vec![0.0_f64; len]; + transpose_f64(&horizontal, width, height, &mut transposed); + drop(horizontal); + let mut column_major = vec![0.0_f32; len]; + vertical_distance_columns(&transposed, height, f64::from(spacing_y).powi(2), &mut column_major); + let mut distances = vec![0.0_f32; len]; + transpose_f32(&column_major, height, width, &mut distances); Ok(Image::try_new_with_metadata( width, height, @@ -191,29 +282,375 @@ pub fn distance_transform_edt_with_spacing( )?) } -fn squared_distance_transform_1d( - input: &[f64], - coordinate_scale_squared: f64, - output: &mut [f64], +fn distance_transform_unit(mask: &[u8], width: usize, height: usize) -> Vec { + let mut workspace = DistanceTransformWorkspace::new(); + workspace.prepare(mask.len()); + let mut distances = vec![0.0_f32; mask.len()]; + distance_transform_unit_into(mask, width, height, &mut workspace, &mut distances); + distances +} + +fn distance_transform_unit_into( + mask: &[u8], + width: usize, + height: usize, + workspace: &mut DistanceTransformWorkspace, + distances: &mut [f32], +) { + horizontal_binary_distance_rows_u16(mask, width, &mut workspace.horizontal); + transpose_u16(&workspace.horizontal, width, height, &mut workspace.transposed); + vertical_distance_columns_u16(&workspace.transposed, height, &mut workspace.column_major); + transpose_f32(&workspace.column_major, height, width, distances); +} + +fn horizontal_binary_distance_rows_u16(mask: &[u8], width: usize, output: &mut [u16]) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let height = mask.len() / width; + let threads = worker_count(height); + if threads == 1 || mask.len() < PARALLEL_THRESHOLD { + horizontal_binary_distance_rows_u16_serial(mask, width, output); + return; + } + let rows_per_worker = height.div_ceil(threads); + mask.par_chunks(rows_per_worker * width) + .zip(output.par_chunks_mut(rows_per_worker * width)) + .for_each(|(worker_mask, worker_output)| { + horizontal_binary_distance_rows_u16_serial(worker_mask, width, worker_output); + }); +} + +fn horizontal_binary_distance_rows_u16_serial(mask: &[u8], width: usize, output: &mut [u16]) { + for (mask_row, output_row) in mask.chunks_exact(width).zip(output.chunks_exact_mut(width)) { + let mut last_background = None; + for x in 0..width { + if mask_row[x] == 0 { + last_background = Some(x); + output_row[x] = 0; + } else if let Some(background) = last_background { + output_row[x] = (x - background) as u16; + } else { + output_row[x] = u16::MAX; + } + } + let mut next_background = None; + for x in (0..width).rev() { + if mask_row[x] == 0 { + next_background = Some(x); + } else if let Some(background) = next_background { + output_row[x] = output_row[x].min((background - x) as u16); + } + } + } +} + +fn transpose_u16(input: &[u16], width: usize, height: usize, output: &mut [u16]) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let threads = worker_count(width); + if threads == 1 || input.len() < PARALLEL_THRESHOLD { + transpose_u16_columns(input, width, height, 0, output); + return; + } + let columns_per_worker = width.div_ceil(threads); + output.par_chunks_mut(columns_per_worker * height).enumerate().for_each( + |(worker_index, worker_output)| { + let x_start = worker_index * columns_per_worker; + transpose_u16_columns(input, width, height, x_start, worker_output); + }, + ); +} + +fn transpose_u16_columns( + input: &[u16], + width: usize, + height: usize, + x_start: usize, + output: &mut [u16], +) { + const TILE: usize = 32; + let column_count = output.len() / height; + for y_start in (0..height).step_by(TILE) { + for local_x_start in (0..column_count).step_by(TILE) { + let y_end = (y_start + TILE).min(height); + let local_x_end = (local_x_start + TILE).min(column_count); + for y in y_start..y_end { + for local_x in local_x_start..local_x_end { + output[local_x * height + y] = input[y * width + x_start + local_x]; + } + } + } + } +} + +fn vertical_distance_columns_u16(columns: &[u16], height: usize, output: &mut [f32]) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let width = columns.len() / height; + let threads = worker_count(width); + if threads == 1 || columns.len() < PARALLEL_THRESHOLD { + vertical_distance_column_range_u16(columns, height, output); + return; + } + let columns_per_worker = width.div_ceil(threads); + columns + .par_chunks(columns_per_worker * height) + .zip(output.par_chunks_mut(columns_per_worker * height)) + .for_each(|(worker_input, worker_output)| { + vertical_distance_column_range_u16(worker_input, height, worker_output); + }); +} + +fn vertical_distance_column_range_u16(columns: &[u16], height: usize, output: &mut [f32]) { + let mut sites = vec![0_usize; height]; + let mut boundaries = vec![0.0_f64; height.saturating_add(1)]; + for (input_column, output_column) in + columns.chunks_exact(height).zip(output.chunks_exact_mut(height)) + { + distance_transform_1d_u16(input_column, output_column, &mut sites, &mut boundaries); + } +} + +fn distance_transform_1d_u16( + input: &[u16], + output: &mut [f32], sites: &mut [usize], boundaries: &mut [f64], ) { - debug_assert_eq!(input.len(), output.len()); - debug_assert!(sites.len() >= input.len()); - debug_assert!(boundaries.len() > input.len()); - if input.is_empty() { + let first = input.iter().position(|&value| value != u16::MAX).expect("background exists"); + let mut envelope_end = 0; + sites[0] = first; + boundaries[0] = f64::NEG_INFINITY; + boundaries[1] = f64::INFINITY; + for q in first + 1..input.len() { + if input[q] == u16::MAX { + continue; + } + let mut intersection = parabola_intersection_u16(input, q, sites[envelope_end]); + while envelope_end > 0 && intersection <= boundaries[envelope_end] { + envelope_end -= 1; + intersection = parabola_intersection_u16(input, q, sites[envelope_end]); + } + envelope_end += 1; + sites[envelope_end] = q; + boundaries[envelope_end] = intersection; + boundaries[envelope_end + 1] = f64::INFINITY; + } + let mut envelope = 0; + for (q, distance) in output.iter_mut().enumerate() { + while boundaries[envelope + 1] < q as f64 { + envelope += 1; + } + let site = sites[envelope]; + let delta = q.abs_diff(site) as u64; + let horizontal = u64::from(input[site]); + *distance = ((delta * delta + horizontal * horizontal) as f32).sqrt(); + } +} + +fn parabola_intersection_u16(input: &[u16], right: usize, left: usize) -> f64 { + let right_coordinate = right as f64; + let left_coordinate = left as f64; + let right_distance = f64::from(input[right]); + let left_distance = f64::from(input[left]); + let right_height = right_coordinate * right_coordinate + right_distance * right_distance; + let left_height = left_coordinate * left_coordinate + left_distance * left_distance; + (right_height - left_height) / (2.0 * (right_coordinate - left_coordinate)) +} + +fn transpose_f64(input: &[f64], width: usize, height: usize, output: &mut [f64]) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let threads = worker_count(width); + if threads == 1 || input.len() < PARALLEL_THRESHOLD { + transpose_f64_columns(input, width, height, 0, output); + return; + } + let columns_per_worker = width.div_ceil(threads); + output.par_chunks_mut(columns_per_worker * height).enumerate().for_each( + |(worker_index, worker_output)| { + let x_start = worker_index * columns_per_worker; + transpose_f64_columns(input, width, height, x_start, worker_output); + }, + ); +} + +fn transpose_f64_columns( + input: &[f64], + width: usize, + height: usize, + x_start: usize, + output: &mut [f64], +) { + const TILE: usize = 32; + let column_count = output.len() / height; + for y_start in (0..height).step_by(TILE) { + for local_x_start in (0..column_count).step_by(TILE) { + let y_end = (y_start + TILE).min(height); + let local_x_end = (local_x_start + TILE).min(column_count); + for y in y_start..y_end { + for local_x in local_x_start..local_x_end { + output[local_x * height + y] = input[y * width + x_start + local_x]; + } + } + } + } +} + +fn transpose_f32(input: &[f32], width: usize, height: usize, output: &mut [f32]) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let threads = worker_count(width); + if threads == 1 || input.len() < PARALLEL_THRESHOLD { + transpose_f32_columns(input, width, height, 0, output); return; } + let columns_per_worker = width.div_ceil(threads); + output.par_chunks_mut(columns_per_worker * height).enumerate().for_each( + |(worker_index, worker_output)| { + let x_start = worker_index * columns_per_worker; + transpose_f32_columns(input, width, height, x_start, worker_output); + }, + ); +} + +fn transpose_f32_columns( + input: &[f32], + width: usize, + height: usize, + x_start: usize, + output: &mut [f32], +) { + const TILE: usize = 32; + let column_count = output.len() / height; + for y_start in (0..height).step_by(TILE) { + for local_x_start in (0..column_count).step_by(TILE) { + let y_end = (y_start + TILE).min(height); + let local_x_end = (local_x_start + TILE).min(column_count); + for y in y_start..y_end { + for local_x in local_x_start..local_x_end { + output[local_x * height + y] = input[y * width + x_start + local_x]; + } + } + } + } +} +fn horizontal_binary_distance_rows( + mask: &[u8], + width: usize, + spacing_squared: f64, + output: &mut [f64], +) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let height = mask.len() / width; + let threads = worker_count(height); + if threads == 1 || mask.len() < PARALLEL_THRESHOLD { + horizontal_binary_distance_rows_serial(mask, width, spacing_squared, output); + return; + } + let rows_per_worker = height.div_ceil(threads); + mask.par_chunks(rows_per_worker * width) + .zip(output.par_chunks_mut(rows_per_worker * width)) + .for_each(|(worker_mask, worker_output)| { + horizontal_binary_distance_rows_serial( + worker_mask, + width, + spacing_squared, + worker_output, + ); + }); +} + +fn horizontal_binary_distance_rows_serial( + mask: &[u8], + width: usize, + spacing_squared: f64, + output: &mut [f64], +) { + for (mask_row, output_row) in mask.chunks_exact(width).zip(output.chunks_exact_mut(width)) { + let mut last_background = None; + for x in 0..width { + if mask_row[x] == 0 { + last_background = Some(x); + output_row[x] = 0.0; + } else if let Some(background) = last_background { + let delta = (x - background) as f64; + output_row[x] = spacing_squared * delta * delta; + } + } + + let mut next_background = None; + for x in (0..width).rev() { + if mask_row[x] == 0 { + next_background = Some(x); + } else if let Some(background) = next_background { + let delta = (background - x) as f64; + output_row[x] = output_row[x].min(spacing_squared * delta * delta); + } + } + } +} + +fn vertical_distance_columns( + columns: &[f64], + height: usize, + spacing_squared: f64, + output: &mut [f32], +) { + const PARALLEL_THRESHOLD: usize = 128 * 1024; + let width = columns.len() / height; + let threads = worker_count(width); + if threads == 1 || columns.len() < PARALLEL_THRESHOLD { + vertical_distance_column_range(columns, height, spacing_squared, output); + return; + } + + let columns_per_worker = width.div_ceil(threads); + columns + .par_chunks(columns_per_worker * height) + .zip(output.par_chunks_mut(columns_per_worker * height)) + .for_each(|(worker_input, worker_output)| { + vertical_distance_column_range(worker_input, height, spacing_squared, worker_output); + }); +} + +fn worker_count(independent_items: usize) -> usize { + rayon::current_num_threads().min(independent_items) +} + +fn vertical_distance_column_range( + columns: &[f64], + height: usize, + spacing_squared: f64, + output: &mut [f32], +) { + let mut sites = vec![0_usize; height]; + let mut boundaries = vec![0.0_f64; height.saturating_add(1)]; + for (input_column, output_column) in + columns.chunks_exact(height).zip(output.chunks_exact_mut(height)) + { + distance_transform_1d_to_f32( + input_column, + spacing_squared, + output_column, + &mut sites, + &mut boundaries, + ); + } +} + +fn distance_transform_1d_to_f32( + input: &[f64], + coordinate_scale_squared: f64, + output: &mut [f32], + sites: &mut [usize], + boundaries: &mut [f64], +) { + debug_assert_eq!(input.len(), output.len()); let Some(first) = input.iter().position(|value| value.is_finite()) else { - output.fill(f64::INFINITY); + output.fill(f32::INFINITY); return; }; - let mut envelope_end = 0_usize; + let mut envelope_end = 0; sites[0] = first; boundaries[0] = f64::NEG_INFINITY; boundaries[1] = f64::INFINITY; - for q in first + 1..input.len() { if !input[q].is_finite() { continue; @@ -231,14 +668,14 @@ fn squared_distance_transform_1d( boundaries[envelope_end + 1] = f64::INFINITY; } - let mut envelope = 0_usize; - for (q, value) in output.iter_mut().enumerate() { + let mut envelope = 0; + for (q, distance) in output.iter_mut().enumerate() { while boundaries[envelope + 1] < q as f64 { envelope += 1; } let site = sites[envelope]; let delta = q as f64 - site as f64; - *value = coordinate_scale_squared.mul_add(delta * delta, input[site]); + *distance = coordinate_scale_squared.mul_add(delta * delta, input[site]).sqrt() as f32; } } @@ -798,8 +1235,9 @@ impl PointMap { mod tests { use super::{ approximate_polygon, connected_components, decode_rle, distance_transform_edt, - distance_transform_edt_with_spacing, encode_rle, find_contours, BinaryMask, Connectivity, - DepthMap, FlowField, PointMap, RleOrder, + distance_transform_edt_u8_into, distance_transform_edt_with_spacing, encode_rle, + find_contours, BinaryMask, Connectivity, DepthMap, DistanceTransformWorkspace, FlowField, + PointMap, RleOrder, }; use spatialrust_image::Image; @@ -871,6 +1309,45 @@ mod tests { } } + #[test] + fn unit_distance_transform_matches_brute_force_on_irregular_mask() { + let (width, height) = (31, 23); + let data = (0..width * height) + .map(|index| u8::from(index % 37 != 0 && index % 101 != 7)) + .collect::>(); + let background = data + .iter() + .enumerate() + .filter(|(_, value)| **value == 0) + .map(|(index, _)| (index % width, index / width)) + .collect::>(); + let mask = BinaryMask::try_new(width, height, data).unwrap(); + let actual = distance_transform_edt(&mask).unwrap(); + for y in 0..height { + for x in 0..width { + let expected = background + .iter() + .map(|&(bx, by)| (x.abs_diff(bx) as f32).hypot(y.abs_diff(by) as f32)) + .fold(f32::INFINITY, f32::min); + assert_eq!(actual[(x, y)][0], expected); + } + } + } + + #[test] + fn reusable_distance_transform_accepts_packed_255_masks() { + let input = [0_u8, 255, 255, 255, 255, 255]; + let mut output = [0.0_f32; 6]; + let mut workspace = DistanceTransformWorkspace::new(); + distance_transform_edt_u8_into(&input, 3, 2, &mut output, &mut workspace).unwrap(); + assert_eq!(output, [0.0, 1.0, 2.0, 1.0, 2.0_f32.sqrt(), 5.0_f32.sqrt()]); + let capacity = workspace.capacity(); + distance_transform_edt_u8_into(&input, 3, 2, &mut output, &mut workspace).unwrap(); + assert_eq!(workspace.capacity(), capacity); + assert!(distance_transform_edt_u8_into(&input, 3, 2, &mut output[..5], &mut workspace,) + .is_err()); + } + #[test] fn distance_transform_defines_empty_and_rejects_missing_background() { let empty = BinaryMask::try_new(0, 0, Vec::new()).unwrap(); diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index f39e35a..e022472 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -606,6 +606,7 @@ to one implicitly, and GPU receipts must retain named upload/readback stages. | 113B | Planned | Explicit reusable scratch storage for multi-pass algorithms | steady-state allocation receipt | | 113C | Planned | Validate dimensions, metadata, overlap, and channel contracts | negative and property tests | | 113D | Planned | Reuse outputs through Python `out=` where supported | object-identity and numerical tests | +| 113E | Complete | Exact EDT caller-owned output and explicit reusable scratch | Rust/Python identity, capacity, and brute-force tests | ### Epic 114 delivery slices @@ -615,6 +616,7 @@ to one implicitly, and GPU receipts must retain named upload/readback stages. | 114B | Planned | Packed `u8` one/three-channel and `f32` internal fast-path selection | dispatch receipt and fallback parity | | 114C | Planned | Preserve generic components, channels, strides, and borders as safe fallbacks | full property suite | | 114D | Planned | Bound worker creation and temporary memory | thread-count and peak-memory receipt | +| 114E | Complete | Exact EDT binary-row fast path, tiled transpose, and bounded pool dispatch | VGA/1080p/4K Criterion and OpenCV receipt | ### Epic 115 delivery slices diff --git a/docs/site/algorithms.html b/docs/site/algorithms.html index 6d212e6..ac52cc0 100644 --- a/docs/site/algorithms.html +++ b/docs/site/algorithms.html @@ -35,7 +35,7 @@

Algorithm catalog

MorphologyErode, dilate, open, close, gradient, top-hat, black-hat with Rect/Cross/Ellipse/Diamond/custom elementsspatialrust-vision · imgproc-morphologyCPU / GPU Image analysisFixed/Otsu/adaptive threshold, histogram, equalization, CLAHE, integral image, Cannyspatialrust-vision · imgproc-analysis, imgproc-cannyCPU Local featuresFAST, Harris, Shi–Tomasi, ORB, descriptor matching, grid selection, pyramidal Lucas–Kanade trackingspatialrust-vision · feature2dCPU - Dense visionExact Euclidean distance transform, connected components, contours, polygon approximation, mask RLE, depth/confidence/flow/point maps, NMS/Soft-NMSspatialrust-vision · dense, detectionCPU + Dense visionExact Euclidean distance transform with reusable workspace/output, connected components, contours, polygon approximation, mask RLE, depth/confidence/flow/point maps, NMS/Soft-NMSspatialrust-vision · dense, detectionCPU parallel Multiview geometryHomography, fundamental/essential matrices, RANSAC, triangulation, relative pose, PnP/PnP-RANSACspatialrust-vision · geometryCPU Stereo and odometryStereo rectification, block matching, disparity-to-depth/XYZ, monocular and RGB-D visual odometryspatialrust-vision · geometry, odometryCPU VideoDense block flow, adaptive background model, multi-object tracker, timestamped pull sourcesspatialrust-vision · videoCPU diff --git a/notes/2026-07-15_exact_edt_acceleration.md b/notes/2026-07-15_exact_edt_acceleration.md new file mode 100644 index 0000000..915cf11 --- /dev/null +++ b/notes/2026-07-15_exact_edt_acceleration.md @@ -0,0 +1,39 @@ +# Exact EDT acceleration receipt — 2026-07-15 + +## Outcome + +The exact unit-spacing Euclidean distance transform now uses a binary-row +nearest-background pass, compact `u16` horizontal distances, cache-tiled +transposes, a parallel Felzenszwalb–Huttenlocher column envelope, and final +`f32` square roots. Anisotropic spacing keeps the general exact `f64` path. + +`DistanceTransformWorkspace`, `distance_transform_edt_into`, and +`distance_transform_edt_u8_into` make scratch/output reuse explicit. The Python +binding exposes the same policy through `DistanceTransformWorkspace` and +`out=`; common packed `0/255` masks no longer require a normalization copy. + +## Correctness + +- canonical VGA, 1080p, and 4K masks: exact fraction `1.0`, max error `0.0` + against OpenCV `DIST_L2/DIST_MASK_PRECISE`; +- irregular unit mask: exact brute-force equality; +- anisotropic property coverage retains the general-spacing implementation; +- non-empty all-foreground masks remain an explicit error. + +## Performance + +Host: Windows 11, 6-core/12-thread Intel CPU, OpenCV 4.10, OpenCL disabled, +CPython 3.12. Timings are machine-specific medians. + +| Measurement | Before | After | +| --- | ---: | ---: | +| Native Criterion 4K allocate | 451.63 ms | about 75 ms | +| Native Criterion 4K reusable output/workspace | unavailable | about 43 ms | +| Python/OpenCV 4K allocate ratio | OpenCV 12.35× | OpenCV 1.63× | +| Python/OpenCV 4K reuse ratio | unavailable | OpenCV 1.07× | + +The native reusable kernel crosses the earlier OpenCV allocate baseline, while +the fully interleaved Python reuse comparison remains a narrow OpenCV win. No +blanket faster-than-OpenCV claim is made. Re-run +`bench/opencv_vision_comparison/performance.py` before quoting host-specific +numbers.