From b32d6384102fe4da1f15b79dc7291549f42dcc76 Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Thu, 16 Jul 2026 09:48:48 +0900 Subject: [PATCH] Add direct fused Sobel paths --- CHANGELOG.md | 9 + README.md | 16 +- .../README.md | 17 + .../performance.py | 199 +++++ crates/spatialrust-py/spatialrust.pyi | 16 +- crates/spatialrust-py/src/lib.rs | 179 ++++- crates/spatialrust-py/tests/test_bindings.py | 36 + crates/spatialrust-vision/benches/filter.rs | 22 +- .../spatialrust-vision/src/advanced_filter.rs | 700 +++++++++++++++++- docs/ROADMAP.md | 1 + docs/site/algorithms.html | 2 +- docs/site/filtering.html | 4 +- notes/2026-07-16_direct_sobel_threshold.md | 73 ++ 13 files changed, 1248 insertions(+), 26 deletions(-) create mode 100644 bench/opencv_sobel_threshold_comparison/README.md create mode 100644 bench/opencv_sobel_threshold_comparison/performance.py create mode 100644 notes/2026-07-16_direct_sobel_threshold.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 493dcb0..978b2d3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,15 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Added +- **Direct and fused 3×3 Sobel (Epic 116E)**: grayscale `u8` first derivatives + now use bounded parallel three-row `i16` rings instead of the generic + full-image `f64` intermediate. Added Rust/Python caller-output APIs for exact + `f32` derivatives, saturated absolute `u8` responses, and fused binary edge + masks. Packed NumPy inputs are borrowed without copying. Standalone Sobel + beats OpenCV 1.88× at 1080p and 2.03× at 4K; fused masks win 3.81×–6.64× + allocated and 2.95×–8.68× with caller-owned output across 300 bit-exact + randomized cases. + - **Allocation-light Canny (Epic 118A/118C)**: `canny()` no longer materializes public gradient, magnitude, and suppression images only to discard them. Added safe strided `canny_into`, reusable `CannyWorkspace`, large-image diff --git a/README.md b/README.md index 03c091e..4368608 100644 --- a/README.md +++ b/README.md @@ -161,7 +161,9 @@ ratio; these are machine-specific measurements, not universal guarantees. | Fused 2× resize → gray, allocate[^fused-gray-2026] | — | **SpatialRust 1.12×** | OpenCV 1.01× | | Fused 2× resize → gray, reuse[^fused-gray-2026] | — | OpenCV 1.90× | OpenCV 1.58× | | Gaussian blur 5×5[^gaussian-2026] | OpenCV 139.02× | OpenCV 3.10× | OpenCV 2.93× | -| Sobel X 3×3 | OpenCV 14.38× | OpenCV 20.31× | OpenCV 23.30× | +| Sobel X 3×3, allocate[^sobel-direct-2026] | OpenCV 1.07× | **SpatialRust 1.88×** | **SpatialRust 2.03×** | +| Fused abs(Sobel X) → binary mask, allocate[^sobel-direct-2026] | **SpatialRust 3.81×** | **SpatialRust 4.87×** | **SpatialRust 6.64×** | +| Fused abs(Sobel X) → binary mask, reuse[^sobel-direct-2026] | **SpatialRust 2.95×** | **SpatialRust 6.63×** | **SpatialRust 8.68×** | | Morphology open 5×5, allocate[^morphology-2026] | OpenCV 60.96× | OpenCV 13.34× | OpenCV 15.27× | | Morphology open 5×5, reuse[^morphology-2026] | OpenCV 60.32× | OpenCV 16.25× | OpenCV 17.78× | | Morphology open 511×511, allocate[^morphology-2026] | OpenCV 2.10× | **SpatialRust 2.61×** | **SpatialRust 2.40×** | @@ -241,6 +243,18 @@ OpenCV also remains faster for standalone `spatialGradient`. See the [focused harness](bench/opencv_sobel_l1_comparison/) and [dated receipt](notes/2026-07-16_paired_sobel_l1_acceleration.md). +[^sobel-direct-2026]: The grayscale `u8` 3×3 first-derivative path replaces + the generic full-image `f64` intermediate with parallel three-row `i16` + rings, writes `f32` directly, and borrows packed NumPy input without copying. + Against OpenCV 4.13, standalone allocation measured 1.134 ms versus 2.137 ms + at 1080p and 3.737 ms versus 7.582 ms at 4K, reversing the former + 20.31×–23.30× deficits while retaining max error zero. VGA remains a narrow + OpenCV win. `sobel_threshold_3x3_u8` additionally fuses signed Sobel, + absolute saturation, and binary threshold; it wins 3.81×–6.64× allocated and + 2.95×–8.68× with caller-owned output. Three hundred randomized X/Y cases are + bit-exact. See the + [focused harness](bench/opencv_sobel_threshold_comparison/). + [^morphology-2026]: Rectangular morphology was remeasured separately with OpenCV 4.13, OpenCL off, with both allocated and caller-owned-output Python API timing scopes. `MorphologyWorkspace` retains all full-image and diff --git a/bench/opencv_sobel_threshold_comparison/README.md b/bench/opencv_sobel_threshold_comparison/README.md new file mode 100644 index 0000000..51793bb --- /dev/null +++ b/bench/opencv_sobel_threshold_comparison/README.md @@ -0,0 +1,17 @@ +# OpenCV fused Sobel threshold comparison + +This harness compares an exact binary edge mask built from a first-order 3x3 +Sobel response. OpenCV uses `Sobel(CV_16S)`, `convertScaleAbs`, then +`threshold(THRESH_BINARY)`. SpatialRust fuses the same steps into one +three-row-ring operation and one `uint8` output. + +```powershell +python bench/opencv_sobel_threshold_comparison/performance.py ` + --output target/opencv-sobel-threshold-performance.json +``` + +OpenCL is disabled, inputs are seeded packed `uint8`, and allocate/reuse calls +are paired and interleaved. Timings are gated by exact pixels for both X and Y +derivatives across 300 randomized cases. Packed NumPy input is borrowed without +a copy; non-contiguous input is explicitly packed. Results are workload- and +host-specific. diff --git a/bench/opencv_sobel_threshold_comparison/performance.py b/bench/opencv_sobel_threshold_comparison/performance.py new file mode 100644 index 0000000..29c509e --- /dev/null +++ b/bench/opencv_sobel_threshold_comparison/performance.py @@ -0,0 +1,199 @@ +"""Reproducible fused Sobel-to-binary-mask comparison with OpenCV.""" + +from __future__ import annotations + +import argparse +import os +import sys +from pathlib import Path + +import cv2 +import numpy as np +import spatialrust as sr + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from opencv_comparison.report import emit_report, environment, make_report, timed_pair + + +PROFILES = { + "vga": (640, 480, 24), + "1080p": (1920, 1080, 16), + "4k": (3840, 2160, 10), +} +THRESHOLD = 96 + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser() + parser.add_argument("--output", type=Path) + parser.add_argument("--profiles", default=",".join(PROFILES)) + parser.add_argument("--warmup", type=int, default=6) + return parser.parse_args() + + +def opencv_allocate(image: np.ndarray, dx: int, dy: int) -> np.ndarray: + signed = cv2.Sobel( + image, + cv2.CV_16S, + dx, + dy, + ksize=3, + borderType=cv2.BORDER_REFLECT_101, + ) + absolute = cv2.convertScaleAbs(signed) + return cv2.threshold(absolute, THRESHOLD, 255, cv2.THRESH_BINARY)[1] + + +def validate_randomized_cases() -> int: + rng = np.random.default_rng(119) + checked = 0 + for case in range(300): + height = int(rng.integers(1, 120)) + width = int(rng.integers(1, 160)) + image = rng.integers(0, 256, (height, width), dtype=np.uint8) + if case % 3 == 0: + image = image[:, ::-1] + packed = np.ascontiguousarray(image) + dx, dy = ((1, 0), (0, 1))[case & 1] + expected = opencv_allocate(packed, dx, dy) + actual = sr.sobel_threshold_image(image, dx, dy, THRESHOLD) + if not np.array_equal(actual, expected): + raise AssertionError(f"random case {case} is not bit-exact") + checked += 1 + return checked + + +def main() -> None: + args = parse_args() + profiles = [value.strip() for value in args.profiles.split(",") if value.strip()] + unknown = sorted(set(profiles) - PROFILES.keys()) + if unknown: + raise ValueError(f"unknown profiles: {', '.join(unknown)}") + if hasattr(cv2, "ocl"): + cv2.ocl.setUseOpenCL(False) + cv2.setNumThreads(os.cpu_count() or 1) + + randomized_cases = validate_randomized_cases() + rng = np.random.default_rng(20_260_716) + results: dict[str, object] = {} + for profile in profiles: + width, height, repeats = PROFILES[profile] + image = rng.integers(0, 256, (height, width), dtype=np.uint8) + signed = np.empty((height, width), dtype=np.int16) + absolute = np.empty((height, width), dtype=np.uint8) + opencv_out = np.empty((height, width), dtype=np.uint8) + spatialrust_out = np.empty((height, width), dtype=np.uint8) + + def cv_allocate() -> np.ndarray: + return opencv_allocate(image, 1, 0) + + def sr_allocate() -> np.ndarray: + return sr.sobel_threshold_image(image, 1, 0, THRESHOLD) + + def cv_reuse() -> np.ndarray: + cv2.Sobel( + image, + cv2.CV_16S, + 1, + 0, + signed, + 3, + 1.0, + 0.0, + cv2.BORDER_REFLECT_101, + ) + cv2.convertScaleAbs(signed, absolute) + return cv2.threshold( + absolute, THRESHOLD, 255, cv2.THRESH_BINARY, opencv_out + )[1] + + def sr_reuse() -> np.ndarray: + return sr.sobel_threshold_image( + image, 1, 0, THRESHOLD, out=spatialrust_out + ) + + expected = cv_allocate() + if not np.array_equal(sr_allocate(), expected): + raise AssertionError(f"{profile} allocated output is not bit-exact") + if cv_reuse() is not opencv_out or sr_reuse() is not spatialrust_out: + raise AssertionError(f"{profile} caller-owned output identity failed") + if not np.array_equal(opencv_out, expected) or not np.array_equal( + spatialrust_out, expected + ): + raise AssertionError(f"{profile} reused output is not bit-exact") + + _, _, cv_timing, sr_timing = timed_pair( + cv_allocate, + sr_allocate, + warmup=args.warmup, + repeats=repeats, + seed=119, + min_sample_time_ms=20.0, + ) + _, _, cv_reuse_timing, sr_reuse_timing = timed_pair( + cv_reuse, + sr_reuse, + warmup=args.warmup, + repeats=repeats, + seed=2119, + min_sample_time_ms=20.0, + ) + cv_ms = float(cv_timing["median"]) + sr_ms = float(sr_timing["median"]) + cv_reuse_ms = float(cv_reuse_timing["median"]) + sr_reuse_ms = float(sr_reuse_timing["median"]) + results[profile] = { + "width": width, + "height": height, + "operation": "abs(Sobel X) > 96 binary mask", + "kernel_size": 3, + "border": "reflect101", + "exact": True, + "opencv_stages": ["Sobel CV_16S", "convertScaleAbs", "threshold"], + "spatialrust_stages": ["fused Sobel threshold"], + "opencv": cv_timing, + "spatialrust": sr_timing, + "spatialrust_speedup": cv_ms / sr_ms, + "faster_implementation": "spatialrust" if sr_ms < cv_ms else "opencv", + "opencv_reuse": cv_reuse_timing, + "spatialrust_reuse": sr_reuse_timing, + "spatialrust_reuse_speedup": cv_reuse_ms / sr_reuse_ms, + "faster_reuse_implementation": ( + "spatialrust" if sr_reuse_ms < cv_reuse_ms else "opencv" + ), + "spatialrust_reuse_vs_opencv_allocate_speedup": cv_ms / sr_reuse_ms, + "faster_spatialrust_reuse_vs_opencv_allocate": ( + "spatialrust" if sr_reuse_ms < cv_ms else "opencv" + ), + } + + receipt = environment( + opencv_version=cv2.__version__, spatialrust_version=sr.__version__ + ) + receipt["opencv_threads"] = cv2.getNumThreads() + receipt["opencv_opencl_enabled"] = bool( + hasattr(cv2, "ocl") and cv2.ocl.useOpenCL() + ) + report = make_report( + suite="opencv-fused-sobel-threshold-performance", + kind="performance", + status="pass", + environment_receipt=receipt, + results={ + "methodology": { + "timing_scope": "allocated and caller-owned-output Python API calls", + "paired_interleaved": True, + "minimum_sample_time_ms": 20.0, + "input": "seeded packed random uint8 grayscale", + "threshold": THRESHOLD, + "randomized_correctness_cases": randomized_cases, + "accuracy": "bit-exact binary mask for alternating X/Y derivatives", + }, + "profiles": results, + }, + ) + emit_report(report, args.output) + + +if __name__ == "__main__": + main() diff --git a/crates/spatialrust-py/spatialrust.pyi b/crates/spatialrust-py/spatialrust.pyi index 58382b3..e3e047e 100644 --- a/crates/spatialrust-py/spatialrust.pyi +++ b/crates/spatialrust-py/spatialrust.pyi @@ -30,7 +30,7 @@ __all__: list[str] = [ "rgbd_to_point_cloud", "depth_to_xyz", "calibrate_pinhole_camera", "calibrate_fisheye_angles", "dense_flow_image", "gray_world_white_balance_image", "stitch_panorama_pair", "filter2d_image", "gaussian_blur_image", - "median_blur_image", "bilateral_filter_image", "sobel_image", "spatial_gradient_image", + "median_blur_image", "bilateral_filter_image", "sobel_image", "sobel_abs_image", "sobel_threshold_image", "spatial_gradient_image", "sobel_l1_magnitude_image", "scharr_image", "laplacian_image", "pyr_down_image", "pyr_up_image", "MorphologyWorkspace", "morphology_image", "threshold_image", "otsu_threshold_image", "adaptive_threshold_image", @@ -183,7 +183,21 @@ def sobel_image( kernel_size: int = ..., scale: float = ..., delta: float = ..., + out: Optional[_F32Array] = ..., ) -> _F32Array: ... +def sobel_abs_image( + image: _U8Array, + dx: int, + dy: int, + out: Optional[_U8Array] = ..., +) -> _U8Array: ... +def sobel_threshold_image( + image: _U8Array, + dx: int, + dy: int, + threshold: int, + out: Optional[_U8Array] = ..., +) -> _U8Array: ... def spatial_gradient_image( image: _U8Array, out_dx: Optional[_I16Array] = ..., diff --git a/crates/spatialrust-py/src/lib.rs b/crates/spatialrust-py/src/lib.rs index 001d66b..0d9524b 100644 --- a/crates/spatialrust-py/src/lib.rs +++ b/crates/spatialrust-py/src/lib.rs @@ -95,8 +95,12 @@ use spatialrust::vision::{ resize_pack_chw_into as resize_pack_chw_into_op, resize_rgb_to_gray as resize_rgb_to_gray_op, resize_rgb_to_gray_into as resize_rgb_to_gray_into_op, rgb_to_gray as rgb_to_gray_op, rgb_to_gray_into as rgb_to_gray_into_op, rgb_to_hsv as rgb_to_hsv_op, scharr as scharr_op, - sobel as sobel_op, sobel_l1_magnitude_u8 as sobel_l1_magnitude_u8_op, - sobel_l1_magnitude_u8_into as sobel_l1_magnitude_u8_into_op, soft_nms as soft_nms_op, + sobel as sobel_op, sobel_3x3_u8 as sobel_3x3_u8_op, sobel_3x3_u8_into as sobel_3x3_u8_into_op, + sobel_abs_3x3_u8 as sobel_abs_3x3_u8_op, sobel_abs_3x3_u8_into as sobel_abs_3x3_u8_into_op, + sobel_l1_magnitude_u8 as sobel_l1_magnitude_u8_op, + sobel_l1_magnitude_u8_into as sobel_l1_magnitude_u8_into_op, + sobel_threshold_3x3_u8 as sobel_threshold_3x3_u8_op, + sobel_threshold_3x3_u8_into as sobel_threshold_3x3_u8_into_op, soft_nms as soft_nms_op, solve_pnp as solve_pnp_op, spatial_gradient_u8 as spatial_gradient_u8_op, spatial_gradient_u8_into as spatial_gradient_u8_into_op, stereo_block_match as stereo_block_match_op, stitch_panorama_pair as stitch_panorama_pair_op, @@ -2171,7 +2175,7 @@ fn bilateral_filter_image<'py>( /// Computes a signed float32 Sobel derivative from a grayscale uint8 image. #[pyfunction] -#[pyo3(signature = (image, dx, dy, kernel_size=3, scale=1.0, delta=0.0))] +#[pyo3(signature = (image, dx, dy, kernel_size=3, scale=1.0, delta=0.0, out=None))] fn sobel_image<'py>( py: Python<'py>, image: PyReadonlyArray2<'_, u8>, @@ -2180,9 +2184,150 @@ fn sobel_image<'py>( kernel_size: usize, scale: f64, delta: f64, + out: Option>>, ) -> PyResult>> { - let image = gray_u8_image_from_numpy(image)?; - let output = sobel_op(image.view(), dx, dy, kernel_size, scale, delta, BorderMode::Reflect101) + let mut packed = Vec::new(); + let image = gray_u8_image_view_from_numpy(&image, &mut packed)?; + if let Some(out) = out { + { + let mut out_rw = out + .try_readwrite() + .map_err(|_| PyValueError::new_err("out must not overlap the Sobel input"))?; + let mut out_array = out_rw.as_array_mut(); + if out_array.shape() != [image.height(), image.width()] { + return Err(PyValueError::new_err(format!( + "out shape must be ({}, {}), found {:?}", + image.height(), + image.width(), + out_array.shape() + ))); + } + let Some(out_slice) = out_array.as_slice_mut() else { + return Err(PyValueError::new_err( + "out must be a contiguous float32 array of shape (H, W)", + )); + }; + if kernel_size == 3 && dx + dy == 1 { + sobel_3x3_u8_into_op( + image, + dx, + dy, + scale, + delta, + BorderMode::Reflect101, + out_slice, + ) + .map_err(to_py_err)?; + } else { + let output = + sobel_op(image, dx, dy, kernel_size, scale, delta, BorderMode::Reflect101) + .map_err(to_py_err)?; + out_slice.copy_from_slice(output.as_slice()); + } + } + return Ok(out); + } + if kernel_size == 3 && dx + dy == 1 { + let output = sobel_3x3_u8_op(image, dx, dy, scale, delta, BorderMode::Reflect101) + .map_err(to_py_err)?; + let array = Array2::from_shape_vec((image.height(), image.width()), output.into_vec()) + .map_err(to_py_err)?; + return Ok(array.into_pyarray_bound(py)); + } + let output = sobel_op(image, dx, dy, kernel_size, scale, delta, BorderMode::Reflect101) + .map_err(to_py_err)?; + let array = Array2::from_shape_vec((image.height(), image.width()), output.into_vec()) + .map_err(to_py_err)?; + Ok(array.into_pyarray_bound(py)) +} + +/// Computes fused absolute 3x3 Sobel response as saturated uint8. +#[pyfunction] +#[pyo3(signature = (image, dx, dy, out=None))] +fn sobel_abs_image<'py>( + py: Python<'py>, + image: PyReadonlyArray2<'_, u8>, + dx: usize, + dy: usize, + out: Option>>, +) -> PyResult>> { + let mut packed = Vec::new(); + let image = gray_u8_image_view_from_numpy(&image, &mut packed)?; + if let Some(out) = out { + { + let mut out_rw = out + .try_readwrite() + .map_err(|_| PyValueError::new_err("out must not overlap the Sobel input"))?; + let mut out_array = out_rw.as_array_mut(); + if out_array.shape() != [image.height(), image.width()] { + return Err(PyValueError::new_err(format!( + "out shape must be ({}, {}), found {:?}", + image.height(), + image.width(), + out_array.shape() + ))); + } + let Some(out_slice) = out_array.as_slice_mut() else { + return Err(PyValueError::new_err( + "out must be a contiguous uint8 array of shape (H, W)", + )); + }; + sobel_abs_3x3_u8_into_op(image, dx, dy, BorderMode::Reflect101, out_slice) + .map_err(to_py_err)?; + } + return Ok(out); + } + let output = sobel_abs_3x3_u8_op(image, dx, dy, BorderMode::Reflect101).map_err(to_py_err)?; + let array = Array2::from_shape_vec((image.height(), image.width()), output.into_vec()) + .map_err(to_py_err)?; + Ok(array.into_pyarray_bound(py)) +} + +/// Computes a fused binary mask from absolute 3x3 Sobel response. +#[pyfunction] +#[pyo3(signature = (image, dx, dy, threshold, out=None))] +fn sobel_threshold_image<'py>( + py: Python<'py>, + image: PyReadonlyArray2<'_, u8>, + dx: usize, + dy: usize, + threshold: u8, + out: Option>>, +) -> PyResult>> { + let mut packed = Vec::new(); + let image = gray_u8_image_view_from_numpy(&image, &mut packed)?; + if let Some(out) = out { + { + let mut out_rw = out + .try_readwrite() + .map_err(|_| PyValueError::new_err("out must not overlap the Sobel input"))?; + let mut out_array = out_rw.as_array_mut(); + if out_array.shape() != [image.height(), image.width()] { + return Err(PyValueError::new_err(format!( + "out shape must be ({}, {}), found {:?}", + image.height(), + image.width(), + out_array.shape() + ))); + } + let Some(out_slice) = out_array.as_slice_mut() else { + return Err(PyValueError::new_err( + "out must be a contiguous uint8 array of shape (H, W)", + )); + }; + sobel_threshold_3x3_u8_into_op( + image, + dx, + dy, + threshold, + BorderMode::Reflect101, + out_slice, + ) + .map_err(to_py_err)?; + } + return Ok(out); + } + let output = sobel_threshold_3x3_u8_op(image, dx, dy, threshold, BorderMode::Reflect101) .map_err(to_py_err)?; let array = Array2::from_shape_vec((image.height(), image.width()), output.into_vec()) .map_err(to_py_err)?; @@ -2198,11 +2343,12 @@ fn spatial_gradient_image<'py>( out_dx: Option>>, out_dy: Option>>, ) -> PyResult<(Bound<'py, PyArray2>, Bound<'py, PyArray2>)> { - let image = gray_u8_image_from_numpy(image)?; + let mut packed = Vec::new(); + let image = gray_u8_image_view_from_numpy(&image, &mut packed)?; match (out_dx, out_dy) { (None, None) => { let (gradient_x, gradient_y) = - spatial_gradient_u8_op(image.view(), BorderMode::Reflect101).map_err(to_py_err)?; + spatial_gradient_u8_op(image, BorderMode::Reflect101).map_err(to_py_err)?; let gradient_x = Array2::from_shape_vec((image.height(), image.width()), gradient_x.into_vec()) .map_err(to_py_err)? @@ -2250,13 +2396,8 @@ fn spatial_gradient_image<'py>( "out_dy must be a contiguous int16 array of shape (H, W)", )); }; - spatial_gradient_u8_into_op( - image.view(), - BorderMode::Reflect101, - dx_slice, - dy_slice, - ) - .map_err(to_py_err)?; + spatial_gradient_u8_into_op(image, BorderMode::Reflect101, dx_slice, dy_slice) + .map_err(to_py_err)?; } Ok((out_dx, out_dy)) } @@ -2272,7 +2413,8 @@ fn sobel_l1_magnitude_image<'py>( image: PyReadonlyArray2<'_, u8>, out: Option>>, ) -> PyResult>> { - let image = gray_u8_image_from_numpy(image)?; + let mut packed = Vec::new(); + let image = gray_u8_image_view_from_numpy(&image, &mut packed)?; if let Some(out) = out { { let mut out_rw = out.try_readwrite().map_err(|_| { @@ -2292,13 +2434,12 @@ fn sobel_l1_magnitude_image<'py>( "out must be a contiguous int16 array of shape (H, W)", )); }; - sobel_l1_magnitude_u8_into_op(image.view(), BorderMode::Reflect101, out_slice) + sobel_l1_magnitude_u8_into_op(image, BorderMode::Reflect101, out_slice) .map_err(to_py_err)?; } return Ok(out); } - let magnitude = - sobel_l1_magnitude_u8_op(image.view(), BorderMode::Reflect101).map_err(to_py_err)?; + let magnitude = sobel_l1_magnitude_u8_op(image, BorderMode::Reflect101).map_err(to_py_err)?; let array = Array2::from_shape_vec((image.height(), image.width()), magnitude.into_vec()) .map_err(to_py_err)?; Ok(array.into_pyarray_bound(py)) @@ -4007,6 +4148,8 @@ fn spatialrust_module(m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(median_blur_image, m)?)?; m.add_function(wrap_pyfunction!(bilateral_filter_image, m)?)?; m.add_function(wrap_pyfunction!(sobel_image, m)?)?; + m.add_function(wrap_pyfunction!(sobel_abs_image, m)?)?; + m.add_function(wrap_pyfunction!(sobel_threshold_image, m)?)?; m.add_function(wrap_pyfunction!(spatial_gradient_image, m)?)?; m.add_function(wrap_pyfunction!(sobel_l1_magnitude_image, m)?)?; m.add_function(wrap_pyfunction!(scharr_image, m)?)?; diff --git a/crates/spatialrust-py/tests/test_bindings.py b/crates/spatialrust-py/tests/test_bindings.py index 5981daa..8a10a68 100644 --- a/crates/spatialrust-py/tests/test_bindings.py +++ b/crates/spatialrust-py/tests/test_bindings.py @@ -640,6 +640,42 @@ def test_spatial_gradient_matches_sobel_and_reuses_outputs(): np.testing.assert_array_equal(out_dy, expected_y) +def test_sobel_reuses_output_and_validates_it(): + gray = np.arange(13 * 17, dtype=np.uint8).reshape(13, 17)[:, ::-1] + expected = sr.sobel_image(gray, 1, 0) + output = np.empty(gray.shape, dtype=np.float32) + assert sr.sobel_image(gray, 1, 0, out=output) is output + np.testing.assert_array_equal(output, expected) + with pytest.raises(ValueError, match="out shape"): + sr.sobel_image(gray, 1, 0, out=np.empty((13, 16), dtype=np.float32)) + with pytest.raises(ValueError, match="contiguous"): + sr.sobel_image(gray, 1, 0, out=np.empty((13, 34), dtype=np.float32)[:, ::2]) + + +def test_absolute_sobel_matches_signed_pipeline_and_reuses_output(): + gray = np.arange(19 * 23, dtype=np.uint8).reshape(19, 23)[:, ::-1] + signed = sr.sobel_image(gray, 1, 0) + expected = np.minimum(np.abs(signed), 255).astype(np.uint8) + actual = sr.sobel_abs_image(gray, 1, 0) + np.testing.assert_array_equal(actual, expected) + output = np.empty(gray.shape, dtype=np.uint8) + assert sr.sobel_abs_image(gray, 1, 0, out=output) is output + np.testing.assert_array_equal(output, expected) + + +def test_thresholded_sobel_matches_absolute_pipeline_and_reuses_output(): + gray = np.arange(19 * 23, dtype=np.uint8).reshape(19, 23)[:, ::-1] + absolute = sr.sobel_abs_image(gray, 1, 0) + expected = np.where(absolute > 96, 255, 0).astype(np.uint8) + actual = sr.sobel_threshold_image(gray, 1, 0, 96) + np.testing.assert_array_equal(actual, expected) + output = np.empty(gray.shape, dtype=np.uint8) + assert sr.sobel_threshold_image(gray, 1, 0, 96, out=output) is output + np.testing.assert_array_equal(output, expected) + with pytest.raises(ValueError): + sr.sobel_threshold_image(gray, 1, 0, 96, out=gray) + + def test_spatial_gradient_rejects_partial_and_invalid_outputs(): gray = np.arange(7 * 9, dtype=np.uint8).reshape(7, 9) output = np.empty(gray.shape, dtype=np.int16) diff --git a/crates/spatialrust-vision/benches/filter.rs b/crates/spatialrust-vision/benches/filter.rs index 523ba2d..5a94300 100644 --- a/crates/spatialrust-vision/benches/filter.rs +++ b/crates/spatialrust-vision/benches/filter.rs @@ -2,8 +2,9 @@ use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criteri use spatialrust_image::Image; use spatialrust_vision::{ bilateral_filter, gaussian_blur, gaussian_blur_u8, gaussian_blur_u8_into, median_blur, - pyr_down, sobel, sobel_l1_magnitude_u8, sobel_l1_magnitude_u8_into, spatial_gradient_u8, - BorderMode, GaussianBlurU8Workspace, + pyr_down, sobel, sobel_3x3_u8, sobel_abs_3x3_u8, sobel_l1_magnitude_u8, + sobel_l1_magnitude_u8_into, sobel_threshold_3x3_u8, spatial_gradient_u8, BorderMode, + GaussianBlurU8Workspace, }; fn benchmark_gaussian(c: &mut Criterion) { @@ -115,6 +116,23 @@ fn benchmark_paired_sobel(c: &mut Criterion) { group.bench_function(BenchmarkId::new("xy_allocate", name), |b| { b.iter(|| spatial_gradient_u8(black_box(image.view()), BorderMode::Reflect101).unwrap()) }); + group.bench_function(BenchmarkId::new("x_f32_direct", name), |b| { + b.iter(|| { + sobel_3x3_u8(black_box(image.view()), 1, 0, 1.0, 0.0, BorderMode::Reflect101) + .unwrap() + }) + }); + group.bench_function(BenchmarkId::new("x_abs_u8", name), |b| { + b.iter(|| { + sobel_abs_3x3_u8(black_box(image.view()), 1, 0, BorderMode::Reflect101).unwrap() + }) + }); + group.bench_function(BenchmarkId::new("x_threshold_u8", name), |b| { + b.iter(|| { + sobel_threshold_3x3_u8(black_box(image.view()), 1, 0, 96, BorderMode::Reflect101) + .unwrap() + }) + }); group.bench_function(BenchmarkId::new("l1_allocate", name), |b| { b.iter(|| { sobel_l1_magnitude_u8(black_box(image.view()), BorderMode::Reflect101).unwrap() diff --git a/crates/spatialrust-vision/src/advanced_filter.rs b/crates/spatialrust-vision/src/advanced_filter.rs index e8df5ab..a5da599 100644 --- a/crates/spatialrust-vision/src/advanced_filter.rs +++ b/crates/spatialrust-vision/src/advanced_filter.rs @@ -117,6 +117,594 @@ pub fn sobel( separable_filter_f32(input, &kernel_x, &kernel_y, delta, border) } +/// Computes a first-order 3x3 Sobel derivative directly into signed `f32` output. +/// +/// This specialization accepts grayscale `u8` input, derivative order `(1, 0)` +/// or `(0, 1)`, and Replicate or Reflect101 borders. It avoids the generic +/// separable filter's full-image `f64` intermediate while preserving the same +/// coefficients, scale, delta, metadata, and strided-input behavior. +pub fn sobel_3x3_u8( + input: ImageView<'_, u8, 1>, + dx: usize, + dy: usize, + scale: f64, + delta: f64, + border: BorderMode, +) -> VisionResult> { + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("Sobel output size overflows".into()))?; + let mut output = vec![0.0; len]; + sobel_3x3_u8_into(input, dx, dy, scale, delta, border, &mut output)?; + Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?) +} + +/// Computes a first-order 3x3 Sobel derivative into caller-owned packed output. +pub fn sobel_3x3_u8_into( + input: ImageView<'_, u8, 1>, + dx: usize, + dy: usize, + scale: f64, + delta: f64, + border: BorderMode, + output: &mut [f32], +) -> VisionResult<()> { + validate_derivative(dx, dy, 3, scale, delta)?; + if dx + dy != 1 { + return Err(VisionError::InvalidParameter( + "specialized 3x3 Sobel requires derivative order (1, 0) or (0, 1)".into(), + )); + } + if !matches!(border, BorderMode::Replicate | BorderMode::Reflect101) { + return Err(VisionError::InvalidParameter( + "specialized 3x3 Sobel supports only Replicate and Reflect101 borders".into(), + )); + } + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("Sobel output size overflows".into()))?; + if output.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "Sobel output needs {len} elements, found {}", + output.len() + ))); + } + if len == 0 { + return Ok(()); + } + + if scale == 1.0 && delta == 0.0 { + sobel_3x3_identity_into(input, dx, border, output); + return Ok(()); + } + + let width = input.width(); + let height = input.height(); + let arch = Arch::new(); + if len >= 262_144 && height > 1 { + let workers = rayon::current_num_threads().min(height); + let rows_per_worker = height.div_ceil(workers); + output.par_chunks_mut(rows_per_worker * width).enumerate().for_each(|(chunk, output)| { + arch.dispatch(|| { + sobel_3x3_rows_dispatch( + input, + dx, + scale, + delta, + border, + chunk * rows_per_worker, + output, + ); + }); + }); + } else { + arch.dispatch(|| { + sobel_3x3_rows_dispatch(input, dx, scale, delta, border, 0, output); + }); + } + Ok(()) +} + +/// Computes an exact absolute first-order 3x3 Sobel response as saturated `u8`. +/// +/// This fuses signed derivative calculation, absolute value, and saturation, +/// avoiding the signed intermediate image required by a two-stage pipeline. +pub fn sobel_abs_3x3_u8( + input: ImageView<'_, u8, 1>, + dx: usize, + dy: usize, + border: BorderMode, +) -> VisionResult> { + let len = input.width().checked_mul(input.height()).ok_or_else(|| { + VisionError::InvalidDimensions("absolute Sobel output size overflows".into()) + })?; + let mut output = vec![0; len]; + sobel_abs_3x3_u8_into(input, dx, dy, border, &mut output)?; + Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?) +} + +/// Computes an exact absolute first-order 3x3 Sobel response into packed output. +pub fn sobel_abs_3x3_u8_into( + input: ImageView<'_, u8, 1>, + dx: usize, + dy: usize, + border: BorderMode, + output: &mut [u8], +) -> VisionResult<()> { + validate_derivative(dx, dy, 3, 1.0, 0.0)?; + if dx + dy != 1 { + return Err(VisionError::InvalidParameter( + "absolute 3x3 Sobel requires derivative order (1, 0) or (0, 1)".into(), + )); + } + if !matches!(border, BorderMode::Replicate | BorderMode::Reflect101) { + return Err(VisionError::InvalidParameter( + "absolute 3x3 Sobel supports only Replicate and Reflect101 borders".into(), + )); + } + let len = input.width().checked_mul(input.height()).ok_or_else(|| { + VisionError::InvalidDimensions("absolute Sobel output size overflows".into()) + })?; + if output.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "absolute Sobel output needs {len} elements, found {}", + output.len() + ))); + } + if len == 0 { + return Ok(()); + } + + sobel_abs_or_threshold_3x3_into(input, dx, border, None, output); + Ok(()) +} + +/// Computes a binary edge mask from absolute 3x3 Sobel response. +/// +/// Pixels whose saturated absolute response is strictly greater than +/// `threshold` become 255; all others become zero. This matches OpenCV's +/// `Sobel(CV_16S)` → `convertScaleAbs` → `THRESH_BINARY` pipeline. +pub fn sobel_threshold_3x3_u8( + input: ImageView<'_, u8, 1>, + dx: usize, + dy: usize, + threshold: u8, + border: BorderMode, +) -> VisionResult> { + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("Sobel mask size overflows".into()))?; + let mut output = vec![0; len]; + sobel_threshold_3x3_u8_into(input, dx, dy, threshold, border, &mut output)?; + Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?) +} + +/// Computes a binary absolute-Sobel edge mask into caller-owned packed output. +pub fn sobel_threshold_3x3_u8_into( + input: ImageView<'_, u8, 1>, + dx: usize, + dy: usize, + threshold: u8, + border: BorderMode, + output: &mut [u8], +) -> VisionResult<()> { + validate_derivative(dx, dy, 3, 1.0, 0.0)?; + if dx + dy != 1 { + return Err(VisionError::InvalidParameter( + "Sobel threshold requires derivative order (1, 0) or (0, 1)".into(), + )); + } + if !matches!(border, BorderMode::Replicate | BorderMode::Reflect101) { + return Err(VisionError::InvalidParameter( + "Sobel threshold supports only Replicate and Reflect101 borders".into(), + )); + } + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("Sobel mask size overflows".into()))?; + if output.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "Sobel mask needs {len} elements, found {}", + output.len() + ))); + } + if len == 0 { + return Ok(()); + } + sobel_abs_or_threshold_3x3_into(input, dx, border, Some(threshold), output); + Ok(()) +} + +fn sobel_abs_or_threshold_3x3_into( + input: ImageView<'_, u8, 1>, + dx: usize, + border: BorderMode, + threshold: Option, + output: &mut [u8], +) { + let width = input.width(); + let height = input.height(); + let len = output.len(); + let workers = + if len >= 262_144 && height > 1 { rayon::current_num_threads().min(height) } else { 1 }; + let rows_per_worker = height.div_ceil(workers); + let mut scratch = vec![0_i16; workers * 3 * width]; + scratch + .par_chunks_mut(3 * width) + .zip(output.par_chunks_mut(rows_per_worker * width)) + .enumerate() + .for_each(|(chunk, (scratch, output))| { + sobel_abs_3x3_stripe( + input, + dx, + border, + threshold, + chunk * rows_per_worker, + scratch, + output, + ); + }); +} + +fn sobel_abs_3x3_stripe( + input: ImageView<'_, u8, 1>, + dx: usize, + border: BorderMode, + threshold: Option, + start_y: usize, + scratch: &mut [i16], + output: &mut [u8], +) { + let width = input.width(); + let height = input.height(); + let mut slots = [0, 1, 2]; + let (top_y, bottom_y) = gradient_neighbors(start_y, height, border); + sobel_horizontal_row(input, top_y, dx, border, ring_row_mut(scratch, width, slots[0])); + sobel_horizontal_row(input, start_y, dx, border, ring_row_mut(scratch, width, slots[1])); + sobel_horizontal_row(input, bottom_y, dx, border, ring_row_mut(scratch, width, slots[2])); + + let stripe_rows = output.len() / width; + for (local_y, output) in output.chunks_mut(width).enumerate() { + let y = start_y + local_y; + let top = ring_row(scratch, width, slots[0]); + let middle = ring_row(scratch, width, slots[1]); + let bottom = ring_row(scratch, width, slots[2]); + if dx == 1 { + for (((output, &top), &middle), &bottom) in + output.iter_mut().zip(top).zip(middle).zip(bottom) + { + *output = threshold_sobel(top + 2 * middle + bottom, threshold); + } + } else { + for ((output, &top), &bottom) in output.iter_mut().zip(top).zip(bottom) { + *output = threshold_sobel(bottom - top, threshold); + } + } + if local_y + 1 < stripe_rows { + slots.rotate_left(1); + let (_, next_bottom_y) = gradient_neighbors(y + 1, height, border); + sobel_horizontal_row( + input, + next_bottom_y, + dx, + border, + ring_row_mut(scratch, width, slots[2]), + ); + } + } +} + +#[inline(always)] +fn saturating_abs_u8(value: i16) -> u8 { + value.unsigned_abs().min(u16::from(u8::MAX)) as u8 +} + +#[inline(always)] +fn threshold_sobel(value: i16, threshold: Option) -> u8 { + let value = saturating_abs_u8(value); + threshold.map_or(value, |threshold| u8::from(value > threshold) * u8::MAX) +} + +fn sobel_3x3_identity_into( + input: ImageView<'_, u8, 1>, + dx: usize, + border: BorderMode, + output: &mut [f32], +) { + let width = input.width(); + let height = input.height(); + let len = width * height; + let workers = + if len >= 262_144 && height > 1 { rayon::current_num_threads().min(height) } else { 1 }; + let rows_per_worker = height.div_ceil(workers); + let mut scratch = vec![0_i16; workers * 3 * width]; + scratch + .par_chunks_mut(3 * width) + .zip(output.par_chunks_mut(rows_per_worker * width)) + .enumerate() + .for_each(|(chunk, (scratch, output))| { + let start_y = chunk * rows_per_worker; + sobel_3x3_identity_stripe(input, dx, border, start_y, scratch, output); + }); +} + +fn sobel_3x3_identity_stripe( + input: ImageView<'_, u8, 1>, + dx: usize, + border: BorderMode, + start_y: usize, + scratch: &mut [i16], + output: &mut [f32], +) { + let width = input.width(); + let height = input.height(); + let mut slots = [0, 1, 2]; + let (top_y, bottom_y) = gradient_neighbors(start_y, height, border); + sobel_horizontal_row(input, top_y, dx, border, ring_row_mut(scratch, width, slots[0])); + sobel_horizontal_row(input, start_y, dx, border, ring_row_mut(scratch, width, slots[1])); + sobel_horizontal_row(input, bottom_y, dx, border, ring_row_mut(scratch, width, slots[2])); + + let stripe_rows = output.len() / width; + for (local_y, output) in output.chunks_mut(width).enumerate() { + let y = start_y + local_y; + let top = ring_row(scratch, width, slots[0]); + let middle = ring_row(scratch, width, slots[1]); + let bottom = ring_row(scratch, width, slots[2]); + if dx == 1 { + for (((output, &top), &middle), &bottom) in + output.iter_mut().zip(top).zip(middle).zip(bottom) + { + *output = (top + 2 * middle + bottom) as f32; + } + } else { + for ((output, &top), &bottom) in output.iter_mut().zip(top).zip(bottom) { + *output = (bottom - top) as f32; + } + } + + if local_y + 1 < stripe_rows { + slots.rotate_left(1); + let next_y = y + 1; + let (_, next_bottom_y) = gradient_neighbors(next_y, height, border); + sobel_horizontal_row( + input, + next_bottom_y, + dx, + border, + ring_row_mut(scratch, width, slots[2]), + ); + } + } +} + +fn sobel_horizontal_row( + input: ImageView<'_, u8, 1>, + y: usize, + dx: usize, + border: BorderMode, + output: &mut [i16], +) { + let width = input.width(); + let row = input.row(y).expect("Sobel horizontal row in bounds"); + if width == 1 { + output[0] = if dx == 1 { 0 } else { 4 * i16::from(row[0]) }; + return; + } + let (left, _) = gradient_neighbors(0, width, border); + output[0] = if dx == 1 { + i16::from(row[1]) - i16::from(row[left]) + } else { + i16::from(row[left]) + 2 * i16::from(row[0]) + i16::from(row[1]) + }; + if dx == 1 { + for (output, (&left, &right)) in + output[1..width - 1].iter_mut().zip(row[..width - 2].iter().zip(&row[2..])) + { + *output = i16::from(right) - i16::from(left); + } + } else { + for (output, ((&left, ¢er), &right)) in output[1..width - 1] + .iter_mut() + .zip(row[..width - 2].iter().zip(&row[1..width - 1]).zip(&row[2..])) + { + *output = i16::from(left) + 2 * i16::from(center) + i16::from(right); + } + } + let x = width - 1; + let (_, right) = gradient_neighbors(x, width, border); + output[x] = if dx == 1 { + i16::from(row[right]) - i16::from(row[x - 1]) + } else { + i16::from(row[x - 1]) + 2 * i16::from(row[x]) + i16::from(row[right]) + }; +} + +fn ring_row(scratch: &[i16], width: usize, slot: usize) -> &[i16] { + &scratch[slot * width..(slot + 1) * width] +} + +fn ring_row_mut(scratch: &mut [i16], width: usize, slot: usize) -> &mut [i16] { + &mut scratch[slot * width..(slot + 1) * width] +} + +#[inline] +fn sobel_3x3_rows_dispatch( + input: ImageView<'_, u8, 1>, + dx: usize, + scale: f64, + delta: f64, + border: BorderMode, + start_y: usize, + output: &mut [f32], +) { + if scale == 1.0 && delta == 0.0 { + if dx == 1 { + sobel_x_3x3_rows(input, border, start_y, output); + } else { + sobel_y_3x3_rows(input, border, start_y, output); + } + } else { + sobel_3x3_rows(input, dx, scale, delta, border, start_y, output); + } +} + +fn sobel_x_3x3_rows( + input: ImageView<'_, u8, 1>, + border: BorderMode, + start_y: usize, + output: &mut [f32], +) { + let width = input.width(); + let height = input.height(); + for (local_y, output) in output.chunks_mut(width).enumerate() { + let y = start_y + local_y; + let (top_y, bottom_y) = gradient_neighbors(y, height, border); + let top = input.row(top_y).expect("Sobel X row in bounds"); + let middle = input.row(y).expect("Sobel X row in bounds"); + let bottom = input.row(bottom_y).expect("Sobel X row in bounds"); + if width == 1 { + output[0] = 0.0; + continue; + } + let (left, _) = gradient_neighbors(0, width, border); + output[0] = sobel_x_value(top, middle, bottom, left, 1) as f32; + for ( + ((output, (top_left, top_right)), (middle_left, middle_right)), + (bottom_left, bottom_right), + ) in output[1..width - 1] + .iter_mut() + .zip(top[..width - 2].iter().zip(&top[2..])) + .zip(middle[..width - 2].iter().zip(&middle[2..])) + .zip(bottom[..width - 2].iter().zip(&bottom[2..])) + { + *output = (i16::from(*top_right) - i16::from(*top_left) + + 2 * (i16::from(*middle_right) - i16::from(*middle_left)) + + i16::from(*bottom_right) + - i16::from(*bottom_left)) as f32; + } + let x = width - 1; + let (_, right) = gradient_neighbors(x, width, border); + output[x] = sobel_x_value(top, middle, bottom, x - 1, right) as f32; + } +} + +fn sobel_y_3x3_rows( + input: ImageView<'_, u8, 1>, + border: BorderMode, + start_y: usize, + output: &mut [f32], +) { + let width = input.width(); + let height = input.height(); + for (local_y, output) in output.chunks_mut(width).enumerate() { + let y = start_y + local_y; + let (top_y, bottom_y) = gradient_neighbors(y, height, border); + let top = input.row(top_y).expect("Sobel Y row in bounds"); + let bottom = input.row(bottom_y).expect("Sobel Y row in bounds"); + if width == 1 { + output[0] = (4 * (i16::from(bottom[0]) - i16::from(top[0]))) as f32; + continue; + } + let (left, _) = gradient_neighbors(0, width, border); + output[0] = sobel_y_value(top, bottom, left, 0, 1) as f32; + for ( + (output, ((top_left, top_center), top_right)), + ((bottom_left, bottom_center), bottom_right), + ) in output[1..width - 1] + .iter_mut() + .zip(top[..width - 2].iter().zip(&top[1..width - 1]).zip(&top[2..])) + .zip(bottom[..width - 2].iter().zip(&bottom[1..width - 1]).zip(&bottom[2..])) + { + *output = (i16::from(*bottom_left) - i16::from(*top_left) + + 2 * (i16::from(*bottom_center) - i16::from(*top_center)) + + i16::from(*bottom_right) + - i16::from(*top_right)) as f32; + } + let x = width - 1; + let (_, right) = gradient_neighbors(x, width, border); + output[x] = sobel_y_value(top, bottom, x - 1, x, right) as f32; + } +} + +#[inline(always)] +fn sobel_x_value(top: &[u8], middle: &[u8], bottom: &[u8], left: usize, right: usize) -> i16 { + i16::from(top[right]) - i16::from(top[left]) + + 2 * (i16::from(middle[right]) - i16::from(middle[left])) + + i16::from(bottom[right]) + - i16::from(bottom[left]) +} + +#[inline(always)] +fn sobel_y_value(top: &[u8], bottom: &[u8], left: usize, center: usize, right: usize) -> i16 { + i16::from(bottom[left]) - i16::from(top[left]) + + 2 * (i16::from(bottom[center]) - i16::from(top[center])) + + i16::from(bottom[right]) + - i16::from(top[right]) +} + +fn sobel_3x3_rows( + input: ImageView<'_, u8, 1>, + dx: usize, + scale: f64, + delta: f64, + border: BorderMode, + start_y: usize, + output: &mut [f32], +) { + let width = input.width(); + let height = input.height(); + for (local_y, output) in output.chunks_mut(width).enumerate() { + let y = start_y + local_y; + let (top_y, bottom_y) = gradient_neighbors(y, height, border); + let top = input.row(top_y).expect("Sobel row in bounds"); + let middle = input.row(y).expect("Sobel row in bounds"); + let bottom = input.row(bottom_y).expect("Sobel row in bounds"); + if width == 1 { + output[0] = scaled_sobel_3x3_pixel(top, middle, bottom, 0, 0, 0, dx, scale, delta); + continue; + } + let (left, _) = gradient_neighbors(0, width, border); + output[0] = scaled_sobel_3x3_pixel(top, middle, bottom, left, 0, 1, dx, scale, delta); + for (x, value) in output.iter_mut().enumerate().take(width - 1).skip(1) { + *value = scaled_sobel_3x3_pixel(top, middle, bottom, x - 1, x, x + 1, dx, scale, delta); + } + let x = width - 1; + let (_, right) = gradient_neighbors(x, width, border); + output[x] = scaled_sobel_3x3_pixel(top, middle, bottom, x - 1, x, right, dx, scale, delta); + } +} + +#[inline(always)] +#[allow(clippy::too_many_arguments)] +fn scaled_sobel_3x3_pixel( + top: &[u8], + middle: &[u8], + bottom: &[u8], + left: usize, + center: usize, + right: usize, + dx: usize, + scale: f64, + delta: f64, +) -> f32 { + let value = if dx == 1 { + i16::from(top[right]) + 2 * i16::from(middle[right]) + i16::from(bottom[right]) + - i16::from(top[left]) + - 2 * i16::from(middle[left]) + - i16::from(bottom[left]) + } else { + i16::from(bottom[left]) + 2 * i16::from(bottom[center]) + i16::from(bottom[right]) + - i16::from(top[left]) + - 2 * i16::from(top[center]) + - i16::from(top[right]) + }; + f64::from(value).mul_add(scale, delta) as f32 +} + /// Computes exact 3×3 Sobel X/Y gradients together for grayscale `u8` input. /// /// This matches OpenCV `spatialGradient`: outputs are signed `i16`, the two @@ -633,8 +1221,9 @@ fn convolve_coefficients(left: &[f64], right: &[f64]) -> Vec { mod tests { use super::{ bilateral_filter, build_gaussian_pyramid, laplacian, median_blur, pyr_down, scharr, sobel, - sobel_l1_magnitude_u8, sobel_l1_magnitude_u8_into, spatial_gradient_u8, - spatial_gradient_u8_into, + sobel_3x3_u8, sobel_3x3_u8_into, sobel_abs_3x3_u8, sobel_abs_3x3_u8_into, + sobel_l1_magnitude_u8, sobel_l1_magnitude_u8_into, sobel_threshold_3x3_u8, + sobel_threshold_3x3_u8_into, spatial_gradient_u8, spatial_gradient_u8_into, }; use crate::BorderMode; use spatialrust_image::{Image, ImageRegion}; @@ -695,6 +1284,113 @@ mod tests { } } + #[test] + fn specialized_sobel_matches_generic_for_strided_input() { + let parent = Image::::try_new( + 11, + 8, + (0..88).map(|index| ((index * 41 + 19) % 256) as u8).collect(), + ) + .unwrap(); + let input = parent.view().subview(ImageRegion::new(1, 1, 9, 6)).unwrap(); + for border in [BorderMode::Replicate, BorderMode::Reflect101] { + for (dx, dy) in [(1, 0), (0, 1)] { + let expected = sobel(input, dx, dy, 3, 0.75, -2.5, border).unwrap(); + let actual = sobel_3x3_u8(input, dx, dy, 0.75, -2.5, border).unwrap(); + assert_eq!(actual, expected); + assert_eq!(actual.metadata(), input.metadata()); + } + } + } + + #[test] + fn specialized_sobel_into_validates_contract() { + let image = Image::::try_new(3, 2, vec![0, 1, 2, 3, 4, 5]).unwrap(); + let mut output = vec![0.0; 6]; + sobel_3x3_u8_into(image.view(), 1, 0, 1.0, 0.0, BorderMode::Reflect101, &mut output) + .unwrap(); + assert!(sobel_3x3_u8_into( + image.view(), + 1, + 0, + 1.0, + 0.0, + BorderMode::Reflect101, + &mut output[..5], + ) + .is_err()); + assert!(sobel_3x3_u8_into( + image.view(), + 1, + 1, + 1.0, + 0.0, + BorderMode::Reflect101, + &mut output, + ) + .is_err()); + assert!(sobel_3x3_u8_into(image.view(), 1, 0, 1.0, 0.0, BorderMode::Wrap, &mut output,) + .is_err()); + } + + #[test] + fn absolute_sobel_matches_saturated_specialized_derivative() { + let parent = Image::::try_new( + 12, + 9, + (0..108).map(|index| ((index * 67 + 13) % 256) as u8).collect(), + ) + .unwrap(); + let input = parent.view().subview(ImageRegion::new(1, 1, 10, 7)).unwrap(); + for border in [BorderMode::Replicate, BorderMode::Reflect101] { + for (dx, dy) in [(1, 0), (0, 1)] { + let derivative = sobel_3x3_u8(input, dx, dy, 1.0, 0.0, border).unwrap(); + let expected = derivative + .as_slice() + .iter() + .map(|value| value.abs().min(255.0) as u8) + .collect::>(); + let actual = sobel_abs_3x3_u8(input, dx, dy, border).unwrap(); + assert_eq!(actual.as_slice(), expected); + let mut reused = vec![0; input.width() * input.height()]; + sobel_abs_3x3_u8_into(input, dx, dy, border, &mut reused).unwrap(); + assert_eq!(reused, expected); + } + } + } + + #[test] + fn thresholded_sobel_matches_absolute_response() { + let image = Image::::try_new( + 17, + 13, + (0..221).map(|index| ((index * 29 + 7) % 256) as u8).collect(), + ) + .unwrap(); + for (dx, dy) in [(1, 0), (0, 1)] { + let absolute = sobel_abs_3x3_u8(image.view(), dx, dy, BorderMode::Reflect101).unwrap(); + let expected = absolute + .as_slice() + .iter() + .map(|&value| u8::from(value > 96) * u8::MAX) + .collect::>(); + let actual = + sobel_threshold_3x3_u8(image.view(), dx, dy, 96, BorderMode::Reflect101).unwrap(); + assert_eq!(actual.as_slice(), expected); + let mut reused = vec![0; expected.len()]; + sobel_threshold_3x3_u8_into( + image.view(), + dx, + dy, + 96, + BorderMode::Reflect101, + &mut reused, + ) + .unwrap(); + assert_eq!(reused, expected); + } + } + #[test] fn paired_spatial_gradient_into_validates_outputs_and_borders() { let image = Image::::try_new(3, 2, vec![0, 1, 2, 3, 4, 5]).unwrap(); diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 8524d46..8a9baae 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -640,6 +640,7 @@ to one implicitly, and GPU receipts must retain named upload/readback stages. | 116B | Complete | Split border handling from contiguous interior loops | five Gaussian border modes plus strided-view tests | | 116C | In progress | Specialized 3x3/5x5/7x7 Gaussian and paired Sobel X/Y passes | exact paired Sobel and accelerated Q8 3×3/5×5 Gaussian complete; high-precision 7×7 fallback remains | | 116D | Complete | Improve Gaussian by at least 10x and Sobel by at least 5x on one canonical large profile | 5×5 Gaussian 20.7× at 1080p and 26.7× at 4K; dated native timing receipt | +| 116E | Complete | Remove standalone 3×3 Sobel's generic `f64` intermediate and fuse common absolute-threshold consumers | exact direct/absolute/mask APIs; direct Sobel beats OpenCV 1.88×/2.03× at 1080p/4K; fused masks win 2.95×–8.68× | ### Epic 117 delivery slices diff --git a/docs/site/algorithms.html b/docs/site/algorithms.html index 2c98945..42bab1a 100644 --- a/docs/site/algorithms.html +++ b/docs/site/algorithms.html @@ -31,7 +31,7 @@

Algorithm catalog

PreprocessingCrop, pad, letterbox, normalize, interleaved-to-CHW, RGB/BGR swap, RGB↔gray/HSV, reusable outputs. Fused bilinear resize-to-gray avoids an intermediate RGB image and wins 1.12× at allocated 1080p→540p. Fused resize-normalize-CHW avoids the HWC intermediate and measured 2.02×–2.33× faster than OpenCV blobFromImage. Both retain exact SpatialRust unfused parity. CHW harness.spatialrust-vision · preprocessCPU Resize and warpNearest, bilinear, bicubic and area resize; reusable fixed-point RGB8 bilinear plans; exact row-parallel half-scale path; remap; affine and perspective warp. VGA caller-output half-scale measured 1.10× faster than OpenCV with exact pixels. Harness.spatialrust-vision · resize, warpCPU / GPU - Filtering2D correlation/convolution, separable/box/Gaussian, accelerated 3×3/5×5 u8 Gaussian with workspace reuse, median, bilateral, Sobel, exact paired gradients and fused L1 magnitude, Scharr, Laplacian, Gaussian pyramidsspatialrust-vision · imgproc-filtersafe CPU explicit GPU API + Filtering2D correlation/convolution, separable/box/Gaussian, accelerated 3×3/5×5 u8 Gaussian, median, bilateral, direct 3×3 Sobel, exact paired gradients, fused absolute/binary Sobel responses and L1 magnitude, Scharr, Laplacian, Gaussian pyramids. Exact direct Sobel beats OpenCV 1.88× at 1080p and 2.03× at 4K; fused threshold masks win 2.95×–8.68× depending on allocation mode. Harness.spatialrust-vision · imgproc-filtersafe CPU explicit GPU API MorphologyErode, dilate, open, close, gradient, top-hat, black-hat; separable sliding Rect fast path plus Cross/Ellipse/Diamond/custom fallbackspatialrust-vision · imgproc-morphologysafe CPU explicit wgpu API 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 diff --git a/docs/site/filtering.html b/docs/site/filtering.html index 3d8d135..2589802 100644 --- a/docs/site/filtering.html +++ b/docs/site/filtering.html @@ -13,7 +13,7 @@
imgproc-filter · safe CPU

Filtering and paired Sobel gradients

-

Generic correlation and separable filters remain available for arbitrary components and channels. Interleaved u8 additionally gets accelerated Gaussian kernels, explicit workspace reuse, exact paired 3×3 Sobel, and fused L1 magnitude.

+

Generic correlation and separable filters remain available for arbitrary components and channels. Interleaved u8 additionally gets accelerated Gaussian kernels, direct 3×3 Sobel, exact paired gradients, fused absolute/binary responses, and fused L1 magnitude.

Gaussian APIPaired-gradient APIFused L1 API
@@ -21,6 +21,7 @@

Filtering and paired Sobel gradients

Contracts

spatial_gradient_u8 matches OpenCV spatialGradient: grayscale u8 input, signed i16 X/Y outputs, 3×3 Sobel coefficients, and Replicate or Reflect101 borders. Its *_into form writes caller-owned packed slices.

sobel_l1_magnitude_u8 computes abs(Gx) + abs(Gy) directly as non-negative i16 values in [0, 2040]. It does not materialize public X/Y intermediates. Python exposes the same operations as spatial_gradient_image and sobel_l1_magnitude_image(..., out=).

+

sobel_3x3_u8 writes a first-order X or Y response directly to f32 from parallel three-row i16 rings. sobel_abs_3x3_u8 saturates the absolute response to u8, while sobel_threshold_3x3_u8 also produces a binary mask in the same pass. Python exposes caller-owned out= forms.

gaussian_blur_u8 accepts 3×3, 5×5, or 7×7 kernels and all five CPU border modes. The accelerated 3×3/5×5 path uses symmetric Q8 coefficients, branch-free unrolled interiors, bounded row parallelism, and a caller-owned GaussianBlurU8Workspace. The 7×7 contract currently uses the high-precision generic fallback.

Large packed images use bounded Rayon row stages and runtime-dispatched safe SIMD context. Strided inputs retain exact row semantics. CPU calls never upload to a device.

@@ -33,6 +34,7 @@

Measured boundary

2.42×

8K allocated fused L1 lead; 300 randomized cases are bit-exact.

OpenCV uses spatialGradient, two absdiff stages, and add. SpatialRust fuses that public workload into one traversal and one result allocation. Reuse is effectively tied at 1080p while OpenCV leads at 4K and 8K, and OpenCV remains faster for standalone paired gradients.

+

The direct CV_32F-compatible Sobel path combines three-row rings with zero-copy packed NumPy input. It is 1.88× faster than OpenCV at 1080p and 2.03× at 4K with exact values; VGA remains a narrow OpenCV win. The fused abs(Sobel X) > 96 mask wins 3.81×–6.64× allocated and 2.95×–8.68× with caller-owned output. Reproduce with the threshold harness.

Reproduce with the focused harness and inspect the dated receipt.

Gaussian progress

diff --git a/notes/2026-07-16_direct_sobel_threshold.md b/notes/2026-07-16_direct_sobel_threshold.md new file mode 100644 index 0000000..e403422 --- /dev/null +++ b/notes/2026-07-16_direct_sobel_threshold.md @@ -0,0 +1,73 @@ +# Epic 116E direct and fused Sobel receipt (2026-07-16) + +## Outcome + +The standard grayscale `u8`, first-order 3×3 Sobel path no longer enters the +generic separable engine or allocates a full `f64` intermediate. It uses +bounded parallel stripes, three reusable `i16` horizontal rows per worker, and +direct `f32` output. Additive fused APIs produce saturated absolute `u8` +responses or binary edge masks without public signed/absolute intermediates. + +## Public APIs + +Rust: + +- `sobel_3x3_u8` / `sobel_3x3_u8_into` +- `sobel_abs_3x3_u8` / `sobel_abs_3x3_u8_into` +- `sobel_threshold_3x3_u8` / `sobel_threshold_3x3_u8_into` + +Python: + +- `sobel_image(..., out=None)` dispatches first-order 3×3 calls to the direct path +- `sobel_abs_image(..., out=None)` +- `sobel_threshold_image(..., threshold, out=None)` + +Rust supports Replicate and Reflect101 borders. Python uses Reflect101. Generic +Sobel remains available for other sizes, orders, components, and channels. + +## OpenCV comparison + +Environment: Windows 11, Intel Family 6 Model 158, 6 cores / 12 logical CPUs, +CPython 3.12.10, OpenCV 4.13.0, OpenCL disabled, 12 OpenCV threads. Inputs are +seeded packed random grayscale `u8`; calls are paired/interleaved and adaptively +batched. + +Standalone `CV_32F` Sobel X allocation: + +| Profile | OpenCV | SpatialRust | Result | +| --- | ---: | ---: | ---: | +| VGA | 0.372 ms | 0.398 ms | OpenCV 1.07× | +| 1080p | 2.137 ms | 1.134 ms | **SpatialRust 1.88×** | +| 4K | 7.582 ms | 3.737 ms | **SpatialRust 2.03×** | + +This reverses historical 20.31× and 23.30× deficits at 1080p and 4K. Packed +NumPy input is borrowed; non-contiguous input is explicitly packed. All +canonical outputs have maximum absolute error zero. + +Fused `abs(Sobel X) > 96` binary mask: + +| Profile | OpenCV allocate | SpatialRust allocate | SpatialRust caller output | Caller output vs OpenCV allocate | +| --- | ---: | ---: | ---: | ---: | +| VGA | 0.304 ms | 0.080 ms | 0.069 ms | **SpatialRust 4.38×** | +| 1080p | 2.545 ms | 0.522 ms | 0.187 ms | **SpatialRust 13.58×** | +| 4K | 9.453 ms | 1.423 ms | 0.635 ms | **SpatialRust 14.90×** | + +OpenCV's allocated pipeline is `Sobel(CV_16S)` → `convertScaleAbs` → binary +`threshold`. Same-mode allocated SpatialRust wins 3.81×, 4.87×, and 6.64×; +same-mode caller-output SpatialRust wins 2.95×, 6.63×, and 8.68×. The larger +cross-mode ratios in the table explain the effect of one fused output versus +three allocated OpenCV stages; they are not substituted for the same-mode +claims. + +## Correctness and validation + +- direct path equals the generic SpatialRust Sobel for X/Y, scale/delta, both + supported borders, metadata, and strided Rust inputs +- Python `out=` identity plus shape/contiguity rejection +- 300 randomized OpenCV cases alternating X/Y, including non-contiguous inputs +- exact binary pixels for all timed profiles +- focused Rust tests, feature Clippy with warnings denied, Python crate check +- reproducible JSON: + `C:\Users\rsasa\Workspace\SpatialRust\target\opencv-sobel-threshold-performance.json` + and + `C:\Users\rsasa\Workspace\SpatialRust\target\opencv-comparison\vision-performance-sobel-direct.json`