diff --git a/Cargo.toml b/Cargo.toml index 78481b6..aff7c0a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -92,6 +92,7 @@ wgpu = "24" criterion = { version = "0.5", features = ["html_reports"] } proptest = "1" rayon = "1.12" +pulp = "0.18.22" image = { version = "0.24.9", default-features = false } kamadak-exif = "0.6.1" exr = { version = "=1.72.0", default-features = false } diff --git a/README.md b/README.md index 1434862..f939a26 100644 --- a/README.md +++ b/README.md @@ -173,6 +173,16 @@ 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 additive paired-gradient path keeps standalone Sobel compatibility while +also exposing exact fused 3×3 L1 magnitude (`abs(Gx) + abs(Gy)`). On a newer +OpenCV 4.13 receipt, the fused allocated Python call is **1.86× faster at +1080p, 2.19× at 4K, and 2.42× at 8K** because SpatialRust writes one result +instead of materializing paired gradients, two absolute-value images, and an +addition result. Caller-owned reuse ties at 1080p and favors OpenCV at 4K/8K; +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). + [^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_l1_comparison/README.md b/bench/opencv_sobel_l1_comparison/README.md new file mode 100644 index 0000000..aacfa68 --- /dev/null +++ b/bench/opencv_sobel_l1_comparison/README.md @@ -0,0 +1,19 @@ +# OpenCV fused Sobel L1 comparison + +This harness compares exact 3×3 grayscale Sobel L1 magnitude, +`abs(Gx) + abs(Gy)`, through public Python APIs. OpenCV uses its paired +`spatialGradient` primitive followed by two `absdiff` stages and `add`; +SpatialRust fuses the same integer operation into one source traversal and one +output write. + +```powershell +python bench/opencv_sobel_l1_comparison/performance.py ` + --output target/opencv-sobel-l1-performance.json +``` + +OpenCL is disabled, OpenCV receives the logical CPU count, input is seeded +packed random `uint8`, and allocate/reuse calls are paired and interleaved. +Every timing is gated by exact `int16` equality plus 300 randomized cases. The +result is a machine-specific fused-workload receipt, not a claim that the +standalone SpatialRust paired-gradient primitive beats OpenCV +`spatialGradient`. diff --git a/bench/opencv_sobel_l1_comparison/performance.py b/bench/opencv_sobel_l1_comparison/performance.py new file mode 100644 index 0000000..f924ded --- /dev/null +++ b/bench/opencv_sobel_l1_comparison/performance.py @@ -0,0 +1,194 @@ +"""Reproducible fused Sobel L1 magnitude 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 = { + "1080p": (1920, 1080, 20), + "4k": (3840, 2160, 14), + "8k": (7680, 4320, 8), +} + + +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_l1_allocate(image: np.ndarray, zero: np.ndarray) -> np.ndarray: + gradient_x, gradient_y = cv2.spatialGradient( + image, ksize=3, borderType=cv2.BORDER_REFLECT_101 + ) + absolute_x = cv2.absdiff(gradient_x, zero) + absolute_y = cv2.absdiff(gradient_y, zero) + return cv2.add(absolute_x, absolute_y) + + +def validate_randomized_cases() -> int: + rng = np.random.default_rng(116) + 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) + zero = np.zeros((height, width), dtype=np.int16) + expected = opencv_l1_allocate(packed, zero) + actual = sr.sobel_l1_magnitude_image(image) + 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) + zero = np.zeros((height, width), dtype=np.int16) + opencv_dx = np.empty((height, width), dtype=np.int16) + opencv_dy = np.empty((height, width), dtype=np.int16) + opencv_abs_x = np.empty((height, width), dtype=np.int16) + opencv_abs_y = np.empty((height, width), dtype=np.int16) + opencv_out = np.empty((height, width), dtype=np.int16) + spatialrust_out = np.empty((height, width), dtype=np.int16) + + def opencv_allocate() -> np.ndarray: + return opencv_l1_allocate(image, zero) + + def spatialrust_allocate() -> np.ndarray: + return sr.sobel_l1_magnitude_image(image) + + def opencv_reuse() -> np.ndarray: + cv2.spatialGradient( + image, + opencv_dx, + opencv_dy, + 3, + cv2.BORDER_REFLECT_101, + ) + cv2.absdiff(opencv_dx, zero, opencv_abs_x) + cv2.absdiff(opencv_dy, zero, opencv_abs_y) + return cv2.add(opencv_abs_x, opencv_abs_y, opencv_out) + + def spatialrust_reuse() -> np.ndarray: + return sr.sobel_l1_magnitude_image(image, spatialrust_out) + + expected = opencv_allocate() + actual = spatialrust_allocate() + if not np.array_equal(actual, expected): + raise AssertionError(f"{profile} allocated output is not bit-exact") + if opencv_reuse() is not opencv_out: + raise AssertionError("OpenCV did not return its caller-owned output") + if spatialrust_reuse() is not spatialrust_out: + raise AssertionError("SpatialRust did not return its caller-owned output") + 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") + + _, _, opencv_timing, spatialrust_timing = timed_pair( + opencv_allocate, + spatialrust_allocate, + warmup=args.warmup, + repeats=repeats, + seed=116, + min_sample_time_ms=20.0, + ) + _, _, opencv_reuse_timing, spatialrust_reuse_timing = timed_pair( + opencv_reuse, + spatialrust_reuse, + warmup=args.warmup, + repeats=repeats, + seed=2116, + min_sample_time_ms=20.0, + ) + opencv_ms = float(opencv_timing["median"]) + spatialrust_ms = float(spatialrust_timing["median"]) + opencv_reuse_ms = float(opencv_reuse_timing["median"]) + spatialrust_reuse_ms = float(spatialrust_reuse_timing["median"]) + results[profile] = { + "width": width, + "height": height, + "operation": "abs(Sobel X) + abs(Sobel Y)", + "kernel_size": 3, + "border": "reflect101", + "dtype": "int16", + "exact": True, + "opencv_stages": ["spatialGradient", "absdiff X", "absdiff Y", "add"], + "spatialrust_stages": ["fused Sobel L1"], + "opencv": opencv_timing, + "spatialrust": spatialrust_timing, + "spatialrust_speedup": opencv_ms / spatialrust_ms, + "faster_implementation": ( + "spatialrust" if spatialrust_ms < opencv_ms else "opencv" + ), + "opencv_reuse": opencv_reuse_timing, + "spatialrust_reuse": spatialrust_reuse_timing, + "spatialrust_reuse_speedup": opencv_reuse_ms / spatialrust_reuse_ms, + "faster_reuse_implementation": ( + "spatialrust" + if spatialrust_reuse_ms < opencv_reuse_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-l1-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", + "randomized_correctness_cases": randomized_cases, + "thread_policy": "logical CPU count for OpenCV; Rayon default for SpatialRust", + "accuracy": "bit-exact int16 L1 magnitude", + }, + "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 4672b12..24af65b 100644 --- a/crates/spatialrust-py/spatialrust.pyi +++ b/crates/spatialrust-py/spatialrust.pyi @@ -30,7 +30,8 @@ __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", "scharr_image", + "median_blur_image", "bilateral_filter_image", "sobel_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", "histogram_image", "equalize_histogram_image", "clahe_image", @@ -48,6 +49,7 @@ _F32Array = NDArray[np.float32] # positions, grids, range images, transforms _F64Array = NDArray[np.float64] _BoolArray = NDArray[np.bool_] _I32Array = NDArray[np.int32] # labels, edge_index +_I16Array = NDArray[np.int16] _U32Array = NDArray[np.uint32] _Vec3 = tuple[float, float, float] _U8Array = NDArray[np.uint8] @@ -181,6 +183,15 @@ def sobel_image( scale: float = ..., delta: float = ..., ) -> _F32Array: ... +def spatial_gradient_image( + image: _U8Array, + out_dx: Optional[_I16Array] = ..., + out_dy: Optional[_I16Array] = ..., +) -> tuple[_I16Array, _I16Array]: ... +def sobel_l1_magnitude_image( + image: _U8Array, + out: Optional[_I16Array] = ..., +) -> _I16Array: ... def scharr_image( image: _U8Array, dx: int, diff --git a/crates/spatialrust-py/src/lib.rs b/crates/spatialrust-py/src/lib.rs index 62eb695..f3fbcbb 100644 --- a/crates/spatialrust-py/src/lib.rs +++ b/crates/spatialrust-py/src/lib.rs @@ -92,7 +92,10 @@ use spatialrust::vision::{ pyr_down as pyr_down_op, pyr_up as pyr_up_op, remap as remap_op, resize as resize_op, resize_into as resize_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, soft_nms as soft_nms_op, solve_pnp as solve_pnp_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, + 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, threshold as threshold_op, AbsolutePose, AdaptiveThresholdMethod, BinaryMask, BorderMode, BoundingBox2, CameraMatrix3, CannyOptions, ConfidenceMap, Connectivity, CornerSelectionOptions, @@ -2146,6 +2149,121 @@ fn sobel_image<'py>( Ok(array.into_pyarray_bound(py)) } +/// Computes exact paired 3x3 Sobel X/Y gradients from grayscale uint8 input. +#[pyfunction] +#[pyo3(signature = (image, out_dx=None, out_dy=None))] +fn spatial_gradient_image<'py>( + py: Python<'py>, + image: PyReadonlyArray2<'_, u8>, + out_dx: Option>>, + out_dy: Option>>, +) -> PyResult<(Bound<'py, PyArray2>, Bound<'py, PyArray2>)> { + let image = gray_u8_image_from_numpy(image)?; + 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)?; + let gradient_x = + Array2::from_shape_vec((image.height(), image.width()), gradient_x.into_vec()) + .map_err(to_py_err)? + .into_pyarray_bound(py); + let gradient_y = + Array2::from_shape_vec((image.height(), image.width()), gradient_y.into_vec()) + .map_err(to_py_err)? + .into_pyarray_bound(py); + Ok((gradient_x, gradient_y)) + } + (Some(out_dx), Some(out_dy)) => { + { + let mut dx_rw = out_dx.try_readwrite().map_err(|_| { + PyValueError::new_err("out_dx must not overlap the gradient input") + })?; + let mut dy_rw = out_dy.try_readwrite().map_err(|_| { + PyValueError::new_err("out_dy must not overlap out_dx or the gradient input") + })?; + let mut dx_array = dx_rw.as_array_mut(); + let mut dy_array = dy_rw.as_array_mut(); + let expected = [image.height(), image.width()]; + if dx_array.shape() != expected { + return Err(PyValueError::new_err(format!( + "out_dx shape must be ({}, {}), found {:?}", + image.height(), + image.width(), + dx_array.shape() + ))); + } + if dy_array.shape() != expected { + return Err(PyValueError::new_err(format!( + "out_dy shape must be ({}, {}), found {:?}", + image.height(), + image.width(), + dy_array.shape() + ))); + } + let Some(dx_slice) = dx_array.as_slice_mut() else { + return Err(PyValueError::new_err( + "out_dx must be a contiguous int16 array of shape (H, W)", + )); + }; + let Some(dy_slice) = dy_array.as_slice_mut() else { + return Err(PyValueError::new_err( + "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)?; + } + Ok((out_dx, out_dy)) + } + _ => Err(PyValueError::new_err("out_dx and out_dy must be provided together")), + } +} + +/// Computes exact fused 3x3 Sobel L1 magnitude from grayscale uint8 input. +#[pyfunction] +#[pyo3(signature = (image, out=None))] +fn sobel_l1_magnitude_image<'py>( + py: Python<'py>, + image: PyReadonlyArray2<'_, u8>, + out: Option>>, +) -> PyResult>> { + let image = gray_u8_image_from_numpy(image)?; + if let Some(out) = out { + { + let mut out_rw = out.try_readwrite().map_err(|_| { + PyValueError::new_err("out must not overlap the Sobel magnitude 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 int16 array of shape (H, W)", + )); + }; + sobel_l1_magnitude_u8_into_op(image.view(), 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 array = Array2::from_shape_vec((image.height(), image.width()), magnitude.into_vec()) + .map_err(to_py_err)?; + Ok(array.into_pyarray_bound(py)) +} + /// Computes a signed float32 Scharr derivative from a grayscale uint8 image. #[pyfunction] #[pyo3(signature = (image, dx, dy, scale=1.0, delta=0.0))] @@ -3700,6 +3818,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!(spatial_gradient_image, m)?)?; + m.add_function(wrap_pyfunction!(sobel_l1_magnitude_image, m)?)?; m.add_function(wrap_pyfunction!(scharr_image, m)?)?; m.add_function(wrap_pyfunction!(laplacian_image, m)?)?; m.add_function(wrap_pyfunction!(pyr_down_image, m)?)?; diff --git a/crates/spatialrust-py/tests/test_bindings.py b/crates/spatialrust-py/tests/test_bindings.py index 3d44a10..0f0ad29 100644 --- a/crates/spatialrust-py/tests/test_bindings.py +++ b/crates/spatialrust-py/tests/test_bindings.py @@ -580,6 +580,57 @@ def test_advanced_filters_and_pyramid_shapes(): assert sr.pyr_up_image(down).shape == (10, 12, 3) +def test_spatial_gradient_matches_sobel_and_reuses_outputs(): + gray = np.arange(13 * 17, dtype=np.uint8).reshape(13, 17)[:, ::-1] + expected_x = sr.sobel_image(gray, 1, 0).astype(np.int16) + expected_y = sr.sobel_image(gray, 0, 1).astype(np.int16) + gradient_x, gradient_y = sr.spatial_gradient_image(gray) + assert gradient_x.dtype == np.int16 + assert gradient_y.dtype == np.int16 + np.testing.assert_array_equal(gradient_x, expected_x) + np.testing.assert_array_equal(gradient_y, expected_y) + + out_dx = np.empty(gray.shape, dtype=np.int16) + out_dy = np.empty(gray.shape, dtype=np.int16) + returned_dx, returned_dy = sr.spatial_gradient_image(gray, out_dx, out_dy) + assert returned_dx is out_dx + assert returned_dy is out_dy + np.testing.assert_array_equal(out_dx, expected_x) + np.testing.assert_array_equal(out_dy, expected_y) + + +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) + with pytest.raises(ValueError, match="provided together"): + sr.spatial_gradient_image(gray, output) + with pytest.raises(ValueError, match="out_dx shape"): + sr.spatial_gradient_image(gray, np.empty((7, 8), dtype=np.int16), output) + with pytest.raises(ValueError, match="out_dy must be a contiguous"): + sr.spatial_gradient_image(gray, output, np.empty((7, 18), dtype=np.int16)[:, ::2]) + + +def test_fused_sobel_l1_matches_paired_gradient_and_reuses_output(): + gray = np.arange(19 * 23, dtype=np.uint8).reshape(19, 23)[:, ::-1] + gradient_x, gradient_y = sr.spatial_gradient_image(gray) + expected = np.abs(gradient_x) + np.abs(gradient_y) + magnitude = sr.sobel_l1_magnitude_image(gray) + assert magnitude.dtype == np.int16 + np.testing.assert_array_equal(magnitude, expected) + + output = np.empty(gray.shape, dtype=np.int16) + assert sr.sobel_l1_magnitude_image(gray, output) is output + np.testing.assert_array_equal(output, expected) + + +def test_fused_sobel_l1_rejects_invalid_output(): + gray = np.arange(7 * 9, dtype=np.uint8).reshape(7, 9) + with pytest.raises(ValueError, match="out shape"): + sr.sobel_l1_magnitude_image(gray, np.empty((7, 8), dtype=np.int16)) + with pytest.raises(ValueError, match="contiguous"): + sr.sobel_l1_magnitude_image(gray, np.empty((7, 18), dtype=np.int16)[:, ::2]) + + def test_morphology_operations_and_noncontiguous_input(): mask = np.zeros((9, 11), dtype=np.uint8) mask[2:7, 3:8] = 255 diff --git a/crates/spatialrust-vision/Cargo.toml b/crates/spatialrust-vision/Cargo.toml index 5d8c291..40aecde 100644 --- a/crates/spatialrust-vision/Cargo.toml +++ b/crates/spatialrust-vision/Cargo.toml @@ -13,7 +13,7 @@ default = [] resize = [] preprocess = ["resize"] warp = ["resize"] -imgproc-filter = [] +imgproc-filter = ["dep:pulp", "dep:rayon"] imgproc-morphology = ["dep:rayon"] imgproc-analysis = [] imgproc-canny = ["imgproc-filter"] @@ -38,6 +38,7 @@ spatialrust-tensor = { workspace = true, optional = true } bytemuck = { workspace = true, optional = true } thiserror.workspace = true rayon = { workspace = true, optional = true } +pulp = { workspace = true, optional = true } [dev-dependencies] criterion.workspace = true diff --git a/crates/spatialrust-vision/benches/filter.rs b/crates/spatialrust-vision/benches/filter.rs index fb1c2fa..6cd5065 100644 --- a/crates/spatialrust-vision/benches/filter.rs +++ b/crates/spatialrust-vision/benches/filter.rs @@ -1,7 +1,8 @@ use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput}; use spatialrust_image::Image; use spatialrust_vision::{ - bilateral_filter, gaussian_blur, median_blur, pyr_down, sobel, BorderMode, + bilateral_filter, gaussian_blur, median_blur, pyr_down, sobel, sobel_l1_magnitude_u8, + sobel_l1_magnitude_u8_into, spatial_gradient_u8, BorderMode, }; fn benchmark_gaussian(c: &mut Criterion) { @@ -64,5 +65,39 @@ fn benchmark_advanced_filters(c: &mut Criterion) { } } -criterion_group!(benches, benchmark_gaussian, benchmark_advanced_filters); +fn benchmark_paired_sobel(c: &mut Criterion) { + for &(name, width, height) in &[("1080p", 1920, 1080), ("4k", 3840, 2160)] { + let image = Image::::try_new( + width, + height, + (0..width * height).map(|index| ((index * 37 + 11) & 255) as u8).collect(), + ) + .unwrap(); + let mut magnitude = vec![0_i16; width * height]; + let mut group = c.benchmark_group("paired_sobel_3x3"); + group.sample_size(10); + group.throughput(Throughput::Elements((width * height) as u64)); + 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("l1_allocate", name), |b| { + b.iter(|| { + sobel_l1_magnitude_u8(black_box(image.view()), BorderMode::Reflect101).unwrap() + }) + }); + group.bench_function(BenchmarkId::new("l1_reuse", name), |b| { + b.iter(|| { + sobel_l1_magnitude_u8_into( + black_box(image.view()), + BorderMode::Reflect101, + black_box(&mut magnitude), + ) + .unwrap() + }) + }); + group.finish(); + } +} + +criterion_group!(benches, benchmark_gaussian, benchmark_advanced_filters, benchmark_paired_sobel); criterion_main!(benches); diff --git a/crates/spatialrust-vision/src/advanced_filter.rs b/crates/spatialrust-vision/src/advanced_filter.rs index 9c3e913..e8df5ab 100644 --- a/crates/spatialrust-vision/src/advanced_filter.rs +++ b/crates/spatialrust-vision/src/advanced_filter.rs @@ -1,5 +1,7 @@ //! Non-linear filters, image derivatives, and Gaussian pyramids. +use pulp::Arch; +use rayon::prelude::*; use spatialrust_image::{Image, ImageView}; use crate::border::fetch; @@ -115,6 +117,332 @@ pub fn sobel( separable_filter_f32(input, &kernel_x, &kernel_y, delta, border) } +/// Computes exact 3×3 Sobel X/Y gradients together for grayscale `u8` input. +/// +/// This matches OpenCV `spatialGradient`: outputs are signed `i16`, the two +/// first derivatives share one source traversal, and only replicated or +/// Reflect101 borders are accepted. CPU storage remains caller-owned and no +/// device transfer is performed. +pub fn spatial_gradient_u8( + input: ImageView<'_, u8, 1>, + border: BorderMode, +) -> VisionResult<(Image, Image)> { + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("spatial gradient size overflows".into()))?; + let mut gradient_x = vec![0; len]; + let mut gradient_y = vec![0; len]; + spatial_gradient_u8_into(input, border, &mut gradient_x, &mut gradient_y)?; + Ok(( + Image::try_new_with_metadata(input.width(), input.height(), gradient_x, input.metadata())?, + Image::try_new_with_metadata(input.width(), input.height(), gradient_y, input.metadata())?, + )) +} + +/// Computes exact paired 3×3 Sobel gradients into caller-owned packed output. +/// +/// Both output slices must contain exactly `input.width() * input.height()` +/// elements and must not overlap each other. Safe Rust borrowing prevents the +/// `u8` input from aliasing either `i16` output. +pub fn spatial_gradient_u8_into( + input: ImageView<'_, u8, 1>, + border: BorderMode, + gradient_x: &mut [i16], + gradient_y: &mut [i16], +) -> VisionResult<()> { + if !matches!(border, BorderMode::Replicate | BorderMode::Reflect101) { + return Err(VisionError::InvalidParameter( + "spatial gradient supports only Replicate and Reflect101 borders".into(), + )); + } + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("spatial gradient size overflows".into()))?; + validate_gradient_output(gradient_x, len, "gradient_x")?; + validate_gradient_output(gradient_y, len, "gradient_y")?; + if len == 0 { + return Ok(()); + } + + let width = input.width(); + let height = input.height(); + let arch = Arch::new(); + if len >= 1_000_000 && height > 1 { + let workers = rayon::current_num_threads().min(height); + let rows_per_worker = height.div_ceil(workers); + gradient_x + .par_chunks_mut(rows_per_worker * width) + .zip(gradient_y.par_chunks_mut(rows_per_worker * width)) + .enumerate() + .for_each(|(chunk, (gradient_x, gradient_y))| { + arch.dispatch(|| { + spatial_gradient_rows( + input, + border, + chunk * rows_per_worker, + gradient_x, + gradient_y, + ); + }); + }); + } else { + arch.dispatch(|| spatial_gradient_rows(input, border, 0, gradient_x, gradient_y)); + } + Ok(()) +} + +/// Computes exact 3×3 Sobel L1 magnitude (`|Gx| + |Gy|`) for grayscale `u8`. +/// +/// The non-negative signed `i16` output ranges from 0 through 2040. Fusing the +/// paired derivatives and magnitude avoids materializing two gradient images. +pub fn sobel_l1_magnitude_u8( + input: ImageView<'_, u8, 1>, + border: BorderMode, +) -> VisionResult> { + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("Sobel magnitude size overflows".into()))?; + let mut magnitude = vec![0; len]; + sobel_l1_magnitude_u8_into(input, border, &mut magnitude)?; + Ok(Image::try_new_with_metadata(input.width(), input.height(), magnitude, input.metadata())?) +} + +/// Computes exact fused 3×3 Sobel L1 magnitude into caller-owned packed output. +pub fn sobel_l1_magnitude_u8_into( + input: ImageView<'_, u8, 1>, + border: BorderMode, + magnitude: &mut [i16], +) -> VisionResult<()> { + if !matches!(border, BorderMode::Replicate | BorderMode::Reflect101) { + return Err(VisionError::InvalidParameter( + "Sobel L1 magnitude supports only Replicate and Reflect101 borders".into(), + )); + } + let len = input + .width() + .checked_mul(input.height()) + .ok_or_else(|| VisionError::InvalidDimensions("Sobel magnitude size overflows".into()))?; + validate_gradient_output(magnitude, len, "magnitude")?; + if len == 0 { + return Ok(()); + } + let width = input.width(); + let height = input.height(); + let arch = Arch::new(); + if len >= 1_000_000 && height > 1 { + let workers = rayon::current_num_threads().min(height); + let rows_per_worker = height.div_ceil(workers); + magnitude.par_chunks_mut(rows_per_worker * width).enumerate().for_each( + |(chunk, magnitude)| { + arch.dispatch(|| { + sobel_l1_rows(input, border, chunk * rows_per_worker, magnitude); + }); + }, + ); + } else { + arch.dispatch(|| sobel_l1_rows(input, border, 0, magnitude)); + } + Ok(()) +} + +fn validate_gradient_output(output: &[i16], len: usize, name: &str) -> VisionResult<()> { + if output.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "{name} needs {len} elements, found {}", + output.len() + ))); + } + Ok(()) +} + +fn spatial_gradient_rows( + input: ImageView<'_, u8, 1>, + border: BorderMode, + start_y: usize, + gradient_x: &mut [i16], + gradient_y: &mut [i16], +) { + let width = input.width(); + let height = input.height(); + for (local_y, (gradient_x, gradient_y)) in + gradient_x.chunks_mut(width).zip(gradient_y.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("gradient row in bounds"); + let middle = input.row(y).expect("gradient row in bounds"); + let bottom = input.row(bottom_y).expect("gradient row in bounds"); + if width == 1 { + write_spatial_gradient_pixel(top, middle, bottom, 0, 0, 0, gradient_x, gradient_y); + continue; + } + let (left, _) = gradient_neighbors(0, width, border); + write_spatial_gradient_pixel(top, middle, bottom, left, 0, 1, gradient_x, gradient_y); + for ((((top, middle), bottom), gradient_x), gradient_y) in top + .windows(3) + .zip(middle.windows(3)) + .zip(bottom.windows(3)) + .zip(gradient_x[1..width - 1].iter_mut()) + .zip(gradient_y[1..width - 1].iter_mut()) + { + let top_left = i16::from(top[0]); + let top_middle = i16::from(top[1]); + let top_right = i16::from(top[2]); + let middle_left = i16::from(middle[0]); + let middle_right = i16::from(middle[2]); + let bottom_left = i16::from(bottom[0]); + let bottom_middle = i16::from(bottom[1]); + let bottom_right = i16::from(bottom[2]); + *gradient_x = top_right + 2 * middle_right + bottom_right + - top_left + - 2 * middle_left + - bottom_left; + *gradient_y = bottom_left + 2 * bottom_middle + bottom_right + - top_left + - 2 * top_middle + - top_right; + } + let x = width - 1; + let (_, right) = gradient_neighbors(x, width, border); + write_spatial_gradient_pixel(top, middle, bottom, x - 1, x, right, gradient_x, gradient_y); + } +} + +fn sobel_l1_rows( + input: ImageView<'_, u8, 1>, + border: BorderMode, + start_y: usize, + magnitude: &mut [i16], +) { + let width = input.width(); + let height = input.height(); + for (local_y, magnitude) in magnitude.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 magnitude row in bounds"); + let middle = input.row(y).expect("Sobel magnitude row in bounds"); + let bottom = input.row(bottom_y).expect("Sobel magnitude row in bounds"); + if width == 1 { + magnitude[0] = sobel_l1_pixel(top, middle, bottom, 0, 0, 0); + continue; + } + let (left, _) = gradient_neighbors(0, width, border); + magnitude[0] = sobel_l1_pixel(top, middle, bottom, left, 0, 1); + for (((top, middle), bottom), magnitude) in top + .windows(3) + .zip(middle.windows(3)) + .zip(bottom.windows(3)) + .zip(magnitude[1..width - 1].iter_mut()) + { + *magnitude = sobel_l1_window(top, middle, bottom); + } + let x = width - 1; + let (_, right) = gradient_neighbors(x, width, border); + magnitude[x] = sobel_l1_pixel(top, middle, bottom, x - 1, x, right); + } +} + +#[inline(always)] +fn sobel_l1_window(top: &[u8], middle: &[u8], bottom: &[u8]) -> i16 { + sobel_l1_values(top[0], top[1], top[2], middle[0], middle[2], bottom[0], bottom[1], bottom[2]) +} + +#[inline(always)] +#[allow(clippy::too_many_arguments)] +fn sobel_l1_pixel( + top: &[u8], + middle: &[u8], + bottom: &[u8], + left: usize, + center: usize, + right: usize, +) -> i16 { + sobel_l1_values( + top[left], + top[center], + top[right], + middle[left], + middle[right], + bottom[left], + bottom[center], + bottom[right], + ) +} + +#[inline(always)] +#[allow(clippy::too_many_arguments)] +fn sobel_l1_values( + top_left: u8, + top_middle: u8, + top_right: u8, + middle_left: u8, + middle_right: u8, + bottom_left: u8, + bottom_middle: u8, + bottom_right: u8, +) -> i16 { + let gradient_x = 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); + let gradient_y = + i16::from(bottom_left) + 2 * i16::from(bottom_middle) + i16::from(bottom_right) + - i16::from(top_left) + - 2 * i16::from(top_middle) + - i16::from(top_right); + gradient_x.abs() + gradient_y.abs() +} + +#[inline(always)] +#[allow(clippy::too_many_arguments)] +fn write_spatial_gradient_pixel( + top: &[u8], + middle: &[u8], + bottom: &[u8], + left: usize, + center: usize, + right: usize, + gradient_x: &mut [i16], + gradient_y: &mut [i16], +) { + let top_left = i16::from(top[left]); + let top_middle = i16::from(top[center]); + let top_right = i16::from(top[right]); + let middle_left = i16::from(middle[left]); + let middle_right = i16::from(middle[right]); + let bottom_left = i16::from(bottom[left]); + let bottom_middle = i16::from(bottom[center]); + let bottom_right = i16::from(bottom[right]); + gradient_x[center] = + top_right + 2 * middle_right + bottom_right - top_left - 2 * middle_left - bottom_left; + gradient_y[center] = + bottom_left + 2 * bottom_middle + bottom_right - top_left - 2 * top_middle - top_right; +} + +fn gradient_neighbors(index: usize, length: usize, border: BorderMode) -> (usize, usize) { + if length <= 1 { + return (0, 0); + } + let left = if index > 0 { + index - 1 + } else if matches!(border, BorderMode::Reflect101) { + 1 + } else { + 0 + }; + let right = if index + 1 < length { + index + 1 + } else if matches!(border, BorderMode::Reflect101) { + length - 2 + } else { + length - 1 + }; + (left, right) +} + /// Computes a 3×3 Scharr first derivative as signed `f32` values. pub fn scharr( input: ImageView<'_, T, CHANNELS>, @@ -305,6 +633,8 @@ 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, }; use crate::BorderMode; use spatialrust_image::{Image, ImageRegion}; @@ -339,6 +669,92 @@ mod tests { assert_eq!(scharr_x[(2, 1)][0], 320.0); } + #[test] + fn paired_spatial_gradient_matches_independent_sobel_for_strided_input() { + let parent = Image::::try_new( + 9, + 6, + (0..54).map(|index| ((index * 37 + 11) % 256) as u8).collect(), + ) + .unwrap(); + let input = parent.view().subview(ImageRegion::new(1, 1, 7, 4)).unwrap(); + for border in [BorderMode::Replicate, BorderMode::Reflect101] { + let (gradient_x, gradient_y) = spatial_gradient_u8(input, border).unwrap(); + let expected_x = sobel(input, 1, 0, 3, 1.0, 0.0, border).unwrap(); + let expected_y = sobel(input, 0, 1, 3, 1.0, 0.0, border).unwrap(); + assert_eq!( + gradient_x.as_slice(), + expected_x.as_slice().iter().map(|value| *value as i16).collect::>() + ); + assert_eq!( + gradient_y.as_slice(), + expected_y.as_slice().iter().map(|value| *value as i16).collect::>() + ); + assert_eq!(gradient_x.metadata(), input.metadata()); + assert_eq!(gradient_y.metadata(), input.metadata()); + } + } + + #[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(); + let mut gradient_x = vec![0; 6]; + let mut gradient_y = vec![0; 6]; + spatial_gradient_u8_into( + image.view(), + BorderMode::Reflect101, + &mut gradient_x, + &mut gradient_y, + ) + .unwrap(); + assert!(spatial_gradient_u8_into( + image.view(), + BorderMode::Reflect101, + &mut gradient_x[..5], + &mut gradient_y, + ) + .is_err()); + assert!(spatial_gradient_u8_into( + image.view(), + BorderMode::Wrap, + &mut gradient_x, + &mut gradient_y, + ) + .is_err()); + } + + #[test] + fn fused_sobel_l1_matches_paired_gradients_for_strided_input() { + let parent = Image::::try_new( + 10, + 7, + (0..70).map(|index| ((index * 53 + 7) % 256) as u8).collect(), + ) + .unwrap(); + let input = parent.view().subview(ImageRegion::new(1, 1, 8, 5)).unwrap(); + for border in [BorderMode::Replicate, BorderMode::Reflect101] { + let (gradient_x, gradient_y) = spatial_gradient_u8(input, border).unwrap(); + let magnitude = sobel_l1_magnitude_u8(input, border).unwrap(); + let expected = gradient_x + .as_slice() + .iter() + .zip(gradient_y.as_slice()) + .map(|(&x, &y)| x.abs() + y.abs()) + .collect::>(); + assert_eq!(magnitude.as_slice(), expected); + assert_eq!(magnitude.metadata(), input.metadata()); + } + } + + #[test] + fn fused_sobel_l1_into_validates_output() { + let image = Image::::try_new(3, 2, vec![0, 1, 2, 3, 4, 5]).unwrap(); + let mut output = vec![0; 6]; + sobel_l1_magnitude_u8_into(image.view(), BorderMode::Replicate, &mut output).unwrap(); + assert!(sobel_l1_magnitude_u8_into(image.view(), BorderMode::Replicate, &mut output[..5]) + .is_err()); + } + #[test] fn laplacian_of_constant_is_zero() { let image = Image::::from_pixel(7, 5, [42]).unwrap(); diff --git a/crates/spatialrust-vision/src/canny.rs b/crates/spatialrust-vision/src/canny.rs index ce44b77..e63da95 100644 --- a/crates/spatialrust-vision/src/canny.rs +++ b/crates/spatialrust-vision/src/canny.rs @@ -4,7 +4,7 @@ use std::collections::VecDeque; use spatialrust_image::{Image, ImageView}; -use crate::{sobel, BorderMode, VisionError, VisionResult}; +use crate::{sobel, spatial_gradient_u8, BorderMode, VisionError, VisionResult}; /// Validated Canny thresholds and gradient settings. #[derive(Clone, Copy, Debug, PartialEq)] @@ -71,23 +71,43 @@ pub fn canny_with_intermediates( options: CannyOptions, ) -> VisionResult { let options = options.validate()?; - let scale = if options.aperture_size == 7 { 1.0 / 16.0 } else { 1.0 }; - let gradient_x = sobel(input, 1, 0, options.aperture_size, scale, 0.0, BorderMode::Replicate)?; - let gradient_y = sobel(input, 0, 1, options.aperture_size, scale, 0.0, BorderMode::Replicate)?; let width = input.width(); let height = input.height(); let len = width .checked_mul(height) .ok_or_else(|| VisionError::InvalidDimensions("Canny image dimensions overflow".into()))?; - let gx = gradient_x - .as_slice() - .iter() - // OpenCV's CV_16S Sobel path uses cvRound, whose supported CPU paths - // round half-way values to even. This matters for aperture 7's 1/16 scale. - .map(|&value| round_i16_ties_even(value)) - .collect::>(); - let gy = - gradient_y.as_slice().iter().map(|&value| round_i16_ties_even(value)).collect::>(); + let (gradient_x, gradient_y, gx, gy) = if options.aperture_size == 3 { + let (gradient_x_i16, gradient_y_i16) = spatial_gradient_u8(input, BorderMode::Replicate)?; + let gx = gradient_x_i16.as_slice().iter().map(|&value| i32::from(value)).collect(); + let gy = gradient_y_i16.as_slice().iter().map(|&value| i32::from(value)).collect(); + let gradient_x = gradient_x_i16.as_slice().iter().map(|&value| f32::from(value)).collect(); + let gradient_y = gradient_y_i16.as_slice().iter().map(|&value| f32::from(value)).collect(); + ( + Image::try_new_with_metadata(width, height, gradient_x, input.metadata())?, + Image::try_new_with_metadata(width, height, gradient_y, input.metadata())?, + gx, + gy, + ) + } else { + let scale = if options.aperture_size == 7 { 1.0 / 16.0 } else { 1.0 }; + let gradient_x = + sobel(input, 1, 0, options.aperture_size, scale, 0.0, BorderMode::Replicate)?; + let gradient_y = + sobel(input, 0, 1, options.aperture_size, scale, 0.0, BorderMode::Replicate)?; + let gx = gradient_x + .as_slice() + .iter() + // OpenCV's CV_16S Sobel path uses cvRound, whose supported CPU paths + // round half-way values to even. This matters for aperture 7's 1/16 scale. + .map(|&value| round_i16_ties_even(value)) + .collect::>(); + let gy = gradient_y + .as_slice() + .iter() + .map(|&value| round_i16_ties_even(value)) + .collect::>(); + (gradient_x, gradient_y, gx, gy) + }; let comparison_magnitude = gx .iter() .zip(&gy) diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index c0de2aa..c6da4b9 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -577,7 +577,7 @@ backend, allocation mode, and accuracy contract. | 113 | Planned | 112 | Caller-owned outputs and reusable workspaces for multi-stage CPU vision without hidden copies | | 114 | Planned | 112–113 | Safe size-aware CPU dispatch for packed fast paths, strided fallbacks, and bounded row/tile parallelism | | 115 | Planned | 113–114 | Accelerated resize and color conversion with precomputed sampling plans and fused preprocessing experiments | -| 116 | Planned | 113–115 | Accelerated separable Gaussian and Sobel engine with cached kernels and shared gradient passes | +| 116 | In progress | 113–115 | Accelerated separable Gaussian and Sobel engine with cached kernels and shared gradient passes | | 117 | Complete | 113–116 | Sliding-window morphology engine with exact OpenCV comparison and generic-mask fallback | | 118 | Planned | 113–117 | Fused Canny fast path that avoids public intermediates unless explicitly requested | | 119 | Planned | 104, 115–118 | Explicit upload-once GPU-resident vision chain with no intermediate readback | @@ -638,7 +638,7 @@ to one implicitly, and GPU receipts must retain named upload/readback stages. | --- | --- | --- | --- | | 116A | Planned | Reuse separable-filter intermediates and cache validated Gaussian kernels | allocation and cache tests | | 116B | Planned | Split border handling from contiguous interior loops | all border-mode property tests | -| 116C | Planned | Specialized 3x3/5x5/7x7 Gaussian and paired Sobel X/Y passes | OpenCV error receipt | +| 116C | In progress | Specialized 3x3/5x5/7x7 Gaussian and paired Sobel X/Y passes | exact 3×3 paired gradients and fused L1 complete; Gaussian and 5×5/7×7 remain | | 116D | Planned | Improve Gaussian by at least 10x and Sobel by at least 5x on one canonical large profile | native timing receipt | ### Epic 117 delivery slices diff --git a/docs/site/algorithms.html b/docs/site/algorithms.html index 9a4a4c2..b3f40cf 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 outputsspatialrust-vision · preprocessCPU Resize and warpNearest, bilinear, bicubic and area resize; remap; affine and perspective warpspatialrust-vision · resize, warpCPU / GPU - Filtering2D correlation/convolution, separable/box/Gaussian, median, bilateral, Sobel, Scharr, Laplacian, Gaussian pyramidsspatialrust-vision · imgproc-filterCPU / GPU + Filtering2D correlation/convolution, separable/box/Gaussian, median, bilateral, Sobel, exact paired gradients and fused L1 magnitude, Scharr, Laplacian, Gaussian pyramidsspatialrust-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 new file mode 100644 index 0000000..bf13d16 --- /dev/null +++ b/docs/site/filtering.html @@ -0,0 +1,40 @@ + + + + + + + Filtering and Sobel · SpatialRust + + + +
+
+
+
imgproc-filter · safe CPU
+

Filtering and paired Sobel gradients

+

Generic correlation and separable filters remain available for arbitrary components and channels. Packed grayscale u8 additionally gets an exact paired 3×3 Sobel traversal and a fused L1 magnitude operation.

+
Paired-gradient APIFused L1 API
+
+ +
+

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=).

+

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.

+
+ +
+

Measured boundary

+
+
1.86×

1080p allocated fused L1 lead over OpenCV 4.13 on the recorded Windows host.

+
2.19×

4K allocated fused L1 lead with exact int16 output.

+
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.

+

Reproduce with the focused harness and inspect the dated receipt.

+
+
+ + + diff --git a/docs/site/vision2.html b/docs/site/vision2.html index 5e0a0ab..ca77ab1 100644 --- a/docs/site/vision2.html +++ b/docs/site/vision2.html @@ -53,8 +53,9 @@

Latest measured outcome

Accuracy preserved

The VGA, 1080p, and 4K canonical masks retain maximum absolute error 0.0 against OpenCV's precise L2 distance transform.

NMS post-processing

SpatialRust measured 3.22×–8.95× faster for NMS, 26.38×–97.25× for batched NMS, and 3.42×–7.40× for linear/Gaussian Soft-NMS. Indices exactly match OpenCV; Soft-NMS scores stay within 1.79e-7.

Structured-mask labeling

Run-length union-find connected components measured 2.17×–3.61× faster than OpenCV SAUF at VGA, 1080p, and 4K. Labels, areas, and boxes are exact across canonical and 320 randomized cases.

+

Fused Sobel L1

Exact 3×3 abs(Gx) + abs(Gy) avoids four OpenCV materialization stages. Allocated Python calls measured 1.86× faster at 1080p, 2.19× at 4K, and 2.42× at 8K across 300 randomized parity cases.

-

These are workload- and host-specific results. The connected-components speed claim covers structured segmentation/document masks; dense random noise is not claimed. The repository receipts contain the reproducible methodology.

+

These are workload- and host-specific results. The Sobel claim covers fused L1 magnitude allocation, not standalone paired gradients; reuse ties at 1080p and favors OpenCV at 4K/8K. The connected-components claim covers structured segmentation/document masks; dense random noise is not claimed. Repository receipts contain the reproducible methodology.

diff --git a/notes/2026-07-16_paired_sobel_l1_acceleration.md b/notes/2026-07-16_paired_sobel_l1_acceleration.md new file mode 100644 index 0000000..a339128 --- /dev/null +++ b/notes/2026-07-16_paired_sobel_l1_acceleration.md @@ -0,0 +1,75 @@ +# Epic 116 paired Sobel and fused L1 receipt (2026-07-16) + +## Outcome + +SpatialRust now exposes OpenCV-compatible paired 3×3 Sobel gradients and an +exact fused L1 magnitude operation for grayscale `u8`. The paired primitive is +used by the aperture-3 Canny path. The fused operation avoids materializing X, +Y, absolute-X, and absolute-Y images when only `abs(Gx) + abs(Gy)` is needed. + +On the recorded host, allocated fused L1 Python calls are **1.86× faster at +1080p**, **2.19× faster at 4K**, and **2.42× faster at 8K** than the equivalent +OpenCV public pipeline. This is a fusion/allocation result, not a claim that +SpatialRust's standalone paired-gradient primitive beats OpenCV +`spatialGradient`. + +## Public API + +Rust: + +- `spatial_gradient_u8` / `spatial_gradient_u8_into` +- `sobel_l1_magnitude_u8` / `sobel_l1_magnitude_u8_into` + +Python: + +- `spatial_gradient_image(image, out_dx=None, out_dy=None)` +- `sobel_l1_magnitude_image(image, out=None)` + +The paired result is signed `i16`. Fused L1 is non-negative signed `i16` in +[0, 2040]. Replicate and Reflect101 borders are supported in Rust; Python uses +Reflect101. Wrong shapes, non-contiguous outputs, partial paired outputs, and +overlap are rejected. + +## OpenCV comparison + +Environment: + +- Windows 11 `10.0.26300`, Intel Family 6 Model 158, 6 cores / 12 logical CPUs +- CPython 3.12.10, OpenCV 4.13.0, OpenCL disabled, 12 OpenCV threads +- seeded packed random `uint8`, paired/interleaved calls +- 300 additional randomized cases, including non-contiguous Python inputs +- exact `int16` equality required before timing + +OpenCV stages are `spatialGradient` → `absdiff(Gx, 0)` → `absdiff(Gy, 0)` → +`add`. SpatialRust uses one fused operation. + +| Profile | Mode | OpenCV | SpatialRust | Result | +| --- | --- | ---: | ---: | ---: | +| 1080p | allocate | 5.508 ms | 2.954 ms | **SpatialRust 1.86×** | +| 1080p | reuse | 2.377 ms | 2.369 ms | effectively tied (SpatialRust 1.004×) | +| 4K | allocate | 22.502 ms | 10.262 ms | **SpatialRust 2.19×** | +| 4K | reuse | 10.839 ms | 12.064 ms | OpenCV 1.11× | +| 8K | allocate | 97.258 ms | 40.238 ms | **SpatialRust 2.42×** | +| 8K | reuse | 39.609 ms | 40.313 ms | OpenCV 1.02× | + +The authoritative JSON was generated as +`target/opencv-sobel-l1-performance.json` by +`bench/opencv_sobel_l1_comparison/performance.py`. + +## Validation + +- exact paired gradients against independent Sobel for strided Rust views +- Replicate and Reflect101 border parity, metadata preservation, width-one and + output-length handling +- fused L1 identity against paired gradients +- Python allocated/caller-owned identity and invalid-output tests +- 300 OpenCV randomized fused-output cases +- `spatialrust-vision` feature tests and Clippy with warnings denied +- release Python extension build and focused binding tests + +## Remaining Epic 116 boundary + +OpenCV remains faster for standalone paired gradients and the existing generic +Gaussian path. Cached Gaussian kernels, reusable separable intermediates, and +specialized 5×5/7×7 passes remain planned; this receipt does not close those +parts of Epic 116.