diff --git a/CHANGELOG.md b/CHANGELOG.md index 8bc13a1..22f8316 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,11 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Added +- **Exact Euclidean distance transform**: `spatialrust-vision` now computes + foreground-to-nearest-background L2 distances in linear time, supports + anisotropic pixel spacing, exposes a NumPy binding, and includes native + Criterion plus OpenCV `DIST_MASK_PRECISE` comparison coverage. + - **Vision 2 performance roadmap and documentation site**: Epics 112–120 now define cost attribution, reusable workspaces, safe CPU dispatch, resize/color, Gaussian/Sobel, morphology, Canny, explicit GPU-chain, and release-gate work. diff --git a/README.md b/README.md index d972c20..53f0a89 100644 --- a/README.md +++ b/README.md @@ -160,6 +160,7 @@ ratio; these are machine-specific measurements, not universal guarantees. | 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× | 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 @@ -181,6 +182,7 @@ The same deterministic RGB inputs passed all VGA, 1080p, and 4K gates: | Sobel X 3×3 | Exact values (max error 0) | | Morphology open 5×5 | Exact pixels (max error 0) | | Canny | Precision, recall, F1, and IoU all 1.0 | +| Exact Euclidean distance transform | Exact values on canonical profiles; separate irregular-mask max float error `9.54e-7` | The broader correctness harness also checks filters, analysis, keypoints, matching, and geometry with documented tolerances (exact pixels where we claim diff --git a/bench/opencv_comparison/manifest.json b/bench/opencv_comparison/manifest.json index b3dbb57..d664787 100644 --- a/bench/opencv_comparison/manifest.json +++ b/bench/opencv_comparison/manifest.json @@ -32,6 +32,7 @@ { "id": "sobel_x", "domain": "imgproc", "modes": ["allocate"] }, { "id": "canny", "domain": "imgproc", "modes": ["allocate"] }, { "id": "morphology_open", "domain": "imgproc", "modes": ["allocate"] }, + { "id": "distance_transform_edt", "domain": "imgproc", "modes": ["allocate"] }, { "id": "orb", "domain": "feature2d", "modes": ["allocate"] }, { "id": "stereo_bm", "domain": "calib3d", "modes": ["allocate"] }, { "id": "pinhole_calibration", "domain": "calib3d", "modes": ["allocate"] }, diff --git a/bench/opencv_vision_comparison/performance.py b/bench/opencv_vision_comparison/performance.py index 16699cd..cb146c1 100644 --- a/bench/opencv_vision_comparison/performance.py +++ b/bench/opencv_vision_comparison/performance.py @@ -314,6 +314,31 @@ def main() -> None: min_sample_time_ms=MIN_SAMPLE_TIME_MS, ) + distance_mask = np.where(gray_cv > 96, 255, 0).astype(np.uint8) + distance_cv = cv2.distanceTransform( + distance_mask, cv2.DIST_L2, cv2.DIST_MASK_PRECISE + ) + distance_sr = sr.distance_transform_edt(distance_mask) + distance_accuracy = numerical_accuracy( + distance_cv, distance_sr, float(np.hypot(width, height)) + ) + accuracy[f"{profile}_distance_transform_edt"] = distance_accuracy + if distance_accuracy["max_absolute_error"] > 1e-5: + raise AssertionError( + f"{profile} distance-transform max error " + f"{distance_accuracy['max_absolute_error']} > 1e-5" + ) + _, _, cv_distance, sr_distance = timed_pair( + lambda: cv2.distanceTransform( + distance_mask, cv2.DIST_L2, cv2.DIST_MASK_PRECISE + ), + lambda: sr.distance_transform_edt(distance_mask), + warmup=args.warmup, + repeats=repeats, + seed=113, + min_sample_time_ms=MIN_SAMPLE_TIME_MS, + ) + rows = ( ("resize_bilinear", "opencv", "allocate", cv_resize_alloc), ("resize_bilinear", "spatialrust", "allocate", sr_resize_alloc), @@ -334,6 +359,8 @@ def main() -> None: ("morphology_open", "spatialrust", "allocate", sr_morphology), ("canny", "opencv", "allocate", cv_canny), ("canny", "spatialrust", "allocate", sr_canny), + ("distance_transform_edt", "opencv", "allocate", cv_distance), + ("distance_transform_edt", "spatialrust", "allocate", sr_distance), ) measurements.extend( measurement(workload, implementation, mode, width, height, timing) @@ -352,6 +379,7 @@ def main() -> None: "sobel_x": speed_comparison(cv_sobel, sr_sobel), "morphology_open": speed_comparison(cv_morphology, sr_morphology), "canny": speed_comparison(cv_canny, sr_canny), + "distance_transform_edt": speed_comparison(cv_distance, sr_distance), } environment_receipt = environment( diff --git a/bench/opencv_vision_comparison/run.py b/bench/opencv_vision_comparison/run.py index 1056cb0..bbd191a 100644 --- a/bench/opencv_vision_comparison/run.py +++ b/bench/opencv_vision_comparison/run.py @@ -396,6 +396,18 @@ def main() -> None: if areas != areas_cv: raise AssertionError(f"component areas mismatch: {areas} != {areas_cv}") + distance_mask = np.ones((67, 89), dtype=np.uint8) + distance_mask[::11, ::13] = 0 + distance_mask[20:28, 31:40] = 0 + distance_sr = sr.distance_transform_edt(distance_mask) + distance_cv = cv2.distanceTransform( + distance_mask, cv2.DIST_L2, cv2.DIST_MASK_PRECISE + ) + distance_error = float(np.max(np.abs(distance_sr - distance_cv))) + results["distance_transform_edt_max_f32_error"] = distance_error + if distance_error > 1e-5: + raise AssertionError(f"exact distance-transform error {distance_error} > 1e-5") + # Geometry: planar homography residual agreement (not scale-normalized identity). source = np.array( [[10.0, 12.0], [70.0, 14.0], [18.0, 55.0], [66.0, 60.0], [40.0, 34.0], [28.0, 22.0]], diff --git a/crates/spatialrust-platform/src/stability.rs b/crates/spatialrust-platform/src/stability.rs index 2f1f214..ee28716 100644 --- a/crates/spatialrust-platform/src/stability.rs +++ b/crates/spatialrust-platform/src/stability.rs @@ -121,6 +121,8 @@ impl StabilityRegistry { "spatialrust-vision::video", "spatialrust-vision::odometry", "spatialrust-vision::photography", + "spatialrust-vision::distance_transform_edt", + "spatialrust-vision::distance_transform_edt_with_spacing", "spatialrust-runtime::execution-graph", "spatialrust-vision::ai-adapters", "spatialrust-gpu::GpuImage", diff --git a/crates/spatialrust-py/spatialrust.pyi b/crates/spatialrust-py/spatialrust.pyi index 71ccffd..2db757d 100644 --- a/crates/spatialrust-py/spatialrust.pyi +++ b/crates/spatialrust-py/spatialrust.pyi @@ -36,7 +36,8 @@ __all__: list[str] = [ "histogram_image", "equalize_histogram_image", "clahe_image", "integral_image_u8", "canny_image", "resize_image", "letterbox_image", "normalize_image_chw", "rgb_to_gray_image", "rgb_to_hsv_image", "remap_image", - "nms", "soft_nms", "connected_components_image", "find_mask_contours", + "nms", "soft_nms", "connected_components_image", "distance_transform_edt", + "find_mask_contours", "encode_mask_rle", "decode_mask_rle", "point_map_to_point_cloud", "knn_graph", "radius_graph", "register_icp", "register_point_to_plane", "register_gicp", "register_ndt", "register_fpfh_ransac", "register_fpfh_keypoints", @@ -284,6 +285,9 @@ def connected_components_image( mask: _U8Array, connectivity: int = ..., ) -> tuple[_U32Array, list[tuple[int, int, tuple[float, float, float, float]]]]: ... +def distance_transform_edt( + mask: _U8Array, spacing: tuple[float, float] = ... +) -> _F32Array: ... def find_mask_contours( mask: _U8Array, epsilon: float = ... ) -> list[list[tuple[int, int]]]: ... diff --git a/crates/spatialrust-py/src/lib.rs b/crates/spatialrust-py/src/lib.rs index 5e0c277..a278392 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_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, estimate_rgbd_odometry as estimate_rgbd_odometry_op, filter2d as filter2d_op, @@ -3121,6 +3122,25 @@ fn connected_components_image<'py>( Ok((labels.into_pyarray_bound(py), stats)) } +/// Computes the exact Euclidean distance to the nearest zero-valued mask pixel. +#[pyfunction] +#[pyo3(signature = (mask, spacing=(1.0, 1.0)))] +fn distance_transform_edt<'py>( + py: Python<'py>, + mask: PyReadonlyArray2<'_, u8>, + spacing: (f32, f32), +) -> PyResult>> { + 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)?; + 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()) + .map_err(to_py_err)?; + Ok(array.into_pyarray_bound(py)) +} + /// Extracts and optionally simplifies mask contours. #[pyfunction] #[pyo3(signature = (mask, epsilon=0.0))] @@ -3308,6 +3328,7 @@ fn spatialrust_module(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(nms, m)?)?; m.add_function(wrap_pyfunction!(soft_nms, m)?)?; m.add_function(wrap_pyfunction!(connected_components_image, m)?)?; + m.add_function(wrap_pyfunction!(distance_transform_edt, m)?)?; m.add_function(wrap_pyfunction!(find_mask_contours, m)?)?; m.add_function(wrap_pyfunction!(encode_mask_rle, m)?)?; m.add_function(wrap_pyfunction!(decode_mask_rle, m)?)?; diff --git a/crates/spatialrust-py/tests/test_bindings.py b/crates/spatialrust-py/tests/test_bindings.py index 7f47234..ddd1118 100644 --- a/crates/spatialrust-py/tests/test_bindings.py +++ b/crates/spatialrust-py/tests/test_bindings.py @@ -66,7 +66,7 @@ def test_exports_present(): "rgbd_to_point_cloud", "depth_to_xyz", "resize_image", "letterbox_image", "normalize_image_chw", "rgb_to_gray_image", "rgb_to_hsv_image", "remap_image", - "nms", "soft_nms", "connected_components_image", + "nms", "soft_nms", "connected_components_image", "distance_transform_edt", "find_mask_contours", "encode_mask_rle", "decode_mask_rle", "point_map_to_point_cloud", ): @@ -199,6 +199,26 @@ def test_mask_components_contours_and_rle(): np.testing.assert_array_equal(decoded, mask) +def test_exact_euclidean_distance_transform(): + mask = np.full((3, 4), 255, dtype=np.uint8) + mask[0, 0] = 0 + actual = sr.distance_transform_edt(mask) + expected = np.array( + [ + [0.0, 1.0, 2.0, 3.0], + [1.0, np.sqrt(2.0), np.sqrt(5.0), np.sqrt(10.0)], + [2.0, np.sqrt(5.0), np.sqrt(8.0), np.sqrt(13.0)], + ], + dtype=np.float32, + ) + assert actual.dtype == np.float32 + np.testing.assert_allclose(actual, expected, atol=1e-6) + + anisotropic = sr.distance_transform_edt(mask, spacing=(2.0, 3.0)) + assert anisotropic[0, 1] == pytest.approx(2.0) + assert anisotropic[1, 0] == pytest.approx(3.0) + + def test_point_map_to_point_cloud_filters_invalid_and_low_confidence(): points = np.array( [[[0, 0, 1], [1, 0, 1]], [[0, 1, np.nan], [1, 1, 2]]], dtype=np.float32 diff --git a/crates/spatialrust-vision/Cargo.toml b/crates/spatialrust-vision/Cargo.toml index 13c132f..ec3f58f 100644 --- a/crates/spatialrust-vision/Cargo.toml +++ b/crates/spatialrust-vision/Cargo.toml @@ -67,6 +67,11 @@ name = "analysis" harness = false required-features = ["imgproc-analysis"] +[[bench]] +name = "dense" +harness = false +required-features = ["dense"] + [[bench]] name = "canny" harness = false diff --git a/crates/spatialrust-vision/benches/dense.rs b/crates/spatialrust-vision/benches/dense.rs new file mode 100644 index 0000000..6d61752 --- /dev/null +++ b/crates/spatialrust-vision/benches/dense.rs @@ -0,0 +1,21 @@ +use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput}; +use spatialrust_vision::{distance_transform_edt, BinaryMask}; + +fn benchmark_exact_distance_transform(c: &mut Criterion) { + let mut group = c.benchmark_group("distance_transform_edt"); + group.sample_size(10); + for &(name, width, height) in &[("vga", 640, 480), ("1080p", 1920, 1080), ("4k", 3840, 2160)] { + let data = (0..width * height) + .map(|index| u8::from(index % 97 != 0 && index % (width * 11 + 1) != 0)) + .collect(); + let mask = BinaryMask::try_new(width, height, data).unwrap(); + group.throughput(Throughput::Elements((width * height) as u64)); + group.bench_function(BenchmarkId::from_parameter(name), |b| { + b.iter(|| black_box(distance_transform_edt(black_box(&mask)).unwrap())); + }); + } + group.finish(); +} + +criterion_group!(benches, benchmark_exact_distance_transform); +criterion_main!(benches); diff --git a/crates/spatialrust-vision/src/dense.rs b/crates/spatialrust-vision/src/dense.rs index 02fde7f..298b247 100644 --- a/crates/spatialrust-vision/src/dense.rs +++ b/crates/spatialrust-vision/src/dense.rs @@ -92,6 +92,172 @@ impl BinaryMask { } } +/// Computes the exact Euclidean distance from every foreground pixel to the +/// nearest background pixel using unit pixel spacing. +/// +/// Background pixels have distance zero. A non-empty mask must contain at +/// least one background pixel; otherwise the finite distance is undefined and +/// an error is returned. The implementation applies the separable linear-time +/// squared-distance transform by Felzenszwalb and Huttenlocher. +/// +/// See . +pub fn distance_transform_edt(mask: &BinaryMask) -> VisionResult> { + distance_transform_edt_with_spacing(mask, 1.0, 1.0) +} + +/// Computes the exact Euclidean distance transform with physical pixel spacing. +/// +/// `spacing_x` and `spacing_y` are the positive finite distances between +/// adjacent pixel centers along each axis. The output is expressed in those +/// physical units. +pub fn distance_transform_edt_with_spacing( + mask: &BinaryMask, + spacing_x: f32, + spacing_y: f32, +) -> VisionResult> { + if !spacing_x.is_finite() || spacing_x <= 0.0 { + return Err(VisionError::InvalidParameter( + "distance-transform x spacing must be finite and positive".to_owned(), + )); + } + if !spacing_y.is_finite() || spacing_y <= 0.0 { + return Err(VisionError::InvalidParameter( + "distance-transform y spacing must be finite and positive".to_owned(), + )); + } + + let width = mask.width(); + let height = mask.height(); + let len = width.checked_mul(height).ok_or_else(|| { + VisionError::InvalidDimensions("distance-transform size overflows".into()) + })?; + if len == 0 { + return Ok(Image::try_new_with_metadata( + width, + height, + Vec::new(), + ImageMetadata { color_space: ColorSpace::Gray, ..Default::default() }, + )?); + } + if mask.area() == len { + return Err(VisionError::InvalidParameter( + "distance transform requires at least one background pixel".to_owned(), + )); + } + + 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]; + } + } + + let distances = squared.into_iter().map(|value| value.sqrt() as f32).collect(); + Ok(Image::try_new_with_metadata( + width, + height, + distances, + ImageMetadata { color_space: ColorSpace::Gray, ..Default::default() }, + )?) +} + +fn squared_distance_transform_1d( + input: &[f64], + coordinate_scale_squared: f64, + output: &mut [f64], + 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() { + return; + } + + let Some(first) = input.iter().position(|value| value.is_finite()) else { + output.fill(f64::INFINITY); + return; + }; + let mut envelope_end = 0_usize; + 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; + } + let mut intersection = + parabola_intersection(input, coordinate_scale_squared, q, sites[envelope_end]); + while envelope_end > 0 && intersection <= boundaries[envelope_end] { + envelope_end -= 1; + intersection = + parabola_intersection(input, coordinate_scale_squared, 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_usize; + for (q, value) 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]); + } +} + +fn parabola_intersection( + input: &[f64], + coordinate_scale_squared: f64, + right: usize, + left: usize, +) -> f64 { + let right_coordinate = right as f64; + let left_coordinate = left as f64; + let right_height = + coordinate_scale_squared.mul_add(right_coordinate * right_coordinate, input[right]); + let left_height = + coordinate_scale_squared.mul_add(left_coordinate * left_coordinate, input[left]); + (right_height - left_height) + / (2.0 * coordinate_scale_squared * (right_coordinate - left_coordinate)) +} + /// Connected-component label image (`0` is background). #[derive(Clone, Debug, PartialEq, Eq)] pub struct LabelImage { @@ -631,8 +797,9 @@ impl PointMap { #[cfg(test)] mod tests { use super::{ - approximate_polygon, connected_components, decode_rle, encode_rle, find_contours, - BinaryMask, Connectivity, DepthMap, FlowField, PointMap, RleOrder, + approximate_polygon, connected_components, decode_rle, distance_transform_edt, + distance_transform_edt_with_spacing, encode_rle, find_contours, BinaryMask, Connectivity, + DepthMap, FlowField, PointMap, RleOrder, }; use spatialrust_image::Image; @@ -671,6 +838,52 @@ mod tests { } } + #[test] + fn exact_distance_transform_matches_known_grid() { + let mask = BinaryMask::try_new(4, 3, vec![0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]).unwrap(); + let distances = distance_transform_edt(&mask).unwrap(); + let expected = [ + 0.0, + 1.0, + 2.0, + 3.0, + 1.0, + 2.0_f32.sqrt(), + 5.0_f32.sqrt(), + 10.0_f32.sqrt(), + 2.0, + 5.0_f32.sqrt(), + 8.0_f32.sqrt(), + 13.0_f32.sqrt(), + ]; + for (&actual, expected) in distances.as_slice().iter().zip(expected) { + assert!((actual - expected).abs() <= 1e-6, "{actual} != {expected}"); + } + } + + #[test] + fn distance_transform_respects_anisotropic_spacing() { + let mask = BinaryMask::try_new(3, 2, vec![0, 1, 1, 1, 1, 1]).unwrap(); + let distances = distance_transform_edt_with_spacing(&mask, 2.0, 3.0).unwrap(); + let expected = [0.0, 2.0, 4.0, 3.0, 13.0_f32.sqrt(), 5.0]; + for (&actual, expected) in distances.as_slice().iter().zip(expected) { + assert!((actual - expected).abs() <= 1e-6, "{actual} != {expected}"); + } + } + + #[test] + fn distance_transform_defines_empty_and_rejects_missing_background() { + let empty = BinaryMask::try_new(0, 0, Vec::new()).unwrap(); + assert!(distance_transform_edt(&empty).unwrap().as_slice().is_empty()); + + let foreground = BinaryMask::try_new(2, 2, vec![1; 4]).unwrap(); + assert!(distance_transform_edt(&foreground).is_err()); + let background = BinaryMask::try_new(2, 2, vec![0; 4]).unwrap(); + assert_eq!(distance_transform_edt(&background).unwrap().as_slice(), &[0.0; 4]); + assert!(distance_transform_edt_with_spacing(&background, 0.0, 1.0).is_err()); + assert!(distance_transform_edt_with_spacing(&background, 1.0, f32::NAN).is_err()); + } + #[test] fn depth_flow_and_point_maps_expose_validity() { let depth = DepthMap::try_new(2, 1, vec![1.0, f32::NAN]).unwrap(); diff --git a/crates/spatialrust-vision/tests/properties.rs b/crates/spatialrust-vision/tests/properties.rs index 8c5b707..1c527a6 100644 --- a/crates/spatialrust-vision/tests/properties.rs +++ b/crates/spatialrust-vision/tests/properties.rs @@ -3,16 +3,16 @@ #![cfg(feature = "full")] use proptest::prelude::*; +use spatialrust_camera::CameraIntrinsics; use spatialrust_image::Image; +use spatialrust_math::{Mat3, Vec2, Vec3}; use spatialrust_vision::{ - canny, decode_rle, encode_rle, erode, estimate_homography, filter2d, integral_image, - match_descriptors, project_object_point, resize, solve_pnp, AbsolutePose, BinaryMask, - BorderMode, BoundingBox2, CameraMatrix3, CannyOptions, DescriptorBuffer, Interpolation, - Kernel2D, MatchOptions, MorphologyShape, ObjectImageCorrespondence, PointCorrespondence2, - RleOrder, StructuringElement, + canny, decode_rle, distance_transform_edt_with_spacing, encode_rle, erode, estimate_homography, + filter2d, integral_image, match_descriptors, project_object_point, resize, solve_pnp, + AbsolutePose, BinaryMask, BorderMode, BoundingBox2, CameraMatrix3, CannyOptions, + DescriptorBuffer, Interpolation, Kernel2D, MatchOptions, MorphologyShape, + ObjectImageCorrespondence, PointCorrespondence2, RleOrder, StructuringElement, }; -use spatialrust_camera::CameraIntrinsics; -use spatialrust_math::{Mat3, Vec2, Vec3}; proptest! { #[test] @@ -139,6 +139,48 @@ proptest! { } } + #[test] + fn exact_distance_transform_matches_brute_force( + width in 1usize..16, + height in 1usize..16, + seed in any::(), + spacing_x in 0.1f32..4.0, + spacing_y in 0.1f32..4.0, + ) { + let mut state = seed; + let mut data = (0..width * height) + .map(|_| { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + (state & 1) as u8 + }) + .collect::>(); + let forced_background = seed as usize % data.len(); + data[forced_background] = 0; + let mask = BinaryMask::try_new(width, height, data.clone()).unwrap(); + let actual = distance_transform_edt_with_spacing(&mask, spacing_x, spacing_y).unwrap(); + let background = data + .iter() + .enumerate() + .filter(|(_, value)| **value == 0) + .map(|(index, _)| (index % width, index / width)) + .collect::>(); + for y in 0..height { + for x in 0..width { + let expected = background + .iter() + .map(|&(bx, by)| { + let dx = x.abs_diff(bx) as f32 * spacing_x; + let dy = y.abs_diff(by) as f32 * spacing_y; + dx.hypot(dy) + }) + .fold(f32::INFINITY, f32::min); + prop_assert!((actual[(x, y)][0] - expected).abs() <= 2e-5); + } + } + } + #[test] fn iou_is_symmetric_and_bounded( ax in -100.0f32..100.0, diff --git a/docs/site/algorithms.html b/docs/site/algorithms.html index e31c50f..6d212e6 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 visionConnected components, contours, polygon approximation, mask RLE, depth/confidence/flow/point maps, NMS/Soft-NMSspatialrust-vision · dense, detectionCPU + Dense visionExact Euclidean distance transform, connected components, contours, polygon approximation, mask RLE, depth/confidence/flow/point maps, NMS/Soft-NMSspatialrust-vision · dense, detectionCPU 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_euclidean_distance_transform.md b/notes/2026-07-15_exact_euclidean_distance_transform.md new file mode 100644 index 0000000..66407b7 --- /dev/null +++ b/notes/2026-07-15_exact_euclidean_distance_transform.md @@ -0,0 +1,64 @@ +# Exact Euclidean distance transform + +This slice adds a safe exact Euclidean distance transform (EDT) to +`spatialrust-vision`'s `dense` feature. For every nonzero binary-mask pixel, it +returns the physical L2 distance to the nearest zero pixel. Unit spacing and +positive finite anisotropic `(x, y)` spacing are supported. + +## Selection and algorithm + +EDT was selected after reviewing the existing catalog because it was absent and +is reusable for mask interior thickness, contour distance, watershed seeds, +feather blending, and signed-distance pipelines. The implementation follows +Felzenszwalb and Huttenlocher's separable lower-envelope transform: a 2D EDT is +computed as two 1D squared-distance transforms and one final square root. Its +time complexity is `O(width * height)` and its working memory is linear in the +image plus the longest axis. + +Primary references: + +- Pedro F. Felzenszwalb and Daniel P. Huttenlocher, "Distance Transforms of + Sampled Functions," Theory of Computing 8 (2012), + . +- OpenCV `distanceTransform` reference semantics and `DIST_MASK_PRECISE`, + . + +The public contract deliberately rejects a non-empty all-foreground mask, +because no finite nearest-background distance exists. Empty masks produce an +empty image. Python accepts conventional `0/255` masks by treating every +nonzero value as foreground; the Rust entry point uses validated `BinaryMask`. + +## Evidence + +- known-grid and anisotropic-spacing unit tests; +- generated-mask property comparison against exhaustive nearest-zero search; +- Python dtype/shape/value coverage; +- OpenCV exact-L2 comparison with maximum error gate `1e-5`; +- VGA/1080p/4K native Criterion and paired Python/OpenCV performance cases. + +### Local Windows receipt + +The release build ran on Windows 11 with CPython 3.12.10, OpenCV 4.10.0, +12 OpenCV threads, and OpenCL disabled. The canonical threshold-derived masks +matched OpenCV exactly at every profile. A separate irregular-mask correctness +case had maximum absolute error `9.536743e-7`. + +| Profile | Native Criterion estimate | Python SpatialRust median | OpenCV median | Outcome | +| --- | ---: | ---: | ---: | --- | +| VGA | 12.585 ms / 24.41 MPix/s | 19.795 ms | 1.867 ms | OpenCV 10.60x | +| 1080p | 100.98 ms / 20.53 MPix/s | 157.172 ms | 12.443 ms | OpenCV 12.63x | +| 4K | 451.63 ms / 18.37 MPix/s | 640.121 ms | 51.813 ms | OpenCV 12.35x | + +These results establish an honest scalar baseline, not a superiority claim. +OpenCV's tuned parallel implementation is faster on this host. Future work can +add caller-owned scratch and bounded row/column parallelism without changing +the exact public contract. + +Run from `C:\Users\rsasa\Workspace\SpatialRust`: + +```powershell +& "$HOME\.cargo\bin\cargo.exe" test -p spatialrust-vision --features full +& "$HOME\.cargo\bin\cargo.exe" bench -p spatialrust-vision --features dense --bench dense +.venv\Scripts\python.exe bench\opencv_vision_comparison\run.py +.venv\Scripts\python.exe bench\opencv_vision_comparison\performance.py +```