Accelerate direct and fused Sobel

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

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

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.

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()

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] = ...,