perf(vision): accelerate Gaussian blur

README.md

@@ -156,7 +156,7 @@ ratio; these are machine-specific measurements, not universal guarantees.
 | Bilinear resize, reuse | OpenCV 27.36× | OpenCV 145.81× | OpenCV 113.87× |
 | RGB to gray, allocate | OpenCV 11.97× | OpenCV 5.98× | OpenCV 12.98× |
 | RGB to gray, reuse | OpenCV 6.01× | OpenCV 13.81× | OpenCV 2.09× |
-| Gaussian blur 5×5 | OpenCV 139.02× | OpenCV 107.91× | OpenCV 167.28× |
+| 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× |
 | 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× |
@@ -173,6 +173,14 @@ 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.
 
+[^gaussian-2026]: The VGA cell retains the Epic 111 historical baseline. The
+  specialized 3×3/5×5 `u8` engine supersedes the 1080p/4K cells on the same
+  Windows host with OpenCV 4.13: 6.188 ms vs 1.995 ms at 1080p and 21.031 ms
+  vs 7.179 ms at 4K. At 8K it measured 87.220 ms vs 24.361 ms (OpenCV 3.58×).
+  Against SpatialRust's generic engine, native Criterion improved allocated
+  5×5 latency by 20.7× at 1080p and 26.7× at 4K; OpenCV still leads the
+  standalone operation.
+
 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
@@ -274,7 +282,7 @@ The same deterministic RGB inputs passed all VGA, 1080p, and 4K gates:
 | Bilinear resize | Exact pixels (max error 0) |
 | RGB to gray | Max error 1/255; MAE 0.1333 / 0.1333 / 0.1329 |
 | AI CHW preprocess | Max float error `5.96e-8` |
-| Gaussian blur 5×5 | Max error 1/255; 96.28% exact pixels |
+| Gaussian blur | Canonical 5×5 profiles exact; 300 randomized 3×3/5×5/7×7 cases max error 2/255 |
 | Sobel X 3×3 | Exact values (max error 0) |
 | Morphology open 5×5 | Exact pixels (max error 0) |
 | Canny | Precision, recall, F1, and IoU all 1.0 |

bench/opencv_gaussian_comparison/README.md

@@ -0,0 +1,17 @@
+# OpenCV specialized Gaussian comparison
+
+This harness compares 5×5 RGB `uint8` Gaussian blur with sigma 1.2 and
+Reflect101 borders through public Python APIs. It measures both allocated and
+caller-owned output calls at 1080p, 4K, and 8K.
+
+```powershell
+python bench/opencv_gaussian_comparison/performance.py `
+  --output target/opencv-gaussian-performance.json
+```
+
+OpenCL is disabled, OpenCV receives the logical CPU count, and paired timings
+use seeded packed input. The timing gate is backed by 300 randomized 3×3,
+5×5, and 7×7 cases, including non-contiguous inputs, with maximum absolute
+error limited to two `uint8` levels. This receipt reports standalone Gaussian
+honestly; it does not imply that SpatialRust already leads OpenCV on this
+kernel.

bench/opencv_gaussian_comparison/performance.py

@@ -0,0 +1,189 @@
+"""Reproducible specialized Gaussian blur 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, 24),
+    "4k": (3840, 2160, 16),
+    "8k": (7680, 4320, 8),
+}
+KERNELS = ((3, 0.8), (5, 1.2), (7, 2.0))
+
+
+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_blur(
+    image: np.ndarray, size: int, sigma: float, out: np.ndarray | None = None
+) -> np.ndarray:
+    return cv2.GaussianBlur(
+        image,
+        (size, size),
+        sigma,
+        dst=out,
+        sigmaY=sigma,
+        borderType=cv2.BORDER_REFLECT_101,
+    )
+
+
+def validate_randomized_cases() -> tuple[int, int]:
+    rng = np.random.default_rng(116)
+    checked = 0
+    max_error = 0
+    for case in range(300):
+        height = int(rng.integers(1, 96))
+        width = int(rng.integers(1, 128))
+        size, sigma = KERNELS[case % len(KERNELS)]
+        image = rng.integers(0, 256, (height, width, 3), dtype=np.uint8)
+        if case % 3 == 0:
+            image = image[:, ::-1]
+        expected = opencv_blur(np.ascontiguousarray(image), size, sigma)
+        actual = sr.gaussian_blur_image(image, size, size, sigma, sigma)
+        error = int(np.abs(expected.astype(np.int16) - actual.astype(np.int16)).max())
+        if error > 2:
+            raise AssertionError(f"random case {case} max error {error} exceeds 2")
+        max_error = max(max_error, error)
+        checked += 1
+    return checked, max_error
+
+
+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, randomized_max_error = 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, 3), dtype=np.uint8)
+        opencv_out = np.empty_like(image)
+        spatialrust_out = np.empty_like(image)
+
+        def opencv_allocate() -> np.ndarray:
+            return opencv_blur(image, 5, 1.2)
+
+        def spatialrust_allocate() -> np.ndarray:
+            return sr.gaussian_blur_image(image, 5, 5, 1.2, 1.2)
+
+        def opencv_reuse() -> np.ndarray:
+            return opencv_blur(image, 5, 1.2, opencv_out)
+
+        def spatialrust_reuse() -> np.ndarray:
+            return sr.gaussian_blur_image(image, 5, 5, 1.2, 1.2, out=spatialrust_out)
+
+        expected = opencv_allocate()
+        actual = spatialrust_allocate()
+        error = np.abs(expected.astype(np.int16) - actual.astype(np.int16))
+        max_error = int(error.max())
+        if max_error > 2:
+            raise AssertionError(f"{profile} max error {max_error} exceeds 2")
+        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, actual
+        ):
+            raise AssertionError(f"{profile} reuse output differs from allocated output")
+
+        _, _, 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": "RGB uint8 Gaussian blur",
+            "kernel_size": 5,
+            "sigma": 1.2,
+            "border": "reflect101",
+            "max_absolute_error": max_error,
+            "exact_fraction": float((error == 0).mean()),
+            "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-specialized-gaussian-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 RGB",
+                "randomized_correctness_cases": randomized_cases,
+                "randomized_max_absolute_error": randomized_max_error,
+                "thread_policy": "logical CPU count for OpenCV; Rayon default for SpatialRust",
+                "accuracy": "maximum absolute uint8 error <= 2",
+            },
+            "profiles": results,
+        },
+    )
+    emit_report(report, args.output)
+
+
+if __name__ == "__main__":
+    main()

crates/spatialrust-py/spatialrust.pyi

@@ -167,6 +167,7 @@ def gaussian_blur_image(
     kernel_height: int,
     sigma_x: float,
     sigma_y: Optional[float] = ...,
+    out: Optional[_U8Array] = ...,
 ) -> _U8Array: ...
 def median_blur_image(image: _U8Array, kernel_size: int) -> _U8Array: ...
 def bilateral_filter_image(

crates/spatialrust-py/src/lib.rs

@@ -82,7 +82,8 @@ use spatialrust::vision::{
     erode_rect_u8_into as erode_rect_u8_into_op,
     estimate_homography_ransac as estimate_homography_ransac_op,
     estimate_rgbd_odometry as estimate_rgbd_odometry_op, filter2d as filter2d_op,
-    find_contours as trace_contours, gaussian_blur as gaussian_blur_op,
+    find_contours as trace_contours, gaussian_blur_u8 as gaussian_blur_u8_op,
+    gaussian_blur_u8_into as gaussian_blur_u8_into_op,
     gray_world_white_balance as gray_world_white_balance_op, histogram_u8 as histogram_u8_op,
     integral_image as integral_image_op, laplacian as laplacian_op, letterbox as letterbox_op,
     match_descriptors as match_descriptors_op, median_blur as median_blur_op,
@@ -99,8 +100,8 @@ use spatialrust::vision::{
     stereo_block_match as stereo_block_match_op, stitch_panorama_pair as stitch_panorama_pair_op,
     threshold as threshold_op, AbsolutePose, AdaptiveThresholdMethod, BinaryMask, BorderMode,
     BoundingBox2, CameraMatrix3, CannyOptions, ConfidenceMap, Connectivity, CornerSelectionOptions,
-    DescriptorBuffer, Detection, DistanceTransformWorkspace, FastOptions, HarrisOptions,
-    Interpolation, Kernel2D, Keypoint2, MaskRle, MatchOptions, MorphologyOperation,
+    DescriptorBuffer, Detection, DistanceTransformWorkspace, FastOptions, GaussianBlurU8Workspace,
+    HarrisOptions, Interpolation, Kernel2D, Keypoint2, MaskRle, MatchOptions, MorphologyOperation,
     MorphologyShape, ObjectImageCorrespondence, OrbOptions, OrbScoreType, PanoramaOptions,
     PerspectiveTransform, PointCorrespondence2, PointMap, RectMorphologyWorkspace,
     RgbdOdometryOptions, RleOrder, RobustEstimationOptions, ShiTomasiOptions, SoftNmsMethod,
@@ -2067,22 +2068,58 @@ fn filter2d_image<'py>(
 
 /// Applies a normalized Gaussian blur to an RGB image using Reflect101 borders.
 #[pyfunction]
-#[pyo3(signature = (image, kernel_width, kernel_height, sigma_x, sigma_y=None))]
+#[pyo3(signature = (image, kernel_width, kernel_height, sigma_x, sigma_y=None, out=None))]
 fn gaussian_blur_image<'py>(
     py: Python<'py>,
     image: PyReadonlyArray3<'_, u8>,
     kernel_width: usize,
     kernel_height: usize,
     sigma_x: f64,
     sigma_y: Option<f64>,
+    out: Option<Bound<'py, PyArray3<u8>>>,
 ) -> PyResult<Bound<'py, PyArray3<u8>>> {
-    let image = rgb_image_from_numpy(image)?;
-    let output = gaussian_blur_op(
-        image.view(),
+    let mut packed = Vec::new();
+    let image = rgb_image_view_from_numpy(&image, &mut packed)?;
+    let sigma_y = sigma_y.unwrap_or(sigma_x);
+    if let Some(out) = out {
+        {
+            let mut out_rw = out
+                .try_readwrite()
+                .map_err(|_| PyValueError::new_err("out must not overlap the Gaussian input"))?;
+            let mut out_array = out_rw.as_array_mut();
+            if out_array.shape() != [image.height(), image.width(), 3] {
+                return Err(PyValueError::new_err(format!(
+                    "out shape must be ({}, {}, 3), 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, 3)",
+                ));
+            };
+            gaussian_blur_u8_into_op(
+                image,
+                kernel_width,
+                kernel_height,
+                sigma_x,
+                sigma_y,
+                BorderMode::Reflect101,
+                out_slice,
+                &mut GaussianBlurU8Workspace::new(),
+            )
+            .map_err(to_py_err)?;
+        }
+        return Ok(out);
+    }
+    let output = gaussian_blur_u8_op(
+        image,
         kernel_width,
         kernel_height,
         sigma_x,
-        sigma_y.unwrap_or(sigma_x),
+        sigma_y,
         BorderMode::Reflect101,
     )
     .map_err(to_py_err)?;

crates/spatialrust-py/tests/test_bindings.py

@@ -558,9 +558,25 @@ def test_filter2d_and_gaussian_preserve_rgb_shape():
     image = np.arange(7 * 9 * 3, dtype=np.uint8).reshape(7, 9, 3)
     identity = sr.filter2d_image(image[:, ::-1], np.array([[1.0]], dtype=np.float64))
     np.testing.assert_array_equal(identity, image[:, ::-1])
-    blurred = sr.gaussian_blur_image(image, 5, 3, 1.2, 0.8)
+    strided = image[:, ::-1]
+    blurred = sr.gaussian_blur_image(strided, 5, 3, 1.2, 0.8)
     assert blurred.shape == image.shape
     assert blurred.dtype == np.uint8
+    output = np.empty_like(image)
+    returned = sr.gaussian_blur_image(strided, 5, 3, 1.2, 0.8, out=output)
+    assert returned is output
+    np.testing.assert_array_equal(output, blurred)
+
+
+def test_gaussian_output_validation():
+    image = np.arange(7 * 9 * 3, dtype=np.uint8).reshape(7, 9, 3)
+    with pytest.raises(ValueError):
+        sr.gaussian_blur_image(image, 5, 5, 1.2, out=np.empty((7, 8, 3), dtype=np.uint8))
+    backing = np.empty((7, 18, 3), dtype=np.uint8)
+    with pytest.raises(ValueError):
+        sr.gaussian_blur_image(image, 5, 5, 1.2, out=backing[:, ::2])
+    with pytest.raises(ValueError):
+        sr.gaussian_blur_image(image, 5, 5, 1.2, out=image)
 
 
 def test_advanced_filters_and_pyramid_shapes():