From 4f5e462ed4615f31a9f98723ce7d360bb239e730 Mon Sep 17 00:00:00 2001 From: rsasaki0109 Date: Thu, 16 Jul 2026 06:29:51 +0900 Subject: [PATCH] perf(vision): accelerate Gaussian blur --- README.md | 12 +- bench/opencv_gaussian_comparison/README.md | 17 + .../opencv_gaussian_comparison/performance.py | 189 +++++++ crates/spatialrust-py/spatialrust.pyi | 1 + crates/spatialrust-py/src/lib.rs | 53 +- crates/spatialrust-py/tests/test_bindings.py | 18 +- crates/spatialrust-vision/benches/filter.rs | 41 +- crates/spatialrust-vision/src/filter.rs | 501 +++++++++++++++++- docs/ROADMAP.md | 8 +- docs/site/algorithms.html | 2 +- docs/site/filtering.html | 12 +- docs/site/vision2.html | 3 +- notes/2026-07-16_gaussian_acceleration.md | 62 +++ 13 files changed, 896 insertions(+), 23 deletions(-) create mode 100644 bench/opencv_gaussian_comparison/README.md create mode 100644 bench/opencv_gaussian_comparison/performance.py create mode 100644 notes/2026-07-16_gaussian_acceleration.md diff --git a/README.md b/README.md index f939a26..3676f06 100644 --- a/README.md +++ b/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 | diff --git a/bench/opencv_gaussian_comparison/README.md b/bench/opencv_gaussian_comparison/README.md new file mode 100644 index 0000000..b32a48c --- /dev/null +++ b/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. diff --git a/bench/opencv_gaussian_comparison/performance.py b/bench/opencv_gaussian_comparison/performance.py new file mode 100644 index 0000000..70a0b15 --- /dev/null +++ b/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() diff --git a/crates/spatialrust-py/spatialrust.pyi b/crates/spatialrust-py/spatialrust.pyi index 24af65b..2cbe37a 100644 --- a/crates/spatialrust-py/spatialrust.pyi +++ b/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( diff --git a/crates/spatialrust-py/src/lib.rs b/crates/spatialrust-py/src/lib.rs index f3fbcbb..a10708d 100644 --- a/crates/spatialrust-py/src/lib.rs +++ b/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,7 +2068,7 @@ 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>, @@ -2075,14 +2076,50 @@ fn gaussian_blur_image<'py>( kernel_height: usize, sigma_x: f64, sigma_y: Option, + out: Option>>, ) -> PyResult>> { - 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)?; diff --git a/crates/spatialrust-py/tests/test_bindings.py b/crates/spatialrust-py/tests/test_bindings.py index 0f0ad29..5c246b6 100644 --- a/crates/spatialrust-py/tests/test_bindings.py +++ b/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(): diff --git a/crates/spatialrust-vision/benches/filter.rs b/crates/spatialrust-vision/benches/filter.rs index 6cd5065..523ba2d 100644 --- a/crates/spatialrust-vision/benches/filter.rs +++ b/crates/spatialrust-vision/benches/filter.rs @@ -1,8 +1,9 @@ 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, sobel_l1_magnitude_u8, - sobel_l1_magnitude_u8_into, spatial_gradient_u8, BorderMode, + 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, }; fn benchmark_gaussian(c: &mut Criterion) { @@ -10,13 +11,47 @@ fn benchmark_gaussian(c: &mut Criterion) { group.sample_size(10); for &(name, width, height) in &[("640p", 640, 480), ("1080p", 1920, 1080), ("4k", 3840, 2160)] { let image = Image::::try_new(width, height, vec![127; width * height * 3]).unwrap(); + let mut output = vec![0_u8; width * height * 3]; + let mut workspace = GaussianBlurU8Workspace::new(); group.throughput(Throughput::Elements((width * height) as u64)); - group.bench_with_input(BenchmarkId::from_parameter(name), &image, |b, image| { + group.bench_with_input(BenchmarkId::new("legacy_allocate", name), &image, |b, image| { b.iter(|| { gaussian_blur(black_box(image.view()), 5, 5, 1.2, 1.2, BorderMode::Reflect101) .unwrap() }); }); + group.bench_with_input( + BenchmarkId::new("specialized_allocate", name), + &image, + |b, image| { + b.iter(|| { + gaussian_blur_u8( + black_box(image.view()), + 5, + 5, + 1.2, + 1.2, + BorderMode::Reflect101, + ) + .unwrap() + }); + }, + ); + group.bench_with_input(BenchmarkId::new("specialized_reuse", name), &image, |b, image| { + b.iter(|| { + gaussian_blur_u8_into( + black_box(image.view()), + 5, + 5, + 1.2, + 1.2, + BorderMode::Reflect101, + black_box(&mut output), + black_box(&mut workspace), + ) + .unwrap() + }); + }); } group.finish(); } diff --git a/crates/spatialrust-vision/src/filter.rs b/crates/spatialrust-vision/src/filter.rs index 57415e7..1452e6d 100644 --- a/crates/spatialrust-vision/src/filter.rs +++ b/crates/spatialrust-vision/src/filter.rs @@ -1,5 +1,7 @@ //! Linear CPU image filters with explicit border and kernel contracts. +use pulp::Arch; +use rayon::prelude::*; use spatialrust_image::{Image, ImageView}; use crate::border::{fetch, map_index}; @@ -242,6 +244,424 @@ pub fn gaussian_blur( separable_filter(input, &x, &y, 0.0, border) } +/// Reusable scratch storage and fixed-point kernel cache for [`gaussian_blur_u8_into`]. +/// +/// The workspace makes the intermediate allocation and kernel construction explicit. A +/// workspace may be reused across calls, but must not be shared by concurrent operations. +#[derive(Clone, Debug, Default)] +pub struct GaussianBlurU8Workspace { + horizontal: Vec, + kernel_x: Vec, + kernel_y: Vec, + kernel_x_key: Option<(usize, u64)>, + kernel_y_key: Option<(usize, u64)>, +} + +impl GaussianBlurU8Workspace { + /// Creates an empty workspace that grows on first use. + #[must_use] + pub const fn new() -> Self { + Self { + horizontal: Vec::new(), + kernel_x: Vec::new(), + kernel_y: Vec::new(), + kernel_x_key: None, + kernel_y_key: None, + } + } + + /// Current intermediate-buffer capacity in scalar channel elements. + #[must_use] + pub fn capacity(&self) -> usize { + self.horizontal.capacity() + } +} + +/// Applies a specialized 3×3, 5×5, or 7×7 Gaussian blur to interleaved `u8` input. +/// +/// The 3×3 and 5×5 paths use normalized Q8 kernels, a `u16` horizontal intermediate, and one +/// final rounded Q16 conversion. The 7×7 contract retains the high-precision separable path; +/// OpenCV comparison is gated to a maximum difference of two `u8` levels across broad inputs. +pub fn gaussian_blur_u8( + input: ImageView<'_, u8, CHANNELS>, + kernel_width: usize, + kernel_height: usize, + sigma_x: f64, + sigma_y: f64, + border: BorderMode, +) -> VisionResult> { + validate_specialized_gaussian_size(kernel_width)?; + validate_specialized_gaussian_size(kernel_height)?; + validate_gaussian_sigma(sigma_x)?; + validate_gaussian_sigma(sigma_y)?; + if kernel_width == 7 || kernel_height == 7 { + return gaussian_blur(input, kernel_width, kernel_height, sigma_x, sigma_y, border); + } + let len = gaussian_output_len::(input.width(), input.height())?; + let mut output = vec![0; len]; + let mut workspace = GaussianBlurU8Workspace::new(); + gaussian_blur_u8_into( + input, + kernel_width, + kernel_height, + sigma_x, + sigma_y, + border, + &mut output, + &mut workspace, + )?; + Ok(Image::try_new_with_metadata(input.width(), input.height(), output, input.metadata())?) +} + +/// Applies the specialized `u8` Gaussian blur into caller-owned packed output and workspace. +#[allow(clippy::too_many_arguments)] +pub fn gaussian_blur_u8_into( + input: ImageView<'_, u8, CHANNELS>, + kernel_width: usize, + kernel_height: usize, + sigma_x: f64, + sigma_y: f64, + border: BorderMode, + output: &mut [u8], + workspace: &mut GaussianBlurU8Workspace, +) -> VisionResult<()> { + let len = gaussian_output_len::(input.width(), input.height())?; + if output.len() != len { + return Err(VisionError::ShapeMismatch(format!( + "Gaussian output needs {len} elements, found {}", + output.len() + ))); + } + validate_specialized_gaussian_size(kernel_width)?; + validate_specialized_gaussian_size(kernel_height)?; + validate_gaussian_sigma(sigma_x)?; + validate_gaussian_sigma(sigma_y)?; + if len == 0 { + return Ok(()); + } + if kernel_width == 7 || kernel_height == 7 { + let high_precision = + gaussian_blur(input, kernel_width, kernel_height, sigma_x, sigma_y, border)?; + output.copy_from_slice(high_precision.as_slice()); + return Ok(()); + } + + prepare_fixed_gaussian_kernel( + kernel_width, + sigma_x, + &mut workspace.kernel_x_key, + &mut workspace.kernel_x, + ); + prepare_fixed_gaussian_kernel( + kernel_height, + sigma_y, + &mut workspace.kernel_y_key, + &mut workspace.kernel_y, + ); + workspace.horizontal.resize(len, 0); + let row_len = input.width() * CHANNELS; + let horizontal = &mut workspace.horizontal; + let kernel_x = workspace.kernel_x.as_slice(); + let arch = Arch::new(); + if len >= 1_000_000 { + horizontal.par_chunks_mut(row_len).enumerate().for_each(|(y, row)| { + arch.dispatch(|| gaussian_horizontal_row(input, y, kernel_x, border, row)); + }); + } else { + for (y, row) in horizontal.chunks_mut(row_len).enumerate() { + arch.dispatch(|| gaussian_horizontal_row(input, y, kernel_x, border, row)); + } + } + + let kernel_y = workspace.kernel_y.as_slice(); + if len >= 1_000_000 { + output.par_chunks_mut(row_len).enumerate().for_each(|(y, row)| { + arch.dispatch(|| { + gaussian_vertical_row( + horizontal, + input.width(), + input.height(), + y, + kernel_y, + border, + row, + ); + }); + }); + } else { + for (y, row) in output.chunks_mut(row_len).enumerate() { + arch.dispatch(|| { + gaussian_vertical_row( + horizontal, + input.width(), + input.height(), + y, + kernel_y, + border, + row, + ); + }); + } + } + Ok(()) +} + +fn gaussian_output_len(width: usize, height: usize) -> VisionResult { + width + .checked_mul(height) + .and_then(|pixels| pixels.checked_mul(CHANNELS)) + .ok_or_else(|| VisionError::InvalidDimensions("Gaussian output size overflows".into())) +} + +fn validate_specialized_gaussian_size(size: usize) -> VisionResult<()> { + if !matches!(size, 3 | 5 | 7) { + return Err(VisionError::InvalidParameter( + "specialized Gaussian kernel size must be 3, 5, or 7".into(), + )); + } + Ok(()) +} + +fn validate_gaussian_sigma(sigma: f64) -> VisionResult<()> { + if !sigma.is_finite() || sigma <= 0.0 { + return Err(VisionError::InvalidParameter( + "Gaussian sigma must be finite and positive".into(), + )); + } + Ok(()) +} + +fn prepare_fixed_gaussian_kernel( + size: usize, + sigma: f64, + cached_key: &mut Option<(usize, u64)>, + cached_kernel: &mut Vec, +) { + let key = (size, sigma.to_bits()); + if *cached_key == Some(key) { + return; + } + let center = (size / 2) as f64; + let denominator = 2.0 * sigma * sigma; + let mut weights = (0..size) + .map(|index| { + let offset = index as f64 - center; + (-(offset * offset) / denominator).exp() + }) + .collect::>(); + let sum = weights.iter().sum::(); + let mut fixed = + weights.drain(..).map(|weight| (weight * 256.0 / sum).round() as i32).collect::>(); + let adjustment = 256 - fixed.iter().sum::(); + fixed[size / 2] += adjustment; + cached_kernel.clear(); + cached_kernel.extend(fixed.into_iter().map(|weight| weight as u16)); + *cached_key = Some(key); +} + +fn gaussian_horizontal_row( + input: ImageView<'_, u8, CHANNELS>, + y: usize, + kernel: &[u16], + border: BorderMode, + output: &mut [u16], +) { + let width = input.width(); + let radius = kernel.len() / 2; + let source = input.row(y).expect("Gaussian source row in bounds"); + let constant = match border { + BorderMode::Constant(pixel) => pixel, + _ => [0; CHANNELS], + }; + if width <= radius * 2 { + for x in 0..width { + gaussian_horizontal_border_pixel(source, width, x, kernel, border, constant, output); + } + return; + } + for x in 0..radius { + gaussian_horizontal_border_pixel(source, width, x, kernel, border, constant, output); + } + for x in width - radius..width { + gaussian_horizontal_border_pixel(source, width, x, kernel, border, constant, output); + } + + let start = radius * CHANNELS; + let end = (width - radius) * CHANNELS; + match kernel { + [outer, center, _] => { + let (outer, center) = (u32::from(*outer), u32::from(*center)); + for index in start..end { + output[index] = (u32::from(source[index]) * center + + (u32::from(source[index - CHANNELS]) + u32::from(source[index + CHANNELS])) + * outer) as u16; + } + } + [outer, inner, center, _, _] => { + let (outer, inner, center) = (u32::from(*outer), u32::from(*inner), u32::from(*center)); + for index in start..end { + output[index] = (u32::from(source[index]) * center + + (u32::from(source[index - CHANNELS]) + u32::from(source[index + CHANNELS])) + * inner + + (u32::from(source[index - 2 * CHANNELS]) + + u32::from(source[index + 2 * CHANNELS])) + * outer) as u16; + } + } + [outer, middle, inner, center, _, _, _] => { + let (outer, middle, inner, center) = + (u32::from(*outer), u32::from(*middle), u32::from(*inner), u32::from(*center)); + for index in start..end { + output[index] = (u32::from(source[index]) * center + + (u32::from(source[index - CHANNELS]) + u32::from(source[index + CHANNELS])) + * inner + + (u32::from(source[index - 2 * CHANNELS]) + + u32::from(source[index + 2 * CHANNELS])) + * middle + + (u32::from(source[index - 3 * CHANNELS]) + + u32::from(source[index + 3 * CHANNELS])) + * outer) as u16; + } + } + _ => unreachable!("specialized Gaussian kernel is validated"), + } +} + +fn gaussian_horizontal_border_pixel( + source: &[u8], + width: usize, + x: usize, + kernel: &[u16], + border: BorderMode, + constant: [u8; CHANNELS], + output: &mut [u16], +) { + let radius = kernel.len() / 2; + for channel in 0..CHANNELS { + let mut sum = 0_u32; + for (tap, &weight) in kernel.iter().enumerate() { + let source_x = x as isize + tap as isize - radius as isize; + let value = map_index(source_x, width, border) + .map_or(constant[channel], |mapped| source[mapped * CHANNELS + channel]); + sum += u32::from(value) * u32::from(weight); + } + output[x * CHANNELS + channel] = sum as u16; + } +} + +#[inline(always)] +fn gaussian_round_u8(sum: u32) -> u8 { + ((sum + 32_768) >> 16).min(255) as u8 +} + +fn gaussian_vertical_border_row( + horizontal: &[u16], + width: usize, + height: usize, + y: usize, + kernel: &[u16], + border: BorderMode, + output: &mut [u8], +) { + let row_len = width * CHANNELS; + let radius = kernel.len() / 2; + let constant = match border { + BorderMode::Constant(pixel) => pixel.map(|value| u32::from(value) * 256), + _ => [0; CHANNELS], + }; + for scalar_x in 0..row_len { + let channel = scalar_x % CHANNELS; + let mut sum = 0_u32; + for (tap, &weight) in kernel.iter().enumerate() { + let source_y = y as isize + tap as isize - radius as isize; + let value = map_index(source_y, height, border).map_or(constant[channel], |mapped| { + u32::from(horizontal[mapped * row_len + scalar_x]) + }); + sum += value * u32::from(weight); + } + output[scalar_x] = gaussian_round_u8(sum); + } +} + +#[allow(clippy::too_many_arguments)] +fn gaussian_vertical_row( + horizontal: &[u16], + width: usize, + height: usize, + y: usize, + kernel: &[u16], + border: BorderMode, + output: &mut [u8], +) { + let row_len = width * CHANNELS; + let radius = kernel.len() / 2; + if y < radius || y + radius >= height { + gaussian_vertical_border_row(horizontal, width, height, y, kernel, border, output); + return; + } + match kernel { + [outer, center, _] => { + let (outer, center) = (u32::from(*outer), u32::from(*center)); + let above = (y - 1) * row_len; + let current = y * row_len; + let below = (y + 1) * row_len; + for index in 0..row_len { + output[index] = gaussian_round_u8( + u32::from(horizontal[current + index]) * center + + (u32::from(horizontal[above + index]) + + u32::from(horizontal[below + index])) + * outer, + ); + } + } + [outer, inner, center, _, _] => { + let (outer, inner, center) = (u32::from(*outer), u32::from(*inner), u32::from(*center)); + let row0 = (y - 2) * row_len; + let row1 = (y - 1) * row_len; + let row2 = y * row_len; + let row3 = (y + 1) * row_len; + let row4 = (y + 2) * row_len; + for index in 0..row_len { + output[index] = gaussian_round_u8( + u32::from(horizontal[row2 + index]) * center + + (u32::from(horizontal[row1 + index]) + + u32::from(horizontal[row3 + index])) + * inner + + (u32::from(horizontal[row0 + index]) + + u32::from(horizontal[row4 + index])) + * outer, + ); + } + } + [outer, middle, inner, center, _, _, _] => { + let (outer, middle, inner, center) = + (u32::from(*outer), u32::from(*middle), u32::from(*inner), u32::from(*center)); + let row0 = (y - 3) * row_len; + let row1 = (y - 2) * row_len; + let row2 = (y - 1) * row_len; + let row3 = y * row_len; + let row4 = (y + 1) * row_len; + let row5 = (y + 2) * row_len; + let row6 = (y + 3) * row_len; + for index in 0..row_len { + output[index] = gaussian_round_u8( + u32::from(horizontal[row3 + index]) * center + + (u32::from(horizontal[row2 + index]) + + u32::from(horizontal[row4 + index])) + * inner + + (u32::from(horizontal[row1 + index]) + + u32::from(horizontal[row5 + index])) + * middle + + (u32::from(horizontal[row0 + index]) + + u32::from(horizontal[row6 + index])) + * outer, + ); + } + } + _ => unreachable!("specialized Gaussian kernel is validated"), + } +} + fn gaussian_kernel(size: usize, sigma: f64) -> VisionResult { if size == 0 || size % 2 == 0 { return Err(VisionError::InvalidParameter( @@ -364,7 +784,8 @@ fn separable_accumulators( #[cfg(test)] mod tests { use super::{ - box_blur, convolve2d, filter2d, gaussian_blur, separable_filter, Kernel1D, Kernel2D, + box_blur, convolve2d, filter2d, gaussian_blur, gaussian_blur_u8, gaussian_blur_u8_into, + separable_filter, GaussianBlurU8Workspace, Kernel1D, Kernel2D, }; use crate::BorderMode; use spatialrust_image::{Image, ImageRegion}; @@ -419,6 +840,84 @@ mod tests { } } + #[test] + fn specialized_gaussian_matches_generic_for_strided_rgb_input() { + let width = 17; + let height = 11; + let stride = width * 3 + 7; + let mut storage = vec![199_u8; stride * height]; + for y in 0..height { + for x in 0..width { + for channel in 0..3 { + storage[y * stride + x * 3 + channel] = + ((x * 31 + y * 17 + channel * 73) & 255) as u8; + } + } + } + let view = spatialrust_image::ImageView::new(width, height, stride, &storage).unwrap(); + for &(size, sigma) in &[(3, 0.8), (5, 1.2), (7, 2.0)] { + for border in [ + BorderMode::Replicate, + BorderMode::Reflect, + BorderMode::Reflect101, + BorderMode::Wrap, + BorderMode::Constant([11, 23, 47]), + ] { + let expected = gaussian_blur(view, size, size, sigma, sigma, border).unwrap(); + let actual = gaussian_blur_u8(view, size, size, sigma, sigma, border).unwrap(); + assert!(expected + .as_slice() + .iter() + .zip(actual.as_slice()) + .all(|(&left, &right)| left.abs_diff(right) <= 1)); + } + } + } + + #[test] + fn specialized_gaussian_reuses_workspace_and_validates_output() { + let image = + Image::::try_new(8, 6, (0..144).map(|value| value as u8).collect()).unwrap(); + let mut output = vec![0; 8 * 6 * 3]; + let mut workspace = GaussianBlurU8Workspace::new(); + gaussian_blur_u8_into( + image.view(), + 5, + 5, + 1.2, + 1.2, + BorderMode::Reflect101, + &mut output, + &mut workspace, + ) + .unwrap(); + let capacity = workspace.capacity(); + gaussian_blur_u8_into( + image.view(), + 5, + 5, + 1.2, + 1.2, + BorderMode::Reflect101, + &mut output, + &mut workspace, + ) + .unwrap(); + assert_eq!(workspace.capacity(), capacity); + assert!(gaussian_blur_u8_into( + image.view(), + 5, + 5, + 1.2, + 1.2, + BorderMode::Reflect101, + &mut output[..10], + &mut workspace, + ) + .is_err()); + assert!(gaussian_blur_u8(image.view(), 9, 5, 1.2, 1.2, BorderMode::Reflect101,).is_err()); + } + #[test] fn invalid_kernels_are_rejected() { assert!(Kernel2D::try_new(0, 1, Vec::new()).is_err()); diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index c6da4b9..c22acf2 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -636,10 +636,10 @@ to one implicitly, and GPU receipts must retain named upload/readback stages. | Slice | Status | Scope | Evidence | | --- | --- | --- | --- | -| 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 | 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 | +| 116A | Complete | Reuse separable-filter intermediates and cache validated Gaussian kernels | workspace capacity and output-reuse tests | +| 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 | ### Epic 117 delivery slices diff --git a/docs/site/algorithms.html b/docs/site/algorithms.html index b3f40cf..ea50cb6 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, 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 with workspace reuse, 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 index bf13d16..3d8d135 100644 --- a/docs/site/filtering.html +++ b/docs/site/filtering.html @@ -13,14 +13,15 @@
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.

- +

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.

+

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

+

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,13 @@

Measured boundary

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.

+

Gaussian progress

+
+
20.7×

1080p allocated 5×5 improvement over the previous generic SpatialRust engine.

+
26.7×

4K allocated improvement; workspace reuse improves the old path by about 65.9×.

+
0

Canonical 5×5 max error against OpenCV; 300 mixed-size randomized cases stay within 2/255.

+
+

Standalone OpenCV Gaussian remains faster: 3.10× at 1080p, 2.93× at 4K, and 3.58× at 8K for allocated Python calls. The next optimization boundary is high-precision 7×7 and useful fused consumers, not a false standalone-win claim. Reproduce with the Gaussian harness.

diff --git a/docs/site/vision2.html b/docs/site/vision2.html index ca77ab1..6c88c34 100644 --- a/docs/site/vision2.html +++ b/docs/site/vision2.html @@ -54,8 +54,9 @@

Latest measured outcome

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.

+

Gaussian engine

Symmetric fixed-point 3×3/5×5 passes, cached kernels, unrolled interiors, and explicit workspace reuse improve the old native 5×5 path by 20.7× at 1080p and 26.7× at 4K. Canonical output is exact; standalone OpenCV remains about 2.93–3.58× faster.

-

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.

+

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. Gaussian improvement compares against SpatialRust's prior generic engine and is not an OpenCV win. 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_gaussian_acceleration.md b/notes/2026-07-16_gaussian_acceleration.md new file mode 100644 index 0000000..7413f67 --- /dev/null +++ b/notes/2026-07-16_gaussian_acceleration.md @@ -0,0 +1,62 @@ +# Specialized Gaussian acceleration + +Date: 2026-07-16 (Asia/Tokyo) + +## Scope + +This slice adds an interleaved `u8` Gaussian surface with 3×3, 5×5, and 7×7 +contracts. The accelerated 3×3/5×5 path uses symmetric normalized Q8 kernels, +a `u16` horizontal intermediate, branch-free unrolled interiors, runtime SIMD +dispatch, and bounded Rayon row stages. `GaussianBlurU8Workspace` owns scratch +and cached kernels explicitly; `gaussian_blur_u8_into` writes caller-owned +packed output. The 7×7 contract retains the high-precision generic path until +its wider accumulator specialization is complete. + +Python `gaussian_blur_image(..., out=)` now avoids copying contiguous input and +validates output shape, contiguity, and aliasing. + +## Native before/after + +Criterion on the same Windows 11, Intel 6-core/12-thread host, RGB `u8`, 5×5, +sigma 1.2, Reflect101: + +| Profile | Previous generic allocate | Specialized allocate | Improvement | Specialized workspace reuse | +| --- | ---: | ---: | ---: | ---: | +| 640×480 | 20.235 ms | 1.630 ms | 12.4× | 1.330 ms | +| 1080p | 120.460 ms | 5.812 ms | 20.7× | 2.488 ms | +| 4K | 666.160 ms | 24.946 ms | 26.7× | 10.115 ms | + +The reuse comparison versus the old allocating path is 48.4× at 1080p and +65.9× at 4K. This closes Epic 116D's 10× Gaussian improvement gate. + +## OpenCV boundary + +Focused Python timings used CPython 3.12.10, OpenCV 4.13.0, 12 threads, +OpenCL disabled, seeded interleaved samples, and minimum 20 ms batches: + +| Profile | Mode | OpenCV | SpatialRust | Outcome | +| --- | --- | ---: | ---: | --- | +| 1080p | allocate | 1.995 ms | 6.188 ms | OpenCV 3.10× | +| 4K | allocate | 7.179 ms | 21.031 ms | OpenCV 2.93× | +| 8K | allocate | 24.361 ms | 87.220 ms | OpenCV 3.58× | +| 1080p | caller output | 1.415 ms | 5.469 ms | OpenCV 3.86× | +| 4K | caller output | 5.191 ms | 20.586 ms | OpenCV 3.97× | +| 8K | caller output | 22.025 ms | 88.229 ms | OpenCV 4.01× | + +The canonical 5×5 profiles were bit-exact. Three hundred randomized 3×3, +5×5, and 7×7 RGB cases, including non-contiguous inputs, stayed within a +maximum absolute error of 2/255. The receipt deliberately does not claim a +standalone OpenCV win. It reduces the gap by over an order of magnitude and +identifies high-precision 7×7 and fused Gaussian consumers as the next work. + +Reproduce with `bench/opencv_gaussian_comparison/performance.py`; the local +JSON receipt is `target/opencv-gaussian-performance-final.json`. + +## Validation + +- five border modes and strided Rust views against the generic implementation +- 3×3, 5×5, and 7×7 sizes with metadata preservation +- workspace capacity reuse and invalid output rejection +- Python allocated/output identity, shape, contiguity, and overlap checks +- 300 randomized OpenCV cases and 1080p/4K/8K paired timings +- focused Rust tests, warnings-denied Clippy, and Criterion compilation