From 3c2194e058cee3cf9955cedf0379a8bde23dedee Mon Sep 17 00:00:00 2001 From: Darshana Sanjeewan Adikari Date: Thu, 12 Dec 2024 15:34:43 +0100 Subject: [PATCH] Clean up Canny more! + Gaussian int. approximation --- Cargo.lock | 3 + crates/edge_detection/Cargo.toml | 1 + crates/edge_detection/benches/bench.rs | 179 ++++++++++--------------- crates/edge_detection/src/canny.rs | 60 +++------ crates/edge_detection/src/conv.rs | 14 +- crates/edge_detection/src/gaussian.rs | 46 +++++++ 6 files changed, 151 insertions(+), 152 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index f38dc0283e..fbb07498d4 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2264,6 +2264,7 @@ dependencies = [ "nalgebra 0.32.6", "num-traits", "pprof", + "simba 0.8.1", "types", ] @@ -5898,11 +5899,13 @@ name = "ransac" version = "0.1.0" dependencies = [ "approx", + "divan", "geometry", "itertools 0.10.5", "linear_algebra", "nalgebra 0.32.6", "ordered-float", + "pprof", "rand 0.8.5", "rand_chacha 0.3.1", ] diff --git a/crates/edge_detection/Cargo.toml b/crates/edge_detection/Cargo.toml index 911cc8e33d..d7561bc8ce 100644 --- a/crates/edge_detection/Cargo.toml +++ b/crates/edge_detection/Cargo.toml @@ -15,6 +15,7 @@ geometry = { workspace = true } linear_algebra = { workspace = true } types = { workspace = true } num-traits = { workspace = true } +simba = { workspace = true } [dev-dependencies] divan = { workspace = true } diff --git a/crates/edge_detection/benches/bench.rs b/crates/edge_detection/benches/bench.rs index 930c1e8ad1..3c33d82df5 100644 --- a/crates/edge_detection/benches/bench.rs +++ b/crates/edge_detection/benches/bench.rs @@ -1,4 +1,6 @@ -use divan::{black_box, AllocProfiler, Bencher}; +use std::{env, fs::File}; + +use divan::{bench, bench_group, black_box, AllocProfiler, Bencher}; use image::GrayImage; use imageproc::{edges::canny, filter::gaussian_blur_f32, gradients::sobel_gradients}; @@ -16,12 +18,28 @@ fn main() { const GAUSSIAN_SIGMA: f32 = 1.4; const EDGE_SOURCE_TYPE: EdgeSourceType = EdgeSourceType::LumaOfYCbCr; -fn get_profiler_guard() -> ProfilerGuard<'static> { - ProfilerGuardBuilder::default() - .frequency(1000) - .blocklist(&["pthread", "vdso"]) - .build() - .unwrap() +fn get_profiler_guard() -> Option> { + if env::var("ENABLE_FLAMEGRAPH").is_ok_and(|v| v == "1") { + ProfilerGuardBuilder::default() + .frequency(1000) + .blocklist(&["pthread", "vdso"]) + .build() + .ok() + } else { + None + } +} + +fn get_flamegraph(file_name: &str, guard: Option>) { + if let Some(report) = guard.map(|guard| guard.report().build().ok()).flatten() { + let file = File::create(format!( + "{}/test_data/output/{}.svg", + env!("CARGO_MANIFEST_DIR"), + file_name + )) + .unwrap(); + report.flamegraph(file).unwrap(); + }; } fn load_test_image() -> YCbCr422Image { @@ -35,7 +53,7 @@ fn get_blurred_source_image(image: &YCbCr422Image) -> GrayImage { gaussian_blur_f32(&edges_source, 3.5) } -#[divan::bench] +#[bench] fn imageproc_sobel_gradients(bencher: Bencher) { let image = load_test_image(); let blurred = get_blurred_source_image(&image); @@ -43,7 +61,7 @@ fn imageproc_sobel_gradients(bencher: Bencher) { bencher.bench_local(move || sobel_gradients(black_box(&blurred))); } -#[divan::bench] +#[bench] fn imageproc_canny(bencher: Bencher) { let image = load_test_image(); let mono = get_edge_source_image(&image, EDGE_SOURCE_TYPE); @@ -51,7 +69,7 @@ fn imageproc_canny(bencher: Bencher) { bencher.bench_local(move || canny(black_box(&mono), 20.0, 50.0)); } -#[divan::bench] +#[bench] fn edge_source_select(bencher: Bencher) { let image = load_test_image(); @@ -59,9 +77,9 @@ fn edge_source_select(bencher: Bencher) { .bench_local(move || get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE))); } -#[divan::bench_group] +#[bench_group] mod blurring { - use divan::{black_box, Bencher}; + use divan::{bench, black_box, Bencher}; use edge_detection::{ gaussian::{ gaussian_blur_box_filter, gaussian_blur_box_filter_nalgebra, @@ -73,7 +91,7 @@ mod blurring { use crate::{load_test_image, EDGE_SOURCE_TYPE, GAUSSIAN_SIGMA}; - #[divan::bench] + #[bench] fn gaussian_blur_with_box_filter(bencher: Bencher) { let image = get_edge_source_image(black_box(&load_test_image()), EDGE_SOURCE_TYPE); bencher.bench_local(move || { @@ -84,7 +102,7 @@ mod blurring { }); } - #[divan::bench] + #[bench] fn gaussian_blur_with_box_filter_nalgebra(bencher: Bencher) { let image = get_edge_source_image(black_box(&load_test_image()), EDGE_SOURCE_TYPE); let transposed_matrix_view = grayimage_to_2d_transposed_matrix_view(&image); @@ -96,7 +114,7 @@ mod blurring { }); } - #[divan::bench] + #[bench] fn gaussian_blur_with_box_filter_nalgebra_i16_input(bencher: Bencher) { let image = get_edge_source_image(black_box(&load_test_image()), EDGE_SOURCE_TYPE); let transposed_matrix_view = grayimage_to_2d_transposed_matrix_view::(&image); @@ -108,7 +126,7 @@ mod blurring { }); } - #[divan::bench] + #[bench] fn gaussian_blur_int_approximation(bencher: Bencher) { let image = get_edge_source_image(black_box(&load_test_image()), EDGE_SOURCE_TYPE); let transposed_matrix_view = grayimage_to_2d_transposed_matrix_view::(&image); @@ -120,7 +138,7 @@ mod blurring { }); } - #[divan::bench] + #[bench] fn imageproc_blurring(bencher: Bencher) { let image = load_test_image(); let edges_source = get_edge_source_image(&image, EDGE_SOURCE_TYPE); @@ -129,11 +147,10 @@ mod blurring { } } -#[divan::bench_group] +#[bench_group] mod sobel_operator { - use std::{env, fs::File}; - use divan::{black_box, Bencher}; + use divan::{bench, black_box, Bencher}; use edge_detection::{ conv::{ direct_convolution, direct_convolution_mut, imgproc_kernel_to_matrix, @@ -146,9 +163,12 @@ mod sobel_operator { use imageproc::gradients::{vertical_sobel, HORIZONTAL_SOBEL, VERTICAL_SOBEL}; use nalgebra::DMatrix; - use crate::{get_blurred_source_image, get_profiler_guard, load_test_image, EDGE_SOURCE_TYPE}; + use crate::{ + get_blurred_source_image, get_flamegraph, get_profiler_guard, load_test_image, + EDGE_SOURCE_TYPE, + }; - #[divan::bench] + #[bench] fn direct_convolution_vertical(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -163,7 +183,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench(min_time = 20)] fn direct_convolution_mut_vertical(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -180,7 +200,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench] fn piecewise_2d_mut_sobel(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -189,11 +209,7 @@ mod sobel_operator { let kernel_vertical = [1, 2, 1]; let kernel_horizontal = [-1, 0, 1]; - let guard = if env::var("ENABLE_FLAMEGRAPH").is_ok_and(|v| v == "1") { - Some(get_profiler_guard()) - } else { - None - }; + let guard = get_profiler_guard(); bencher.bench_local(move || { let mut out = vec![0i16; transposed_matrix_view.len()]; black_box(piecewise_2d_convolution_mut::<3, u8, i32, i16>( @@ -204,17 +220,10 @@ mod sobel_operator { )); }); - if let Some(report) = guard.map(|guard| guard.report().build().ok()).flatten() { - let file = File::create(format!( - "{}/test_data/output/piecewise_horiz.svg", - env!("CARGO_MANIFEST_DIR") - )) - .unwrap(); - report.flamegraph(file).unwrap(); - }; + get_flamegraph("piecewise_horiz", guard); } - #[divan::bench] + #[bench] fn piecewise_vertical_mut_sobel(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -232,7 +241,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench] fn piecewise_horizontal_mut_sobel(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -250,7 +259,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench] fn direct_convolution_horizontal(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -265,7 +274,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench] fn direct_convolution_sobel_vertical_wrapper(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -278,7 +287,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench] fn direct_convolution_sobel_vertical_wrapper_i16_input(bencher: Bencher) { let image = load_test_image(); let gray = get_edge_source_image(black_box(&image), black_box(EDGE_SOURCE_TYPE)); @@ -292,7 +301,7 @@ mod sobel_operator { }); } - #[divan::bench] + #[bench] fn imageproc_sobel_vertical(bencher: Bencher) { let image = load_test_image(); let blurred = get_blurred_source_image(&image); @@ -301,12 +310,10 @@ mod sobel_operator { } } -#[divan::bench_group] +#[bench_group] mod edge_points { - use std::{env, fs::File}; - - use divan::{black_box, Bencher}; + use divan::{bench, black_box, Bencher}; use edge_detection::{ canny::non_maximum_suppression, @@ -320,17 +327,15 @@ mod edge_points { }; use nalgebra::DMatrix; - use crate::{get_profiler_guard, load_test_image, EDGE_SOURCE_TYPE, GAUSSIAN_SIGMA}; + use crate::{ + get_flamegraph, get_profiler_guard, load_test_image, EDGE_SOURCE_TYPE, GAUSSIAN_SIGMA, + }; - #[divan::bench] + #[bench(min_time = 20)] fn our_canny(bencher: Bencher) { let image = load_test_image(); - let guard = if env::var("ENABLE_FLAMEGRAPH").is_ok_and(|v| v == "1") { - Some(get_profiler_guard()) - } else { - None - }; + let guard = get_profiler_guard(); bencher.bench_local(move || { black_box(get_edges_canny( black_box(3.5), @@ -340,17 +345,10 @@ mod edge_points { EDGE_SOURCE_TYPE, )) }); - if let Some(report) = guard.map(|guard| guard.report().build().ok()).flatten() { - let file = File::create(format!( - "{}/test_data/output/edges_our_canny.svg", - env!("CARGO_MANIFEST_DIR") - )) - .unwrap(); - report.flamegraph(file).unwrap(); - }; + get_flamegraph("edges_our_canny", guard); } - #[divan::bench] + #[bench] fn imageproc_sobel_vertical(bencher: Bencher) { let image = load_test_image(); @@ -364,15 +362,11 @@ mod edge_points { }); } - #[divan::bench] + #[bench] fn direct_convolution_sobel_both_axes(bencher: Bencher) { let image = load_test_image(); - let guard = if env::var("ENABLE_FLAMEGRAPH").is_ok_and(|v| v == "1") { - Some(get_profiler_guard()) - } else { - None - }; + let guard = get_profiler_guard(); bencher.bench_local(move || { black_box(get_edges_sobel_nalgebra( black_box(3.5), @@ -382,18 +376,11 @@ mod edge_points { EDGE_SOURCE_TYPE, )) }); - if let Some(report) = guard.map(|guard| guard.report().build().ok()).flatten() { - let file = File::create(format!( - "{}/test_data/output/edges_sobel.svg", - env!("CARGO_MANIFEST_DIR") - )) - .unwrap(); - report.flamegraph(file).unwrap(); - }; + get_flamegraph("edges_sobel", guard); } - // #[divan::bench] - #[divan::bench(min_time = 10)] + // #[bench] + #[bench(min_time = 10)] fn non_maximum_suppression_our_impl(bencher: Bencher) { let image = load_test_image(); @@ -404,11 +391,7 @@ mod edge_points { let gradients_y_transposed = sobel_operator_vertical::<3, i16>(&blurred); let gradients_x_transposed = sobel_operator_horizontal::<3, i16>(&blurred); - let guard = if env::var("ENABLE_FLAMEGRAPH").is_ok_and(|v| v == "1") { - Some(get_profiler_guard()) - } else { - None - }; + let guard = get_profiler_guard(); // let magnitudes = gradients_x_transposed.zip_map(&gradients_y_transposed, |x, y| { // (x * x) as i32 + (y * y) as i32 @@ -423,17 +406,10 @@ mod edge_points { 20, )); }); - if let Some(report) = guard.map(|guard| guard.report().build().ok()).flatten() { - let file = File::create(format!( - "{}/test_data/output/non_maximum_suppression_our_impl.svg", - env!("CARGO_MANIFEST_DIR") - )) - .unwrap(); - report.flamegraph(file).unwrap(); - }; + get_flamegraph("non_maximum_suppression_our_impl", guard); } - #[divan::bench] + #[bench] fn nms_synthetic(bencher: Bencher) { let angles = (0..360).map(|deg| (deg as f32).to_radians()); let (width, height) = (200, 100); @@ -485,15 +461,11 @@ mod edge_points { }); } - #[divan::bench] + #[bench] fn imageproc_canny(bencher: Bencher) { let image = load_test_image(); - let guard = if env::var("ENABLE_FLAMEGRAPH").is_ok_and(|v| v == "1") { - Some(get_profiler_guard()) - } else { - None - }; + let guard = get_profiler_guard(); bencher.bench_local(move || { black_box(get_edges_canny_imageproc( black_box(3.5), @@ -503,13 +475,6 @@ mod edge_points { EDGE_SOURCE_TYPE, )) }); - if let Some(report) = guard.map(|guard| guard.report().build().ok()).flatten() { - let file = File::create(format!( - "{}/test_data/output/edges_canny.svg", - env!("CARGO_MANIFEST_DIR") - )) - .unwrap(); - report.flamegraph(file).unwrap(); - }; + get_flamegraph("edges_canny", guard); } } diff --git a/crates/edge_detection/src/canny.rs b/crates/edge_detection/src/canny.rs index c0a1882158..401cf71786 100644 --- a/crates/edge_detection/src/canny.rs +++ b/crates/edge_detection/src/canny.rs @@ -17,11 +17,6 @@ pub fn canny( let sigma = gaussian_sigma.unwrap_or(1.4); const SOBEL_KERNEL_SIZE: usize = 3; - // TODO remove this - // let blurred = gaussian_blur_f32(image, sigma); - - // let converted = grayimage_to_2d_transposed_matrix_view::(&blurred); - let input = grayimage_to_2d_transposed_matrix_view::(image); let converted = gaussian_blur_box_filter_nalgebra(&input, sigma); @@ -56,26 +51,15 @@ pub(crate) fn get_gradient_magnitude( #[derive(Debug, Clone, Copy, PartialEq, Eq)] #[repr(u8)] -enum Octant { - FirstOctant0Deg, - SecondOctant45Deg, - ThirdOctant90Deg, - FourthOctant135Deg, +enum OctantWithDegName { + FirstOctant0, + SecondOctant45, + ThirdOctant90, + FourthOctant135, } #[inline(always)] -fn approximate_direction_integer_only(y: i16, x: i16) -> Octant { - // match (y, x) { - // (0, _) => return 0, - // (_, 0) => return 90, - // _ => {} - // } - // if y == 0 { - // return 0; - // } else if x == 0 { - // return 90; - // } - +fn approximate_direction_integer_only(y: i16, x: i16) -> OctantWithDegName { // This trick is taken from OpenCV's Canny implementation // The idea is to perform the tan22.5 and tan67.5 boundary calculations with integers only avoiding float calculations and atan2, etc // To grab the perpendicular two pixels to edge direction, only 4 of the 8 possible directions are needed (the opposide ones of each direction yields the same.) @@ -103,16 +87,16 @@ fn approximate_direction_integer_only(y: i16, x: i16) -> Octant { // check if the inequalities hold if y_shifted < tan_22_5_mul_x || y == 0 { - Octant::FirstOctant0Deg + OctantWithDegName::FirstOctant0 } else if y_shifted >= tan_67_5_mul_x { - Octant::ThirdOctant90Deg + OctantWithDegName::ThirdOctant90 } else { // let only_one_is_negative = y.is_positive() ^ x.is_positive(); // if only_one_is_negative { if y.signum() == x.signum() { - Octant::SecondOctant45Deg + OctantWithDegName::SecondOctant45 } else { - Octant::FourthOctant135Deg + OctantWithDegName::FourthOctant135 } } } @@ -156,18 +140,18 @@ pub fn non_maximum_suppression( gradients_y_slice[index], gradients_x_slice[index], ) { - Octant::FirstOctant0Deg => { + OctantWithDegName::FirstOctant0 => { pixel > flat_slice[index - 1] && pixel > flat_slice[index + 1] } - Octant::SecondOctant45Deg => { + OctantWithDegName::SecondOctant45 => { pixel > flat_slice[previous_column_point - 1] && pixel > flat_slice[next_column_point + 1] } - Octant::ThirdOctant90Deg => { + OctantWithDegName::ThirdOctant90 => { pixel > flat_slice[previous_column_point] && pixel > flat_slice[next_column_point] } - Octant::FourthOctant135Deg => { + OctantWithDegName::FourthOctant135 => { pixel > flat_slice[previous_column_point + 1] && pixel > flat_slice[next_column_point - 1] } @@ -242,7 +226,7 @@ mod tests { use std::i16; - use super::{approximate_direction_integer_only, Octant}; + use super::{approximate_direction_integer_only, OctantWithDegName}; #[test] fn non_maximum_suppression_direction_approximation() { @@ -282,12 +266,12 @@ mod tests { } } - fn approximate_direction(y: i16, x: i16) -> Octant { + fn approximate_direction(y: i16, x: i16) -> OctantWithDegName { if y == 0 { - return Octant::FirstOctant0Deg; + return OctantWithDegName::FirstOctant0; } if x == 0 { - return Octant::ThirdOctant90Deg; + return OctantWithDegName::ThirdOctant90; } let mut angle = (y as f32).atan2(x as f32).to_degrees(); @@ -296,13 +280,13 @@ mod tests { } // Clamp angle. if !(22.5..157.5).contains(&angle) { - Octant::FirstOctant0Deg + OctantWithDegName::FirstOctant0 } else if (22.5..67.5).contains(&angle) { - Octant::SecondOctant45Deg + OctantWithDegName::SecondOctant45 } else if (67.5..112.5).contains(&angle) { - Octant::ThirdOctant90Deg + OctantWithDegName::ThirdOctant90 } else if (112.5..157.5).contains(&angle) { - Octant::FourthOctant135Deg + OctantWithDegName::FourthOctant135 } else { unreachable!() } diff --git a/crates/edge_detection/src/conv.rs b/crates/edge_detection/src/conv.rs index 96943d18c6..442d28037e 100644 --- a/crates/edge_detection/src/conv.rs +++ b/crates/edge_detection/src/conv.rs @@ -26,7 +26,7 @@ where let mut result = DMatrix::::zeros(image_rows, image_cols); - direct_convolution_mut_try_again(image, &mut result.as_mut_slice(), kernel); + direct_convolution_mut_try_again(image, result.as_mut_slice(), kernel); result } @@ -223,7 +223,7 @@ pub fn piecewise_vertical_convolution_mut( transposed_image: &DMatrix, - mut dst: &mut [OutputType], + dst: &mut [OutputType], piecewise_kernel_horizontal: &[KType; KSIZE], piecewise_kernel_vertical: &[KType; KSIZE], ) where @@ -239,15 +239,15 @@ pub fn piecewise_2d_convolution_mut( - &transposed_image, - &mut dst, - &piecewise_kernel_horizontal, + transposed_image, + dst, + piecewise_kernel_horizontal, ); piecewise_vertical_convolution_mut::( &DMatrix::from_column_slice(transposed_image.nrows(), transposed_image.ncols(), dst), - &mut dst, - &piecewise_kernel_vertical, + dst, + piecewise_kernel_vertical, ); // let result_view = DMatrixView::from_slice(&out, image.nrows(), image.ncols()); diff --git a/crates/edge_detection/src/gaussian.rs b/crates/edge_detection/src/gaussian.rs index 6e01278bae..710500d07a 100644 --- a/crates/edge_detection/src/gaussian.rs +++ b/crates/edge_detection/src/gaussian.rs @@ -1,3 +1,4 @@ +use core::f32; use num_traits::{AsPrimitive, PrimInt}; use std::ops::{Div, Mul, MulAssign}; @@ -30,6 +31,36 @@ pub fn gaussian_blur_box_filter(image: &GrayImage, sigma: f32) -> ImageBuffer f32 { + ((2.0 * f32::consts::PI).sqrt() * r).recip() * (-x.powi(2) / (2.0 * r.powi(2))).exp() +} + +pub fn gaussian_blur_try_2_nalgebra( + image: &DMatrix, + _sigma: f32, +) -> DMatrix +where + InputType: Into + + AsPrimitive + + AsPrimitive + + PrimInt + + Scalar + + Mul + + MulAssign, + KernelType: PrimInt, + // OutputType: Into, +{ + // let kernel = imgproc_kernel_to_matrix::<3>(&GAUSSIAN_BLUR_3x3); + + let kernel = SMatrix::::from_row_slice(&[1, 2, 1, 2, 4, 2, 1, 2, 1]); + + let max = 16; + direct_convolution::<3, InputType, KernelType, OutputType>(image, &kernel) / max +} + type KernelType = i32; type OutputType = i16; /// Gaussian smoothing approximation with box filters @@ -117,6 +148,7 @@ mod tests { let converted = grayimage_to_2d_transposed_matrix_view(&luma8); let blurred = gaussian_blur_box_filter_nalgebra::(&converted, 3.5); + let blurred_int_approximation = gaussian_blur_try_2_nalgebra(&converted, 3.5); GrayImage::from_raw( image.width(), image.height(), @@ -127,5 +159,19 @@ mod tests { "{crate_dir}/test_data/output/gaussian_box_filter_nalgebra.png" )) .expect("The image saving should not fail"); + + GrayImage::from_raw( + image.width(), + image.height(), + blurred_int_approximation + .iter() + .map(|&v| v as u8) + .collect::>(), + ) + .unwrap() + .save(format!( + "{crate_dir}/test_data/output/gaussian_box_filter_nalgebra_int_approx.png" + )) + .expect("The image saving should not fail"); } }