Repository: rust-cv/cv (language: rust)

Issue:
Optimize the Akaze feature detector
While working on fixing #63, I'm going through all of Akaze's source code quite deeply, and so to not loose the understanding, I'm collecting here the low-hanging fruit optimization opportunities I'm discovering. These are:
-  [x] In `create_nonlinear_scale_space`, the `Lx` and `Ly` images could be computed in parallel, as these only access `Lsmooth` in read. This should be useful for the first evolutions when images are of a significant size.
- [x] Descriptor could be computed just after angle is calculated, in order to optimize cache locality when accessing the evolutions.
- [x] In angle and descriptor computation, window is iterator with x major, while in memory images are y major. This is bad for  cache locality.
- [x] Improved filter functions with better loop structures.
- [x] SIMD in filters.
- In `find_scale_space_extrema`, there is a O(n²) loop for testing duplicates. This could be optimized to O(n) using some form of spacing hashing.

A note here: I did not consider micro-optimizations, such as parallelizing loops within simple image computations. If so, more could be done, but the benefit is not fully clear due to the overhead cost of spawning threads. 


Discussion / hints:
@vadixidav and @xd009642 on my AMD Ryzen 9 7950X I get a factor 2x performance improvement with changes in #68, essentially by using Rayon in obvious places. I believe the code can be improved further, but serious profiling is needed to do so. Modern processors are so dependent on thermal conditions that it is hard to just compare execution timings.

For the record, I tried to replace all `azip` with their parallel versions and the performance decreased, so that is not included in #68 at the moment.

Note: the PR will be more clear once #65 is merged.
@stephanemagnenat I previously found that it was more efficient from a throughput perspective to perform akaze extraction on multiple frames in parallel rather than parallelizing the individual frames, which optimizes for latency. If you add rayon parallelization, please make rayon an optional dependency and gate out the rayon code when it isn't enabled, using ndarray or similar instead. For robotics, it may be preferable to use rayon for low latency odometry, whereas for mapping purposes it may be preferable to boost throughput.
Of course, rayon is an optional dependency. I always have first the normal code without rayon, and then the rayon-specific code (see #68). Also, my feeling at the moment is that rayon should be used smartly, for example to not spawn threads when processing tiny images.

For my application, which is related to real-time vision, latency is key. I will probably also do some more detailed profiling, because one target is WASM, where rayon is generally not available, and currently the performance of the normal code is much worst than the equivalent OpenCV CPU code.
Some profiling shows that the hottest function is `imageproc::filter::separable_filter` (~50 % of time) which is not vectorized, even in release. So now we need to find a middle ground between the incorrect previous version and the correct-but-slow version in imageproc.

Fix patch:
diff --git a/akaze/Cargo.toml b/akaze/Cargo.toml
index 42f774b6..96a27334 100644
--- a/akaze/Cargo.toml
+++ b/akaze/Cargo.toml
@@ -26,9 +26,12 @@ float-ord = { version = "0.3.2", default-features = false }
 space = "0.17.0"
 bitarray = "0.9.0"
 thiserror = { version = "1.0.40", default-features = false }
+wide = "0.7"
 serde = { version = "1.0", default-features = false, features = ["derive"], optional = true }
-imageproc = "0.23.0"
+rayon = { version = "1.7.0", optional = true }
 
+[target.'cfg(target_arch = "wasm32")'.dependencies]
+web-sys = { version = "0.3.64", features = ["Window", "Performance"], optional = true }
 
 [dev-dependencies]
 eight-point = { version = "0.8.0", path = "../eight-point" }
@@ -40,6 +43,7 @@ criterion = "0.3.5"
 pretty_env_logger = "0.4.0"
 image = "0.24"
 bitarray = { version = "0.9.0", features = ["space"] }
+imageproc = "0.23.0"
 
 [[bench]]
 name = "criterion"
diff --git a/akaze/benches/criterion.rs b/akaze/benches/criterion.rs
index 6dc62aee..f6b8337e 100644
--- a/akaze/benches/criterion.rs
+++ b/akaze/benches/criterion.rs
@@ -1,19 +1,52 @@
 use akaze::Akaze;
-use bitarray::BitArray;
 use criterion::{criterion_group, criterion_main, Criterion};
-use std::path::Path;
 
-fn image_to_kps(path: impl AsRef<Path>) -> (Vec<akaze::KeyPoint>, Vec<BitArray<64>>) {
-    Akaze::sparse().extract_path(path).unwrap()
+fn load_image() -> akaze::image::GrayFloatImage {
+    akaze::image::GrayFloatImage::from_dynamic(&image::open("../res/0000000000.png").unwrap())
 }
 
 fn extract(c: &mut Criterion) {
-    c.bench_function("extract", |b| b.iter(|| image_to_kps("res/0000000000.png")));
+    let image = load_image();
+    let akaze = Akaze::sparse();
+    c.bench_function("extract", |b| {
+        b.iter(|| akaze.extract_from_gray_float_image(&image))
+    });
 }
 
 criterion_group!(
-    name = benches;
+    name = akaze;
     config = Criterion::default().sample_size(10);
     targets = extract
 );
-criterion_main!(benches);
+
+fn bench_horizontal_filter(c: &mut Criterion) {
+    let image = load_image();
+    let small_kernel = akaze::image::gaussian_kernel(1.0, 7);
+    c.bench_function("horizontal_filter_small_kernel", |b| {
+        b.iter(|| akaze::image::horizontal_filter(&image.0, &small_kernel))
+    });
+    let large_kernel = akaze::image::gaussian_kernel(10.0, 71);
+    c.bench_function("horizontal_filter_large_kernel", |b| {
+        b.iter(|| akaze::image::horizontal_filter(&image.0, &large_kernel))
+    });
+}
+
+fn bench_vertical_filter(c: &mut Criterion) {
+    let image = load_image();
+    let small_kernel = akaze::image::gaussian_kernel(1.0, 7);
+    c.bench_function("vertical_filter_small_kernel", |b| {
+        b.iter(|| akaze::image::vertical_filter(&image.0, &small_kernel))
+    });
+    let large_kernel = akaze::image::gaussian_kernel(10.0, 71);
+    c.bench_function("vertical_filter_large_kernel", |b| {
+        b.iter(|| akaze::image::vertical_filter(&image.0, &large_kernel))
+    });
+}
+
+criterion_group!(
+    name = akaze_image;
+    config = Criterion::default().sample_size(10);
+    targets = bench_horizontal_filter, bench_vertical_filter
+);
+
+criterion_main!(akaze, akaze_image);
diff --git a/akaze/src/derivatives.rs b/akaze/src/derivatives.rs
index 337affea..16ad052c 100644
--- a/akaze/src/derivatives.rs
+++ b/akaze/src/derivatives.rs
@@ -1,21 +1,13 @@
-use crate::image::GrayFloatImage;
+use crate::image::{separable_filter, GrayFloatImage};
 
 pub fn simple_scharr_horizontal(image: &GrayFloatImage) -> GrayFloatImage {
     // similar to cv::Scharr with xorder=1, yorder=0, scale=1, delta=0
-    GrayFloatImage(imageproc::filter::separable_filter(
-        &image.0,
-        &[-1., 0., 1.],
-        &[3., 10., 3.],
-    ))
+    GrayFloatImage(separable_filter(&image.0, &[-1., 0., 1.], &[3., 10., 3.]))
 }
 
 pub fn simple_scharr_vertical(image: &GrayFloatImage) -> GrayFloatImage {
     // similar to cv::Scharr with xorder=0, yorder=1, scale=1, delta=0
-    GrayFloatImage(imageproc::filter::separable_filter(
-        &image.0,
-        &[3., 10., 3.],
-        &[-1., 0., 1.],
-    ))
+    GrayFloatImage(separable_filter(&image.0, &[3., 10., 3.], &[-1., 0., 1.]))
 }
 
 /// Compute the Scharr derivative horizontally
@@ -34,11 +26,7 @@ pub fn scharr_horizontal(image: &GrayFloatImage, sigma_size: u32) -> GrayFloatIm
     }
     let main_kernel = computer_scharr_kernel(sigma_size, FilterOrder::Main);
     let off_kernel = computer_scharr_kernel(sigma_size, FilterOrder::Off);
-    GrayFloatImage(imageproc::filter::separable_filter(
-        &image.0,
-        &main_kernel,
-        &off_kernel,
-    ))
+    GrayFloatImage(separable_filter(&image.0, &main_kernel, &off_kernel))
 }
 
 /// Compute the Scharr derivative vertically
@@ -57,11 +45,7 @@ pub fn scharr_vertical(image: &GrayFloatImage, sigma_size: u32) -> GrayFloatImag
     }
     let main_kernel = computer_scharr_kernel(sigma_size, FilterOrder::Main);
     let off_kernel = computer_scharr_kernel(sigma_size, FilterOrder::Off);
-    GrayFloatImage(imageproc::filter::separable_filter(
-        &image.0,
-        &off_kernel,
-        &main_kernel,
-    ))
+    GrayFloatImage(separable_filter(&image.0, &off_kernel, &main_kernel))
 }
 
 #[derive(Copy, Clone, Debug, PartialEq)]
diff --git a/akaze/src/descriptors.rs b/akaze/src/descriptors.rs
index 4c6971fa..7d3c30dc 100644
--- a/akaze/src/descriptors.rs
+++ b/akaze/src/descriptors.rs
@@ -1,6 +1,9 @@
 use crate::{Akaze, Error, EvolutionStep, KeyPoint};
 use bitarray::BitArray;
 
+#[cfg(feature = "rayon")]
+use rayon::prelude::*;
+
 impl Akaze {
     /// Extract descriptors from keypoints/an evolution
     ///
@@ -15,15 +18,30 @@ impl Akaze {
         evolutions: &[EvolutionStep],
         keypoints: &[KeyPoint],
     ) -> (Vec<KeyPoint>, Vec<BitArray<64>>) {
-        keypoints
-            .iter()
-            .filter_map(|&keypoint| {
-                Some((
-                    keypoint,
-                    self.get_mldb_descriptor(&keypoint, evolutions).ok()?,
-                ))
-            })
-            .unzip()
+        #[cfg(not(feature = "rayon"))]
+        {
+            keypoints
+                .iter()
+                .filter_map(|&keypoint| {
+                    Some((
+                        keypoint,
+                        self.get_mldb_descriptor(&keypoint, evolutions).ok()?,
+                    ))
+                })
+                .unzip()
+        }
+        #[cfg(feature = "rayon")]
+        {
+            keypoints
+                .par_iter()
+                .filter_map(|&keypoint| {
+                    Some((
+                        keypoint,
+                        self.get_mldb_descriptor(&keypoint, evolutions).ok()?,
+                    ))
+                })
+                .unzip()
+        }
     }
 
     /// Computes the rotation invariant M-LDB binary descriptor (maximum descriptor length)
diff --git a/akaze/src/detector_response.rs b/akaze/src/detector_response.rs
index fa34599d..696c03a1 100644
--- a/akaze/src/detector_response.rs
+++ b/akaze/src/detector_response.rs
@@ -1,15 +1,26 @@
 use crate::{derivatives, evolution::EvolutionStep, image::GrayFloatImage, Akaze};
 use ndarray::azip;
 
+#[cfg(feature = "rayon")]
+use rayon::prelude::*;
+
 impl Akaze {
     fn compute_multiscale_derivatives(&self, evolutions: &mut [EvolutionStep]) {
-        for evolution in evolutions.iter_mut() {
+        let process_evolution = |evolution: &mut EvolutionStep| {
             // The image decreases in size by a factor which is 2^octave.
             let ratio = 2.0f64.powi(evolution.octave as i32);
             // The scale of the edge filter.
             let sigma_size = f64::round(evolution.esigma * self.derivative_factor / ratio) as u32;
             compute_multiscale_derivatives_for_evolution(evolution, sigma_size);
+        };
+        #[cfg(not(feature = "rayon"))]
+        for evolution in evolutions.iter_mut() {
+            process_evolution(evolution);
         }
+        #[cfg(feature = "rayon")]
+        evolutions.into_par_iter().for_each(|evolution| {
+            process_evolution(evolution);
+        });
     }
 
     /// Compute the detector response - the determinant of the Hessian - and save the result
@@ -21,7 +32,7 @@ impl Akaze {
     #[allow(non_snake_case, clippy::suspicious_operation_groupings)]
     pub fn detector_response(&self, evolutions: &mut [EvolutionStep]) {
         self.compute_multiscale_derivatives(evolutions);
-        for evolution in evolutions.iter_mut() {
+        let process_evolution = |evolution: &mut EvolutionStep| {
             let ratio = f64::powi(2.0, evolution.octave as i32);
             let sigma_size = f64::round(evolution.esigma * self.derivative_factor / ratio);
             let sigma_size_quat = sigma_size.powi(4) as f32;
@@ -34,14 +45,41 @@ impl Akaze {
             ) {
                 *Ldet = (Lxx * Lyy - Lxy * Lxy) * sigma_size_quat;
             });
+        };
+        #[cfg(not(feature = "rayon"))]
+        for evolution in evolutions.iter_mut() {
+            process_evolution(evolution);
         }
+        #[cfg(feature = "rayon")]
+        evolutions.into_par_iter().for_each(|evolution| {
+            process_evolution(evolution);
+        });
     }
 }
 
 fn compute_multiscale_derivatives_for_evolution(evolution: &mut EvolutionStep, sigma_size: u32) {
-    evolution.Lx = derivatives::scharr_horizontal(&evolution.Lsmooth, sigma_size);
-    evolution.Ly = derivatives::scharr_vertical(&evolution.Lsmooth, sigma_size);
-    evolution.Lxx = derivatives::scharr_horizontal(&evolution.Lx, sigma_size);
-    evolution.Lyy = derivatives::scharr_vertical(&evolution.Ly, sigma_size);
-    evolution.Lxy = derivatives::scharr_vertical(&evolution.Lx, sigma_size);
+    #[cfg(not(feature = "rayon"))]
+    {
+        evolution.Lx = derivatives::scharr_horizontal(&evolution.Lsmooth, sigma_size);
+        evolution.Ly = derivatives::scharr_vertical(&evolution.Lsmooth, sigma_size);
+        evolution.Lxx = derivatives::scharr_horizontal(&evolution.Lx, sigma_size);
+        evolution.Lyy = derivatives::scharr_vertical(&evolution.Ly, sigma_size);
+        evolution.Lxy = derivatives::scharr_vertical(&evolution.Lx, sigma_size);
+    }
+    #[cfg(feature = "rayon")]
+    {
+        (evolution.Lx, evolution.Ly) = rayon::join(
+            || derivatives::scharr_horizontal(&evolution.Lsmooth, sigma_size),
+            || derivatives::scharr_vertical(&evolution.Lsmooth, sigma_size),
+        );
+        (evolution.Lxx, (evolution.Lyy, evolution.Lxy)) = rayon::join(
+            || derivatives::scharr_horizontal(&evolution.Lx, sigma_size),
+            || {
+                rayon::join(
+                    || derivatives::scharr_vertical(&evolution.Ly, sigma_size),
+                    || derivatives::scharr_vertical(&evolution.Lx, sigma_size),
+                )
+            },
+        )
+    }
 }
diff --git a/akaze/src/image.rs b/akaze/src/image.rs
index 5eaf6e7f..7151749e 100644
--- a/akaze/src/image.rs
+++ b/akaze/src/image.rs
@@ -1,10 +1,12 @@
 use derive_more::{Deref, DerefMut};
 use image::{DynamicImage, ImageBuffer, Luma, Pixel};
-use imageproc::filter::separable_filter_equal;
 use log::*;
 use ndarray::{azip, s, Array2, ArrayView2, ArrayViewMut2};
 use nshare::{MutNdarray2, RefNdarray2};
 use std::f32;
+use wide::f32x4;
+
+type GrayImageBuffer = ImageBuffer<Luma<f32>, Vec<f32>>;
 
 /// The image type we use in this library.
 ///
@@ -31,7 +33,7 @@ use std::f32;
 /// like using a separable filter, and using the filters implemented
 /// here ended up speeding up everything a lot.
 #[derive(Debug, Clone, Deref, DerefMut)]
-pub struct GrayFloatImage(pub ImageBuffer<Luma<f32>, Vec<f32>>);
+pub struct GrayFloatImage(pub GrayImageBuffer);
 
 impl GrayFloatImage {
     /// Create a unit float image from the image crate's DynamicImage type.
@@ -177,6 +179,146 @@ impl GrayFloatImage {
     }
 }
 
+pub fn horizontal_filter(image: &GrayImageBuffer, kernel: &[f32]) -> GrayImageBuffer {
+    // Validate kernel size.
+    let kernel_size = kernel.len();
+    debug_assert!(kernel_size % 2 == 1);
+    let kernel_half_size = kernel_size / 2;
+    // Prepare output.
+    let width = image.width() as usize;
+    let height = image.height() as usize;
+    let mut output = vec![0.0; width * height];
+    // Create SIMD kernel, padded with 0.
+    let kernel_simd = kernel
+        .chunks(4)
+        .map(|chunk| {
+            let data = [
+                #[allow(clippy::get_first)]
+                chunk.get(0).copied().unwrap_or(0.0),
+                chunk.get(1).copied().unwrap_or(0.0),
+                chunk.get(2).copied().unwrap_or(0.0),
+                chunk.get(3).copied().unwrap_or(0.0),
+            ];
+            f32x4::new(data)
+        })
+        .collect::<Vec<_>>();
+    let kernel_simd_size = 4 * (kernel_size + 3) / 4;
+    let kernel_simd_extra_elements = kernel_simd_size - kernel_size;
+    // Process each row independently.
+    let row_in_it = image.as_raw().chunks_exact(width);
+    let row_out_it = output.chunks_exact_mut(width);
+    let mut scratch = vec![0f32; width + kernel_half_size * 2 + kernel_simd_extra_elements];
+    for (row_in, row_out) in row_in_it.zip(row_out_it) {
+        // Prefill extended buffer with center and edge values.
+        scratch[0..kernel_half_size].fill(row_in[0]);
+        scratch[kernel_half_size..kernel_half_size + width].copy_from_slice(row_in);
+        scratch[kernel_half_size + width..2 * kernel_half_size + width].fill(row_in[width - 1]);
+        scratch[2 * kernel_half_size + width..].fill(0.);
+        // Apply kernel.
+        scratch
+            .windows(kernel_simd_size)
+            .zip(row_out)
+            .for_each(|(window, output)| {
+                *output = window
+                    .chunks_exact(4)
+                    .map(|chunk| f32x4::new([chunk[0], chunk[1], chunk[2], chunk[3]]))
+                    .zip(kernel_simd.iter())
+                    .fold(f32x4::splat(0.), |acc, (a, b)| a.mul_add(*b, acc))
+                    .reduce_add()
+            });
+    }
+    GrayImageBuffer::from_raw(width as u32, height as u32, output).unwrap()
+}
+
+pub fn vertical_filter(image: &GrayImageBuffer, kernel: &[f32]) -> GrayImageBuffer {
+    // Validate kernel size.
+    let kernel_size = kernel.len();
+    debug_assert!(kernel_size % 2 == 1);
+    let kernel_half_size = kernel_size / 2;
+    // Prepare output.
+    let width = image.width() as usize;
+    let height = image.height() as usize;
+    let mut output = vec![0.0; width * height];
+    // Create SIMD kernel, padded with 0.
+    let kernel_simd = kernel
+        .chunks(4)
+        .map(|chunk| {
+            let data = [
+                #[allow(clippy::get_first)]
+                chunk.get(0).copied().unwrap_or(0.0),
+                chunk.get(1).copied().unwrap_or(0.0),
+                chunk.get(2).copied().unwrap_or(0.0),
+                chunk.get(3).copied().unwrap_or(0.0),
+            ];
+            f32x4::new(data)
+        })
+        .collect::<Vec<_>>();
+    let kernel_simd_size = 4 * (kernel_size + 3) / 4;
+    let kernel_simd_extra_elements = kernel_simd_size - kernel_size;
+    // We use a scratch buffer of L1 cache width (64 bytes) to optimize memory access.
+    const SCRATCH_WIDTH: usize = 16;
+    let scratch_height = height + kernel_half_size * 2 + kernel_simd_extra_elements;
+    let mut scratch = vec![0f32; SCRATCH_WIDTH * scratch_height];
+    let image = image.as_raw();
+    for x_s in (0..width).step_by(SCRATCH_WIDTH) {
+        // Fill the scratch buffer with the column values.
+        // First paddings.
+        let x_e: usize = (x_s + SCRATCH_WIDTH).min(width);
+        for x in x_s..x_e {
+            let scratch_col_start = (x - x_s) * scratch_height;
+            for i in 0..kernel_half_size {
+                scratch[scratch_col_start + i] = image[x];
+            }
+            let scratch_end_col_start = scratch_col_start + kernel_half_size + height;
+            let image_last_row_start = (height - 1) * width;
+            for i in 0..kernel_half_size {
+                scratch[scratch_end_col_start + i] = image[image_last_row_start + x];
+            }
+            for i in 0..kernel_simd_extra_elements {
+                scratch[scratch_end_col_start + kernel_half_size + i] = 0.;
+            }
+        }
+        // Then main content.
+        for y in 0..height {
+            let image_row_start = y * width;
+            for x in x_s..x_e {
+                scratch[(x - x_s) * scratch_height + y + kernel_half_size] =
+                    image[image_row_start + x];
+            }
+        }
+        // Apply kernel.
+        let col_count = x_e - x_s;
+        scratch
+            .chunks(scratch_height)
+            .take(col_count)
+            .enumerate()
+            .for_each(|(dx, col)| {
+                let x = x_s + dx;
+                col.windows(kernel_simd_size)
+                    .enumerate()
+                    .for_each(|(y, window)| {
+                        let value = window
+                            .chunks_exact(4)
+                            .map(|chunk| f32x4::new([chunk[0], chunk[1], chunk[2], chunk[3]]))
+                            .zip(kernel_simd.iter())
+                            .fold(f32x4::splat(0.), |acc, (a, b)| a.mul_add(*b, acc))
+                            .reduce_add();
+                        output[y * width + x] = value;
+                    });
+            });
+    }
+    GrayImageBuffer::from_raw(width as u32, height as u32, output).unwrap()
+}
+
+pub fn separable_filter(
+    image: &GrayImageBuffer,
+    h_kernel: &[f32],
+    v_kernel: &[f32],
+) -> GrayImageBuffer {
+    let h = horizontal_filter(image, h_kernel);
+    vertical_filter(&h, v_kernel)
+}
+
 /// The Gaussian function.
 ///
 /// # Arguments
@@ -188,14 +330,15 @@ fn gaussian(x: f32, r: f32) -> f32 {
     ((2.0 * f32::consts::PI).sqrt() * r).recip() * (-x.powi(2) / (2.0 * r.powi(2))).exp()
 }
 
-/// Generate a Gaussina kernel.
+/// Generate a Gaussian kernel.
 ///
 /// # Arguments
 /// * `r` - sigma.
 /// * `kernel_size` - The size of the kernel.
 /// # Return value
 /// The kernel (a vector).
-fn gaussian_kernel(r: f32, kernel_size: usize) -> Vec<f32> {
+pub fn gaussian_kernel(r: f32, kernel_size: usize) -> Vec<f32> {
+    assert!(kernel_size % 2 == 1, "kernel_size must be odd");
     let mut kernel = vec![0f32; kernel_size];
     let half_width = (kernel_size / 2) as i32;
     let mut sum = 0f32;
@@ -222,12 +365,13 @@ pub fn gaussian_blur(image: &GrayFloatImage, r: f32) -> GrayFloatImage {
     let kernel_radius = (2.0 * r).ceil() as usize;
     let kernel_size = kernel_radius * 2 + 1;
     let kernel = gaussian_kernel(r, kernel_size);
-    GrayFloatImage(separable_filter_equal(image, &kernel))
+    GrayFloatImage(separable_filter(image, &kernel, &kernel))
 }
 
 #[cfg(test)]
 mod tests {
     use super::gaussian_kernel;
+
     #[test]
     fn gaussian_kernel_correct() {
         // test against known correct kernel
@@ -246,4 +390,24 @@ mod tests {
             assert!(f32::abs(*i - *j) < 0.0001);
         }
     }
+
+    #[test]
+    fn horizontal_filter() {
+        let image = image::open("../res/0000000000.png").unwrap();
+        let image = super::GrayFloatImage::from_dynamic(&image);
+        let kernel = gaussian_kernel(3.0, 7);
+        let filtered_ours = super::horizontal_filter(&image.0, &kernel);
+        let filtered_imageproc = imageproc::filter::horizontal_filter(&image.0, &kernel);
+        imageproc::assert_pixels_eq_within!(filtered_ours, filtered_imageproc, 0.0001);
+    }
+
+    #[test]
+    fn vertical_filter() {
+        let image = image::open("../res/0000000000.png").unwrap();
+        let image = super::GrayFloatImage::from_dynamic(&image);
+        let kernel = gaussian_kernel(3.0, 7);
+        let filtered_ours = super::vertical_filter(&image.0, &kernel);
+        let filtered_imageproc = imageproc::filter::vertical_filter(&image.0, &kernel);
+        imageproc::assert_pixels_eq_within!(filtered_ours, filtered_imageproc, 0.0001);
+    }
 }
diff --git a/akaze/src/lib.rs b/akaze/src/lib.rs
index 659b0b21..de9b1043 100644
--- a/akaze/src/lib.rs
+++ b/akaze/src/lib.rs
@@ -4,7 +4,7 @@ mod descriptors;
 mod detector_response;
 mod evolution;
 mod fed_tau;
-mod image;
+pub mod image;
 mod nonlinear_diffusion;
 mod scale_space_extrema;
 
@@ -18,6 +18,36 @@ use log::*;
 use nonlinear_diffusion::pm_g2;
 use std::{cmp::Reverse, path::Path};
 
+#[cfg(not(target_arch = "wasm32"))]
+use std::time::Instant;
+#[cfg(target_arch = "wasm32")]
+struct Instant(f64);
+#[cfg(target_arch = "wasm32")]
+impl Instant {
+    fn now() -> Self {
+        #[cfg(feature = "web-sys")]
+        {
+            Self(web_sys::window().unwrap().performance().unwrap().now())
+        }
+        #[cfg(not(feature = "web-sys"))]
+        {
+            Self(0.)
+        }
+    }
+    fn elapsed(&self) -> std::time::Duration {
+        #[cfg(feature = "web-sys")]
+        {
+            std::time::Duration::from_secs_f64(
+                (web_sys::window().unwrap().performance().unwrap().now() - self.0) * 0.001,
+            )
+        }
+        #[cfg(not(feature = "web-sys"))]
+        {
+            std::time::Duration::from_secs(0)
+        }
+    }
+}
+
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
 
@@ -199,16 +229,27 @@ impl Akaze {
             }
             evolutions[i].Lsmooth = gaussian_blur(&evolutions[i].Lt, 1.0f32);
             trace!("Gaussian blur finished.");
-            evolutions[i].Lx = derivatives::simple_scharr_horizontal(&evolutions[i].Lsmooth);
-            trace!("Computing derivative Lx done.");
-            evolutions[i].Ly = derivatives::simple_scharr_vertical(&evolutions[i].Lsmooth);
-            trace!("Computing derivative Ly done.");
+            #[cfg(not(feature = "rayon"))]
+            {
+                evolutions[i].Lx = derivatives::simple_scharr_horizontal(&evolutions[i].Lsmooth);
+                trace!("Computing derivative Lx done.");
+                evolutions[i].Ly = derivatives::simple_scharr_vertical(&evolutions[i].Lsmooth);
+                trace!("Computing derivative Ly done.");
+            }
+            #[cfg(feature = "rayon")]
+            {
+                (evolutions[i].Lx, evolutions[i].Ly) = rayon::join(
+                    || derivatives::simple_scharr_horizontal(&evolutions[i].Lsmooth),
+                    || derivatives::simple_scharr_vertical(&evolutions[i].Lsmooth),
+                );
+            }
             evolutions[i].Lflow = pm_g2(&evolutions[i].Lx, &evolutions[i].Ly, contrast_factor);
             trace!("Lflow finished.");
-            for j in 0..evolutions[i].fed_tau_steps.len() {
+            let evolution = &mut evolutions[i];
+            for j in 0..evolution.fed_tau_steps.len() {
                 trace!("Starting diffusion step.");
-                let step_size = evolutions[i].fed_tau_steps[j];
-                nonlinear_diffusion::calculate_step(&mut evolutions[i], step_size as f32);
+                let step_size = evolution.fed_tau_steps[j];
+                nonlinear_diffusion::calculate_step(evolution, step_size as f32);
                 trace!("Diffusion step finished with step size {}", step_size);
             }
         }
@@ -223,9 +264,13 @@ impl Akaze {
     /// The resulting keypoints.
     ///
     pub fn find_image_keypoints(&self, evolutions: &mut [EvolutionStep]) -> Vec<KeyPoint> {
+        let start = Instant::now();
         self.detector_response(evolutions);
-        trace!("Computing detector response finished.");
-        self.detect_keypoints(evolutions)
+        info!("Computed detector response in: {:?}", start.elapsed());
+        let start = Instant::now();
+        let keypoints = self.detect_keypoints(evolutions);
+        info!("Detected keypoints in: {:?}", start.elapsed());
+        keypoints
     }
 
     /// Extract features using the Akaze feature extractor.
@@ -263,18 +308,31 @@ impl Akaze {
         &self,
         float_image: &GrayFloatImage,
     ) -> (Vec<KeyPoint>, Vec<BitArray<64>>) {
+        trace!("Allocating evolutions.");
+        let start = Instant::now();
         let mut evolutions =
             self.allocate_evolutions(float_image.0.width(), float_image.0.height());
+        info!("Allocated evolutions in: {:?}", start.elapsed());
+        trace!("Creating non-linear space.");
+        let start = Instant::now();
         self.create_nonlinear_scale_space(&mut evolutions, float_image);
+        info!("Created non-linear scale space in: {:?}", start.elapsed());
         trace!("Finding image keypoints.");
         let mut keypoints = self.find_image_keypoints(&mut evolutions);
         trace!("Sorting keypoints by response and truncating the worst keypoints based on the set maximum");
+        let start = Instant::now();
         keypoints.sort_unstable_by_key(|kp| Reverse(FloatOrd(kp.response)));
         keypoints.truncate(self.maximum_features);
+        info!(
+            "Keypoints sorted in: {:?}, {} left.",
+            start.elapsed(),
+            keypoints.len()
+        );
         trace!("Extracting descriptors.");
+        let start = Instant::now();
         let (keypoints, descriptors) = self.extract_descriptors(&evolutions, &keypoints);
+        info!("Extracted keypoints in: {:?}", start.elapsed());
         trace!("Computing descriptors finished.");
-        info!("Extracted {} features", keypoints.len());
         (keypoints, descriptors)
     }
 
diff --git a/akaze/src/scale_space_extrema.rs b/akaze/src/scale_space_extrema.rs
index f863ddc9..f368baf4 100644
--- a/akaze/src/scale_space_extrema.rs
+++ b/akaze/src/scale_space_extrema.rs
@@ -2,6 +2,9 @@ use crate::{evolution::EvolutionStep, Akaze, KeyPoint};
 use log::*;
 use std::f32::consts::PI;
 
+#[cfg(feature = "rayon")]
+use rayon::prelude::*;
+
 impl Akaze {
     /// Compute scale space extrema to get the detector response.
     ///
@@ -11,7 +14,7 @@ impl Akaze {
     fn find_scale_space_extrema(&self, evolutions: &mut [EvolutionStep]) -> Vec<KeyPoint> {
         let mut keypoint_cache: Vec<KeyPoint> = vec![];
         let smax = 10.0f32 * f32::sqrt(2.0f32);
-        for (e_id, evolution) in evolutions.iter_mut().enumerate() {
+        for (e_id, evolution) in evolutions.iter().enumerate() {
             let w = evolution.Ldet.width();
             let h = evolution.Ldet.height();
             // maintain 9 iterators, one for the current pixel and one
@@ -239,12 +242,12 @@ fn compute_main_orientation(keypoint: &mut KeyPoint, evolutions: &[EvolutionStep
     let cv_fast_atan2_equiv = |y: f32, x: f32| (y.atan2(x) + 2. * PI).rem_euclid(2. * PI);
     // Calculate derivatives responses for points within radius of 6*scale
     let mut idx = 0;
-    for i in -6..=6 {
-        for j in -6..=6 {
+    for j in -6..=6 {
+        for i in -6..=6 {
             if i * i + j * j < 36 {
                 let iy = f32::round(yf + (j as f32) * s) as usize;
                 let ix = f32::round(xf + (i as f32) * s) as usize;
-                let gweight = GAUSS25[id[(i + 6) as usize]][id[(j + 6) as usize]];
+                let gweight = GAUSS25[id[(j + 6) as usize]][id[(i + 6) as usize]];
                 res_x[idx] = gweight * evolutions[level].Lx.get(ix, iy);
                 res_y[idx] = gweight * evolutions[level].Ly.get(ix, iy);
                 angs[idx] = cv_fast_atan2_equiv(res_y[idx], res_x[idx]);
@@ -295,8 +298,7 @@ fn do_subpixel_refinement(
     in_keypoints: &[KeyPoint],
     evolutions: &[EvolutionStep],
 ) -> Vec<KeyPoint> {
-    let mut result: Vec<KeyPoint> = vec![];
-    for keypoint in in_keypoints.iter() {
+    let process_keypoint = |keypoint: &KeyPoint| {
         let ratio = f32::powf(2.0f32, keypoint.octave as f32);
         let x = f32::round(keypoint.point.0 / ratio) as usize;
         let y = f32::round(keypoint.point.1 / ratio) as usize;
@@ -338,16 +340,23 @@ fn do_subpixel_refinement(
             );
             assert_eq!(keypoint_clone.angle, 0.);
             keypoint_clone.size *= 2.;
-            result.push(keypoint_clone);
+            compute_main_orientation(&mut keypoint_clone, evolutions);
+            Some(keypoint_clone)
+        } else {
+            None
         }
-    }
+    };
+    #[cfg(not(feature = "rayon"))]
+    let result: Vec<_> = in_keypoints.iter().filter_map(process_keypoint).collect();
+    #[cfg(feature = "rayon")]
+    let result: Vec<_> = in_keypoints
+        .par_iter()
+        .filter_map(process_keypoint)
+        .collect();
     debug!(
         "{}/{} remain after subpixel refinement.",
         result.len(),
         in_keypoints.len()
     );
-    for keypoint in result.iter_mut() {
-        compute_main_orientation(keypoint, evolutions);
-    }
     result
 }
diff --git a/cv-sfm/Cargo.toml b/cv-sfm/Cargo.toml
index e03963bf..3606fafd 100644
--- a/cv-sfm/Cargo.toml
+++ b/cv-sfm/Cargo.toml
@@ -8,12 +8,12 @@ edition = "2021"
 serde-serialize = [
     "serde",
     "cv-core/serde-serialize",
-    "bitarray/serde",
     "cv-pinhole/serde-serialize",
     "bitarray/serde",
     "slotmap/serde",
     "hgg/serde",
     "hamming-lsh/serde",
+    "akaze/serde",
 ]
 
 [dependencies]


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 4382afd1..026df7e4 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -35,4 +35,10 @@ jobs:
         uses: actions-rs/cargo@v1
         with:
           command: test
-          args: --features cv/serde-serialize
\ No newline at end of file
+          args: --features cv/serde-serialize
+      
+      - name: Run cargo test for Akaze with rayon feature
+        uses: actions-rs/cargo@v1
+        with:
+          command: test
+          args: --features akaze/rayon
\ No newline at end of file