diff --git a/README.md b/README.md index a5b6c09..a9ef571 100644 --- a/README.md +++ b/README.md @@ -158,7 +158,8 @@ ratio; these are machine-specific measurements, not universal guarantees. | 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× | | Sobel X 3×3 | OpenCV 14.38× | OpenCV 20.31× | OpenCV 23.30× | -| Morphology open 5×5 | OpenCV 805.78× | OpenCV 578.52× | OpenCV 774.40× | +| Morphology open 5×5[^morphology-2026] | OpenCV 65.6× | OpenCV 13.1× | OpenCV 18.3× | +| Morphology open 511×511[^morphology-2026] | OpenCV 2.40× | **SpatialRust 2.46×** | **SpatialRust 2.06×** | | Canny | OpenCV 10.66× | OpenCV 12.54× | OpenCV 12.65× | | Exact Euclidean distance transform, allocate | OpenCV 1.99× | OpenCV 1.85× | OpenCV 1.45× | | Exact Euclidean distance transform, reuse | OpenCV 1.02× | OpenCV 1.06× | **SpatialRust 1.07×** | @@ -170,6 +171,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. +[^morphology-2026]: Rectangular morphology was remeasured separately with + OpenCV 4.13, OpenCL off, and the same allocated Python-API timing scope. + The separable sliding min/max path is bit-exact across 980 randomized + operation cases. It removes the old 500–900× gap for 5×5 kernels and wins + on the named 511×511 large-background workload, but OpenCV remains much + faster for small rectangles. See the [focused harness](bench/opencv_morphology_comparison/) + and [receipt](notes/2026-07-15_rectangular_morphology_acceleration.md). + The EDT fast path is exact on the canonical masks and reduced the native 4K allocation benchmark from 451.63 ms to about 75 ms. With caller-owned output and [`DistanceTransformWorkspace`](https://rsasaki0109.github.io/SpatialRust/spatialrust_vision/struct.DistanceTransformWorkspace.html), diff --git a/bench/opencv_morphology_comparison/README.md b/bench/opencv_morphology_comparison/README.md new file mode 100644 index 0000000..d638c80 --- /dev/null +++ b/bench/opencv_morphology_comparison/README.md @@ -0,0 +1,15 @@ +# OpenCV rectangular morphology comparison + +This focused harness compares bit-exact grayscale rectangular opening through +the public Python APIs. It covers the common 5×5 case and a 511×511 +background-estimation/document workload where the window-area-independent +SpatialRust engine can overtake OpenCV on large images. + +```powershell +python bench/opencv_morphology_comparison/performance.py ` + --output target/opencv-morphology-performance.json +``` + +OpenCL is disabled, input is seeded packed random `uint8`, calls are paired and +interleaved, and every timing is gated by exact output equality. The report is +a machine-specific receipt, not a blanket speed guarantee. diff --git a/bench/opencv_morphology_comparison/performance.py b/bench/opencv_morphology_comparison/performance.py new file mode 100644 index 0000000..e6e955a --- /dev/null +++ b/bench/opencv_morphology_comparison/performance.py @@ -0,0 +1,166 @@ +"""Reproducible rectangular morphology comparison with OpenCV.""" + +from __future__ import annotations + +import argparse +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, 30), + "1080p": (1920, 1080, 20), + "4k": (3840, 2160, 12), +} +KERNELS = (5, 511) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser() + parser.add_argument("--output", type=Path) + parser.add_argument("--profiles", default=",".join(PROFILES)) + parser.add_argument("--kernels", default=",".join(map(str, KERNELS))) + parser.add_argument("--warmup", type=int, default=6) + return parser.parse_args() + + +def validate_randomized_cases() -> int: + rng = np.random.default_rng(117) + operations = { + "erode": cv2.MORPH_ERODE, + "dilate": cv2.MORPH_DILATE, + "open": cv2.MORPH_OPEN, + "close": cv2.MORPH_CLOSE, + "gradient": cv2.MORPH_GRADIENT, + "tophat": cv2.MORPH_TOPHAT, + "blackhat": cv2.MORPH_BLACKHAT, + } + checked = 0 + for case in range(140): + height = int(rng.integers(1, 80)) + width = int(rng.integers(1, 100)) + kernel_width = int(rng.integers(1, 18)) + kernel_height = int(rng.integers(1, 18)) + iterations = int(rng.integers(0, 3)) + image = rng.integers(0, 256, (height, width), dtype=np.uint8) + kernel = np.ones((kernel_height, kernel_width), dtype=np.uint8) + for name, code in operations.items(): + expected = cv2.morphologyEx( + image, + code, + kernel, + iterations=iterations, + borderType=cv2.BORDER_REPLICATE, + ) + actual = sr.morphology_image( + image, name, kernel_width, kernel_height, "rect", iterations + ) + if not np.array_equal(actual, expected): + raise AssertionError(f"random case {case} failed for {name}") + checked += 1 + return checked + + +def main() -> None: + args = parse_args() + profiles = [value.strip() for value in args.profiles.split(",") if value.strip()] + kernels = [int(value) for value in args.kernels.split(",") if value.strip()] + unknown = sorted(set(profiles) - PROFILES.keys()) + if unknown: + raise ValueError(f"unknown profiles: {', '.join(unknown)}") + if any(value < 1 for value in kernels): + raise ValueError("kernel sizes must be positive") + if hasattr(cv2, "ocl"): + cv2.ocl.setUseOpenCL(False) + + randomized_cases = validate_randomized_cases() + rng = np.random.default_rng(20_260_715) + results: dict[str, object] = {} + for profile in profiles: + width, height, repeats = PROFILES[profile] + image = rng.integers(0, 256, (height, width), dtype=np.uint8) + profile_results: dict[str, object] = {} + for kernel_size in kernels: + kernel = np.ones((kernel_size, kernel_size), dtype=np.uint8) + + def opencv_open() -> np.ndarray: + return cv2.morphologyEx( + image, + cv2.MORPH_OPEN, + kernel, + borderType=cv2.BORDER_REPLICATE, + ) + + def spatialrust_open() -> np.ndarray: + return sr.morphology_image( + image, "open", kernel_size, kernel_size, "rect", 1 + ) + + expected = opencv_open() + actual = spatialrust_open() + if not np.array_equal(actual, expected): + raise AssertionError(f"{profile}/{kernel_size} is not bit-exact") + _, _, opencv_timing, spatialrust_timing = timed_pair( + opencv_open, + spatialrust_open, + warmup=args.warmup, + repeats=repeats, + seed=117 + kernel_size, + min_sample_time_ms=20.0, + ) + opencv_ms = float(opencv_timing["median"]) + spatialrust_ms = float(spatialrust_timing["median"]) + profile_results[str(kernel_size)] = { + "width": width, + "height": height, + "kernel_width": kernel_size, + "kernel_height": kernel_size, + "operation": "open", + "iterations": 1, + "border": "replicate", + "exact": True, + "opencv": opencv_timing, + "spatialrust": spatialrust_timing, + "spatialrust_speedup": opencv_ms / spatialrust_ms, + "faster_implementation": ( + "spatialrust" if spatialrust_ms < opencv_ms else "opencv" + ), + } + results[profile] = profile_results + + 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-rectangular-morphology-performance", + kind="performance", + status="pass", + environment_receipt=receipt, + results={ + "methodology": { + "timing_scope": "Python API call returning an allocated uint8 image", + "paired_interleaved": True, + "minimum_sample_time_ms": 20.0, + "input": "seeded packed random uint8 grayscale", + "randomized_correctness_cases": randomized_cases, + "thread_policy": "library defaults; OpenCV thread count recorded", + }, + "profiles": results, + }, + ) + emit_report(report, args.output) + + +if __name__ == "__main__": + main() diff --git a/crates/spatialrust-py/src/lib.rs b/crates/spatialrust-py/src/lib.rs index aa1a350..6879f02 100644 --- a/crates/spatialrust-py/src/lib.rs +++ b/crates/spatialrust-py/src/lib.rs @@ -75,31 +75,32 @@ use spatialrust::vision::{ clahe as clahe_op, connected_components_u8 as label_components_u8, decode_rle as decode_mask_runs, detect_and_describe_orb as detect_and_describe_orb_op, detect_fast as detect_fast_op, detect_harris as detect_harris_op, - detect_shi_tomasi as detect_shi_tomasi_op, + detect_shi_tomasi as detect_shi_tomasi_op, dilate_rect_u8 as dilate_rect_u8_op, distance_transform_edt_u8_into as distance_transform_edt_u8_into_op, distance_transform_edt_with_spacing as distance_transform_edt_op, encode_rle as encode_mask_runs, equalize_histogram as equalize_histogram_op, - estimate_homography_ransac as estimate_homography_ransac_op, + erode_rect_u8 as erode_rect_u8_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, 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, - morphology_ex as morphology_ex_op, nms as nms_op, otsu_threshold_u8 as otsu_threshold_u8_op, - pack_chw as pack_chw_op, pack_chw_into as pack_chw_into_op, - point_map_to_point_cloud as point_map_to_cloud, 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, 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, MorphologyShape, ObjectImageCorrespondence, - OrbOptions, OrbScoreType, PanoramaOptions, PerspectiveTransform, PointCorrespondence2, - PointMap, RgbdOdometryOptions, RleOrder, RobustEstimationOptions, ShiTomasiOptions, - SoftNmsMethod, StereoBmOptions, StructuringElement, ThresholdType, + morphology_ex as morphology_ex_op, morphology_rect_u8 as morphology_rect_u8_op, nms as nms_op, + otsu_threshold_u8 as otsu_threshold_u8_op, pack_chw as pack_chw_op, + pack_chw_into as pack_chw_into_op, point_map_to_point_cloud as point_map_to_cloud, + 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, + 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, + MorphologyShape, ObjectImageCorrespondence, OrbOptions, OrbScoreType, PanoramaOptions, + PerspectiveTransform, PointCorrespondence2, PointMap, RgbdOdometryOptions, RleOrder, + RobustEstimationOptions, ShiTomasiOptions, SoftNmsMethod, StereoBmOptions, StructuringElement, + ThresholdType, }; use spatialrust::vision::{dense_flow_block_match as dense_flow_native, DenseFlowOptions}; use spatialrust::voxelize::{ @@ -704,6 +705,24 @@ fn gray_u8_image_from_numpy(array: PyReadonlyArray2<'_, u8>) -> PyResult( + array: &'a PyReadonlyArray2<'py, u8>, + packed: &'a mut Vec, +) -> PyResult> { + let shape = array.shape(); + if shape.len() != 2 { + return Err(PyValueError::new_err("expected an (H, W) uint8 grayscale array")); + } + let (height, width) = (shape[0], shape[1]); + let data = if let Ok(slice) = array.as_slice() { + slice + } else { + packed.extend(array.as_array().iter().copied()); + packed.as_slice() + }; + ImageView::new(width, height, width, data).map_err(to_py_err) +} + /// Fits pinhole intrinsics from known camera-space points and image pixels. #[pyfunction] #[pyo3(signature = (camera_points, pixels, width, height, huber_delta=2.0, max_iterations=12))] @@ -2201,7 +2220,8 @@ fn morphology_image<'py>( shape: &str, iterations: usize, ) -> PyResult>> { - let image = gray_u8_image_from_numpy(image)?; + let mut packed = Vec::new(); + let image = gray_u8_image_view_from_numpy(&image, &mut packed)?; let shape = match shape.to_ascii_lowercase().as_str() { "rect" | "rectangle" => MorphologyShape::Rect, "cross" => MorphologyShape::Cross, @@ -2216,43 +2236,77 @@ fn morphology_image<'py>( let element = StructuringElement::try_new(shape, kernel_width, kernel_height).map_err(to_py_err)?; let operation = operation.to_ascii_lowercase(); + let rect = shape == MorphologyShape::Rect; let output = match operation.as_str() { - "erode" => { - spatialrust::vision::erode(image.view(), &element, iterations, BorderMode::Replicate) - } - "dilate" => { - spatialrust::vision::dilate(image.view(), &element, iterations, BorderMode::Replicate) - } + "erode" if rect => erode_rect_u8_op(image, &element, iterations, BorderMode::Replicate), + "dilate" if rect => dilate_rect_u8_op(image, &element, iterations, BorderMode::Replicate), + "erode" => spatialrust::vision::erode(image, &element, iterations, BorderMode::Replicate), + "dilate" => spatialrust::vision::dilate(image, &element, iterations, BorderMode::Replicate), + "open" if rect => morphology_rect_u8_op( + image, + MorphologyOperation::Open, + &element, + iterations, + BorderMode::Replicate, + ), + "close" if rect => morphology_rect_u8_op( + image, + MorphologyOperation::Close, + &element, + iterations, + BorderMode::Replicate, + ), + "gradient" if rect => morphology_rect_u8_op( + image, + MorphologyOperation::Gradient, + &element, + iterations, + BorderMode::Replicate, + ), + "tophat" | "top-hat" if rect => morphology_rect_u8_op( + image, + MorphologyOperation::TopHat, + &element, + iterations, + BorderMode::Replicate, + ), + "blackhat" | "black-hat" if rect => morphology_rect_u8_op( + image, + MorphologyOperation::BlackHat, + &element, + iterations, + BorderMode::Replicate, + ), "open" => morphology_ex_op( - image.view(), + image, MorphologyOperation::Open, &element, iterations, BorderMode::Replicate, ), "close" => morphology_ex_op( - image.view(), + image, MorphologyOperation::Close, &element, iterations, BorderMode::Replicate, ), "gradient" => morphology_ex_op( - image.view(), + image, MorphologyOperation::Gradient, &element, iterations, BorderMode::Replicate, ), "tophat" | "top-hat" => morphology_ex_op( - image.view(), + image, MorphologyOperation::TopHat, &element, iterations, BorderMode::Replicate, ), "blackhat" | "black-hat" => morphology_ex_op( - image.view(), + image, MorphologyOperation::BlackHat, &element, iterations, diff --git a/crates/spatialrust-py/tests/test_bindings.py b/crates/spatialrust-py/tests/test_bindings.py index 076b02c..108e568 100644 --- a/crates/spatialrust-py/tests/test_bindings.py +++ b/crates/spatialrust-py/tests/test_bindings.py @@ -589,6 +589,16 @@ def test_morphology_operations_and_noncontiguous_input(): assert output.dtype == np.uint8 +def test_rectangular_morphology_fast_path_accepts_noncontiguous_input(): + image = np.arange(13 * 17, dtype=np.uint8).reshape(13, 17) + contiguous = image[:, ::-1].copy() + noncontiguous = image[:, ::-1] + for operation in ("erode", "dilate", "open", "close", "gradient", "tophat", "blackhat"): + expected = sr.morphology_image(contiguous, operation, 4, 6, "rect", 2) + actual = sr.morphology_image(noncontiguous, operation, 4, 6, "rect", 2) + np.testing.assert_array_equal(actual, expected) + + def test_threshold_histogram_clahe_and_integral_contracts(): image = np.arange(9 * 11, dtype=np.uint8).reshape(9, 11)[:, ::-1] assert sr.threshold_image(image, 40).shape == image.shape diff --git a/crates/spatialrust-vision/Cargo.toml b/crates/spatialrust-vision/Cargo.toml index 36ab142..5d8c291 100644 --- a/crates/spatialrust-vision/Cargo.toml +++ b/crates/spatialrust-vision/Cargo.toml @@ -14,7 +14,7 @@ resize = [] preprocess = ["resize"] warp = ["resize"] imgproc-filter = [] -imgproc-morphology = [] +imgproc-morphology = ["dep:rayon"] imgproc-analysis = [] imgproc-canny = ["imgproc-filter"] feature2d = ["imgproc-filter", "resize"] diff --git a/crates/spatialrust-vision/benches/morphology.rs b/crates/spatialrust-vision/benches/morphology.rs index 51e1705..88543f2 100644 --- a/crates/spatialrust-vision/benches/morphology.rs +++ b/crates/spatialrust-vision/benches/morphology.rs @@ -1,7 +1,8 @@ use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput}; use spatialrust_image::Image; use spatialrust_vision::{ - morphology_ex, BorderMode, MorphologyOperation, MorphologyShape, StructuringElement, + morphology_ex, morphology_rect_u8, BorderMode, MorphologyOperation, MorphologyShape, + StructuringElement, }; fn benchmark_morphology(c: &mut Criterion) { @@ -25,6 +26,36 @@ fn benchmark_morphology(c: &mut Criterion) { }); } group.finish(); + + for kernel_size in [5, 511] { + let element = + StructuringElement::try_new(MorphologyShape::Rect, kernel_size, kernel_size).unwrap(); + let mut group = + c.benchmark_group(format!("morphology_open_gray8_rect_{kernel_size}x{kernel_size}")); + group.sample_size(10); + for &(name, width, height) in + &[("vga", 640, 480), ("1080p", 1920, 1080), ("4k", 3840, 2160)] + { + let pixels = (0..width * height) + .map(|index| ((index * 73 + index / width * 29) & 255) as u8) + .collect(); + let image = Image::::try_new(width, height, pixels).unwrap(); + group.throughput(Throughput::Elements((width * height) as u64)); + group.bench_with_input(BenchmarkId::from_parameter(name), &image, |b, image| { + b.iter(|| { + morphology_rect_u8( + black_box(image.view()), + MorphologyOperation::Open, + &element, + 1, + BorderMode::Replicate, + ) + .unwrap() + }); + }); + } + group.finish(); + } } criterion_group!(benches, benchmark_morphology); diff --git a/crates/spatialrust-vision/src/morphology.rs b/crates/spatialrust-vision/src/morphology.rs index c3c0bd3..d20fd96 100644 --- a/crates/spatialrust-vision/src/morphology.rs +++ b/crates/spatialrust-vision/src/morphology.rs @@ -1,8 +1,9 @@ //! CPU mathematical morphology with explicit structuring elements and borders. +use rayon::prelude::*; use spatialrust_image::{Image, ImageView}; -use crate::border::fetch; +use crate::border::{fetch, map_index}; use crate::{BorderMode, PixelComponent, VisionError, VisionResult}; /// Built-in structuring-element geometry. @@ -191,12 +192,351 @@ pub fn morphology_ex( } } +/// Erodes a grayscale `u8` image with a rectangular element in linear time. +/// +/// This specialized path is independent of the rectangle area and preserves +/// the generic implementation's anchor, border, iteration, and metadata +/// semantics. +pub fn erode_rect_u8( + input: ImageView<'_, u8, 1>, + element: &StructuringElement, + iterations: usize, + border: BorderMode, +) -> VisionResult> { + rect_sequence(input, element, border, &[(Extreme::Minimum, iterations)]) +} + +/// Dilates a grayscale `u8` image with a rectangular element in linear time. +pub fn dilate_rect_u8( + input: ImageView<'_, u8, 1>, + element: &StructuringElement, + iterations: usize, + border: BorderMode, +) -> VisionResult> { + rect_sequence(input, element, border, &[(Extreme::Maximum, iterations)]) +} + +/// Applies composite grayscale `u8` morphology using the rectangular fast path. +pub fn morphology_rect_u8( + input: ImageView<'_, u8, 1>, + operation: MorphologyOperation, + element: &StructuringElement, + iterations: usize, + border: BorderMode, +) -> VisionResult> { + validate_rect(element)?; + match operation { + MorphologyOperation::Open => rect_sequence( + input, + element, + border, + &[(Extreme::Minimum, iterations), (Extreme::Maximum, iterations)], + ), + MorphologyOperation::Close => rect_sequence( + input, + element, + border, + &[(Extreme::Maximum, iterations), (Extreme::Minimum, iterations)], + ), + MorphologyOperation::Gradient => { + let high = dilate_rect_u8(input, element, iterations, border)?; + let low = erode_rect_u8(input, element, iterations, border)?; + subtract_u8(high, low) + } + MorphologyOperation::TopHat => { + let original = pack(input)?; + let opened = + morphology_rect_u8(input, MorphologyOperation::Open, element, iterations, border)?; + subtract_u8(original, opened) + } + MorphologyOperation::BlackHat => { + let original = pack(input)?; + let closed = + morphology_rect_u8(input, MorphologyOperation::Close, element, iterations, border)?; + subtract_u8(closed, original) + } + } +} + #[derive(Clone, Copy)] enum Extreme { Minimum, Maximum, } +fn validate_rect(element: &StructuringElement) -> VisionResult<()> { + if element.mask.iter().all(|&active| active) { + Ok(()) + } else { + Err(VisionError::InvalidParameter( + "rectangular morphology fast path requires a full rectangular mask".into(), + )) + } +} + +fn rect_sequence( + input: ImageView<'_, u8, 1>, + element: &StructuringElement, + border: BorderMode, + stages: &[(Extreme, usize)], +) -> VisionResult> { + validate_rect(element)?; + let metadata = input.metadata(); + let (width, height) = (input.width(), input.height()); + let mut current = pack(input)?.into_vec(); + let mut workspace = RectWorkspace::default(); + for &(extreme, iterations) in stages { + for _ in 0..iterations { + workspace.apply( + ¤t, + width, + height, + element.width, + element.height, + element.anchor_x, + element.anchor_y, + border, + extreme, + ); + std::mem::swap(&mut current, &mut workspace.output); + } + } + Ok(Image::try_new_with_metadata(width, height, current, metadata)?) +} + +#[derive(Default)] +struct RectWorkspace { + horizontal: Vec, + transposed: Vec, + filtered: Vec, + output: Vec, + padded: Vec, + prefix: Vec, + suffix: Vec, +} + +impl RectWorkspace { + #[allow(clippy::too_many_arguments)] + fn apply( + &mut self, + input: &[u8], + width: usize, + height: usize, + kernel_width: usize, + kernel_height: usize, + anchor_x: usize, + anchor_y: usize, + border: BorderMode, + extreme: Extreme, + ) { + let len = width * height; + if len >= 1_000_000 { + self.apply_parallel( + input, + width, + height, + kernel_width, + kernel_height, + anchor_x, + anchor_y, + border, + extreme, + ); + return; + } + self.horizontal.resize(len, 0); + for y in 0..height { + let start = y * width; + filter_line( + &input[start..start + width], + &mut self.horizontal[start..start + width], + kernel_width, + anchor_x, + border, + extreme, + &mut self.padded, + &mut self.prefix, + &mut self.suffix, + ); + } + + self.transposed.resize(len, 0); + transpose_blocked(&self.horizontal, &mut self.transposed, width, height); + self.filtered.resize(len, 0); + for x in 0..width { + let start = x * height; + filter_line( + &self.transposed[start..start + height], + &mut self.filtered[start..start + height], + kernel_height, + anchor_y, + border, + extreme, + &mut self.padded, + &mut self.prefix, + &mut self.suffix, + ); + } + self.output.resize(len, 0); + transpose_blocked(&self.filtered, &mut self.output, height, width); + } + + #[allow(clippy::too_many_arguments)] + fn apply_parallel( + &mut self, + input: &[u8], + width: usize, + height: usize, + kernel_width: usize, + kernel_height: usize, + anchor_x: usize, + anchor_y: usize, + border: BorderMode, + extreme: Extreme, + ) { + let len = width * height; + self.horizontal.resize(len, 0); + self.horizontal.par_chunks_mut(width).enumerate().for_each_init( + LineBuffers::default, + |buffers, (y, output)| { + let start = y * width; + filter_line( + &input[start..start + width], + output, + kernel_width, + anchor_x, + border, + extreme, + &mut buffers.padded, + &mut buffers.prefix, + &mut buffers.suffix, + ); + }, + ); + + self.transposed.resize(len, 0); + self.transposed.par_chunks_mut(height).enumerate().for_each(|(x, output)| { + for (y, value) in output.iter_mut().enumerate() { + *value = self.horizontal[y * width + x]; + } + }); + self.filtered.resize(len, 0); + self.filtered.par_chunks_mut(height).enumerate().for_each_init( + LineBuffers::default, + |buffers, (x, output)| { + let start = x * height; + filter_line( + &self.transposed[start..start + height], + output, + kernel_height, + anchor_y, + border, + extreme, + &mut buffers.padded, + &mut buffers.prefix, + &mut buffers.suffix, + ); + }, + ); + self.output.resize(len, 0); + self.output.par_chunks_mut(width).enumerate().for_each(|(y, output)| { + for (x, value) in output.iter_mut().enumerate() { + *value = self.filtered[x * height + y]; + } + }); + } +} + +#[derive(Default)] +struct LineBuffers { + padded: Vec, + prefix: Vec, + suffix: Vec, +} + +#[allow(clippy::too_many_arguments)] +fn filter_line( + input: &[u8], + output: &mut [u8], + kernel: usize, + anchor: usize, + border: BorderMode, + extreme: Extreme, + padded: &mut Vec, + prefix: &mut Vec, + suffix: &mut Vec, +) { + if input.is_empty() { + return; + } + let padded_len = input.len() + kernel - 1; + padded.resize(padded_len, 0); + prefix.resize(padded_len, 0); + suffix.resize(padded_len, 0); + let constant = match border { + BorderMode::Constant([value]) => value, + _ => 0, + }; + for (index, value) in padded.iter_mut().enumerate() { + let source = index as isize - anchor as isize; + *value = map_index(source, input.len(), border).map_or(constant, |mapped| input[mapped]); + } + + for block_start in (0..padded_len).step_by(kernel) { + let block_end = (block_start + kernel).min(padded_len); + prefix[block_start] = padded[block_start]; + for index in block_start + 1..block_end { + prefix[index] = extreme_u8(prefix[index - 1], padded[index], extreme); + } + suffix[block_end - 1] = padded[block_end - 1]; + for index in (block_start..block_end - 1).rev() { + suffix[index] = extreme_u8(suffix[index + 1], padded[index], extreme); + } + } + for (index, value) in output.iter_mut().enumerate() { + *value = extreme_u8(suffix[index], prefix[index + kernel - 1], extreme); + } +} + +#[inline(always)] +fn extreme_u8(left: u8, right: u8, extreme: Extreme) -> u8 { + match extreme { + Extreme::Minimum => left.min(right), + Extreme::Maximum => left.max(right), + } +} + +fn transpose_blocked(input: &[u8], output: &mut [u8], width: usize, height: usize) { + const BLOCK: usize = 32; + for y0 in (0..height).step_by(BLOCK) { + for x0 in (0..width).step_by(BLOCK) { + let y1 = (y0 + BLOCK).min(height); + let x1 = (x0 + BLOCK).min(width); + for y in y0..y1 { + for x in x0..x1 { + output[x * height + y] = input[y * width + x]; + } + } + } + } +} + +fn subtract_u8(left: Image, right: Image) -> VisionResult> { + if left.width() != right.width() || left.height() != right.height() { + return Err(VisionError::ShapeMismatch( + "morphology subtraction dimensions must match".into(), + )); + } + let (width, height, metadata) = (left.width(), left.height(), left.metadata()); + let output = left + .into_vec() + .into_iter() + .zip(right.into_vec()) + .map(|(a, b)| a.saturating_sub(b)) + .collect(); + Ok(Image::try_new_with_metadata(width, height, output, metadata)?) +} + fn repeat_extreme( input: ImageView<'_, T, CHANNELS>, element: &StructuringElement, @@ -333,7 +673,8 @@ fn fill_ellipse(mask: &mut [bool], width: usize, height: usize) { #[cfg(test)] mod tests { use super::{ - dilate, erode, morphology_ex, MorphologyOperation, MorphologyShape, StructuringElement, + dilate, dilate_rect_u8, erode, erode_rect_u8, morphology_ex, morphology_rect_u8, + MorphologyOperation, MorphologyShape, StructuringElement, }; use crate::BorderMode; use spatialrust_image::{Image, ImageRegion}; @@ -403,4 +744,91 @@ mod tests { assert!(StructuringElement::try_from_mask(2, 2, 0, 0, vec![false; 4]).is_err()); assert!(StructuringElement::try_from_mask(2, 2, 0, 0, vec![true; 3]).is_err()); } + + #[test] + fn rectangular_u8_fast_path_matches_generic_for_anchors_borders_and_iterations() { + let mut state = 0x9e37_79b9_u32; + let pixels = (0..63) + .map(|_| { + state = state.wrapping_mul(1_664_525).wrapping_add(1_013_904_223); + (state >> 24) as u8 + }) + .collect(); + let image = Image::::try_new(9, 7, pixels).unwrap(); + let elements = [ + StructuringElement::try_new_with_anchor(MorphologyShape::Rect, 1, 1, 0, 0).unwrap(), + StructuringElement::try_new_with_anchor(MorphologyShape::Rect, 4, 2, 0, 1).unwrap(), + StructuringElement::try_new_with_anchor(MorphologyShape::Rect, 5, 3, 4, 0).unwrap(), + StructuringElement::try_new_with_anchor(MorphologyShape::Rect, 11, 9, 5, 4).unwrap(), + ]; + let borders = [ + BorderMode::Constant([37]), + BorderMode::Replicate, + BorderMode::Reflect, + BorderMode::Reflect101, + BorderMode::Wrap, + ]; + + for element in &elements { + for &border in &borders { + for iterations in 0..=2 { + assert_eq!( + erode_rect_u8(image.view(), element, iterations, border).unwrap(), + erode(image.view(), element, iterations, border).unwrap() + ); + assert_eq!( + dilate_rect_u8(image.view(), element, iterations, border).unwrap(), + dilate(image.view(), element, iterations, border).unwrap() + ); + } + } + } + } + + #[test] + fn rectangular_u8_composites_match_generic() { + let image = Image::::try_new( + 7, + 5, + (0..35).map(|index| ((index * 47 + index * index * 3) % 256) as u8).collect(), + ) + .unwrap(); + let element = + StructuringElement::try_new_with_anchor(MorphologyShape::Rect, 4, 3, 1, 2).unwrap(); + for operation in [ + MorphologyOperation::Open, + MorphologyOperation::Close, + MorphologyOperation::Gradient, + MorphologyOperation::TopHat, + MorphologyOperation::BlackHat, + ] { + for iterations in 0..=2 { + assert_eq!( + morphology_rect_u8( + image.view(), + operation, + &element, + iterations, + BorderMode::Replicate, + ) + .unwrap(), + morphology_ex( + image.view(), + operation, + &element, + iterations, + BorderMode::Replicate, + ) + .unwrap() + ); + } + } + } + + #[test] + fn rectangular_fast_path_rejects_sparse_masks() { + let image = Image::::from_pixel(3, 3, [1]).unwrap(); + let cross = StructuringElement::try_new(MorphologyShape::Cross, 3, 3).unwrap(); + assert!(erode_rect_u8(image.view(), &cross, 1, BorderMode::Replicate).is_err()); + } } diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 3932ea3..02d269c 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -578,7 +578,7 @@ backend, allocation mode, and accuracy contract. | 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 | -| 117 | Planned | 113–116 | Sliding-window morphology engine with exact OpenCV comparison and generic-mask fallback | +| 117 | In progress | 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 | | 120 | Planned | 112–119 | Vision 2 cross-platform correctness, speed, memory, allocation, and transfer release gate | @@ -645,10 +645,10 @@ to one implicitly, and GPU receipts must retain named upload/readback stages. | Slice | Status | Scope | Evidence | | --- | --- | --- | --- | -| 117A | Planned | Separable sliding min/max for rectangular elements | naive-reference property tests | -| 117B | Planned | Dedicated small rectangular paths and generic Cross/Ellipse/Diamond/custom-mask fallback | shape and anchor parity | +| 117A | Complete | Separable sliding min/max for rectangular elements | generic-reference property tests and OpenCV receipt | +| 117B | Complete | Packed rectangular `u8` dispatch and generic Cross/Ellipse/Diamond/custom-mask fallback | shape, border, stride, iteration, and anchor parity | | 117C | Planned | Ping-pong workspace for iterations and composite operations | allocation and alias tests | -| 117D | Planned | Improve morphology by at least 20x on one canonical large profile | bit-exact OpenCV and timing receipt | +| 117D | Complete | Improve morphology by at least 20x on one canonical large profile | 43.8× 4K 5×5 baseline improvement; bit-exact 511×511 OpenCV wins | ### Epic 118 delivery slices diff --git a/docs/site/algorithms.html b/docs/site/algorithms.html index 6950744..9a4a4c2 100644 --- a/docs/site/algorithms.html +++ b/docs/site/algorithms.html @@ -32,7 +32,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 - MorphologyErode, dilate, open, close, gradient, top-hat, black-hat with Rect/Cross/Ellipse/Diamond/custom elementsspatialrust-vision · imgproc-morphologyCPU / GPU + 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 Dense visionExact Euclidean distance transform with reusable workspace/output; row-major run-length/union-find connected components; contours, polygon approximation, mask RLE, and dense spatial maps. Structured VGA/1080p/4K masks measured 2.17×–3.61× faster than OpenCV SAUF with exact labels, areas, and boxes across canonical plus 320 randomized cases. Harness.spatialrust-vision · denseCPU diff --git a/docs/site/morphology.html b/docs/site/morphology.html new file mode 100644 index 0000000..417d495 --- /dev/null +++ b/docs/site/morphology.html @@ -0,0 +1,52 @@ + + + + + + + Morphology · SpatialRust + + + +
+
+
+
Vision algorithm guide
+

Mathematical morphology

+

Exact erosion, dilation, opening, closing, gradient, top-hat, and black-hat with explicit elements, anchors, iterations, and border behavior.

+
+ +
+

Choose the operation

+
+

Erode / dilate

Local minimum and maximum. Use them for shrinking or expanding bright foreground regions and for conservative mask cleanup.

+

Open / close

Opening removes small bright structures; closing fills small dark gaps. Iterations follow OpenCV's erode×N then dilate×N ordering.

+

Derived contrasts

Gradient extracts boundaries; top-hat and black-hat isolate bright or dark detail relative to the selected neighborhood scale.

+
+
+ +
+

CPU dispatch

+

Packed grayscale u8 rectangles use two separable one-dimensional min/max passes. Each line is split into prefix/suffix blocks, so work is linear in image size plus border padding rather than image size times rectangle area. Large images use bounded Rayon row/column stages; transpose buffers keep the column pass cache-friendly.

+

Cross, ellipse, diamond, custom masks, other component types, channels, and strided views retain the safe generic implementation. Public CPU calls never upload to a GPU. Device morphology is a separate explicit wgpu API.

+
+ +
+

Accuracy and speed boundary

+

The OpenCV 4.13 comparison uses random packed uint8, replicated borders, allocated Python outputs, OpenCL off, and paired interleaved calls. All published timings first require bit-exact output; 980 additional randomized operation cases also pass.

+
+
2.46×

SpatialRust lead for 1080p 511×511 opening on the recorded 12-thread Windows host.

+
2.06×

SpatialRust lead for 4K 511×511 opening on the same host.

+
small kernels

OpenCV still leads 5×5. The result is a large-window crossover, not a blanket morphology claim.

+
+

Reproduce with the focused harness and inspect the dated receipt.

+
+ +
+

Further reading

+

The implementation follows the van Herk/Gil–Werman family of constant-comparison sliding extrema. Semantics are checked against the OpenCV morphology API; the optimization direction is described by Kimmel and Gil and related streaming min/max work.

+
+
+ + + diff --git a/notes/2026-07-15_rectangular_morphology_acceleration.md b/notes/2026-07-15_rectangular_morphology_acceleration.md new file mode 100644 index 0000000..c0e9ee3 --- /dev/null +++ b/notes/2026-07-15_rectangular_morphology_acceleration.md @@ -0,0 +1,82 @@ +# Rectangular morphology acceleration receipt — 2026-07-15 + +## Outcome + +SpatialRust now dispatches packed grayscale `u8` rectangular morphology to a +safe separable sliding min/max engine. The result is bit-exact with OpenCV and +reduces the recorded 4K 5×5 opening from 4599.72 ms to 105.05 ms, a **43.8× +improvement** over the previous SpatialRust Python path. + +OpenCV remains decisively faster for small rectangles. For a named large-scale +background-estimation/document workload, SpatialRust crosses over: 511×511 +opening is **2.46× faster at 1080p** and **2.06× faster at 4K** on this host. + +## Algorithm + +Rectangles are separable, so a two-dimensional minimum or maximum is computed +as horizontal and vertical one-dimensional windows. Each padded line is split +into kernel-sized blocks with forward prefix and backward suffix extrema; a +window result combines one suffix and one prefix value. Work therefore does +not multiply by the rectangular kernel area. Large images run bounded parallel +row/column stages with explicit transpose buffers. No `unsafe` code or hidden +device transfer is used. + +The existing generic implementation remains the exact fallback for other +component types, channels, strides, and Cross/Ellipse/Diamond/custom masks. + +References: + +- OpenCV morphology API and iteration semantics: +- OpenCV row/column morphology implementation: +- Kimmel and Gil, *Efficient Implementation of Min/Max Filters*: +- Lemire, *Streaming Maximum-Minimum Filter Using No More than Three Comparisons per Element*: + +## Reproduction environment + +- Windows 11 `10.0.26300`, AMD64 +- Intel Family 6 Model 158, 6 cores / 12 logical CPUs +- CPython 3.12.10 +- OpenCV 4.13.0, 12 reported threads, OpenCL disabled +- SpatialRust 1.0.0 release wheel +- seeded packed random grayscale `uint8`; `BORDER_REPLICATE` +- six warmups; paired/interleaved order; calls batched to at least 20 ms +- 30 VGA, 20 1080p, and 12 4K samples per kernel + +Run: + +```powershell +python bench/opencv_morphology_comparison/performance.py ` + --output target/opencv-morphology-performance.json +``` + +## Python API medians + +| Profile | Kernel | OpenCV | SpatialRust | Result | +| --- | ---: | ---: | ---: | ---: | +| VGA | 5×5 | 0.135 ms | 8.858 ms | OpenCV 65.6× | +| VGA | 511×511 | 8.100 ms | 19.433 ms | OpenCV 2.40× | +| 1080p | 5×5 | 1.570 ms | 20.598 ms | OpenCV 13.1× | +| 1080p | 511×511 | 55.543 ms | 22.559 ms | **SpatialRust 2.46×** | +| 4K | 5×5 | 5.741 ms | 105.051 ms | OpenCV 18.3× | +| 4K | 511×511 | 207.788 ms | 100.889 ms | **SpatialRust 2.06×** | + +The timing scope includes the Python call and allocated output. OpenCV's +kernel array is prepared outside its call; SpatialRust's public binding still +constructs and validates its rectangular element inside the timed call. + +## Correctness gates + +- exact output for every timed profile and kernel; +- 980 seeded randomized comparisons covering all seven operations, dimensions + from 1 pixel upward, odd/even rectangles, zero/two iterations, and replicated + borders; +- Rust parity against the generic implementation for explicit asymmetric + anchors, all five border modes, oversized kernels, and multiple iterations; +- contiguous and non-contiguous NumPy input parity. + +## Scope boundary + +This is not a general claim that SpatialRust morphology is faster than OpenCV. +The measured win is limited to large rectangular windows on sufficiently large +images. OpenCV's tuned SIMD small-window engine remains substantially faster; +Cross/Ellipse/Diamond/custom masks use SpatialRust's generic correctness path.