rsasaki0109 · rsasaki0109 · Jul 15, 2026 · Jul 15, 2026
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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<Bound<'py, PyArray2<f32>>>,
+    workspace: Option<PyRefMut<'_, PyDistanceTransformWorkspace>>,
 ) -> PyResult<Bound<'py, PyArray2<f32>>> {
+    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::<Vec<_>>();
+                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::<PyImageMetadata>()?;
     m.add_class::<PyTensor>()?;
     m.add_class::<PyKeypoint2>()?;
+    m.add_class::<PyDistanceTransformWorkspace>()?;
     m.add_class::<PyOnnxRuntimeSession>()?;
     m.add_class::<PyDlpackTensorView>()?;
     m.add_class::<PyPointCloud>()?;

diff --git 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
@@ -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
@@ -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();
 }