perf(vision): fuse paired Sobel gradients

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 }

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

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

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

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,