Fix akaze orientation (from upstream branch)

akaze/Cargo.toml

@@ -25,6 +25,8 @@ ndarray = { version = "0.15.4", default-features = false }
 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 }
+imageproc = "0.23.0"
 
 
 [dev-dependencies]

akaze/src/contrast_factor.rs

@@ -22,8 +22,8 @@ pub fn compute_contrast_factor(
     let mut num_points: f64 = 0.0;
     let mut histogram = vec![0; num_bins];
     let gaussian = gaussian_blur(image, gradient_histogram_scale as f32);
-    let Lx = crate::derivatives::scharr_horizontal(&gaussian, 1);
-    let Ly = crate::derivatives::scharr_vertical(&gaussian, 1);
+    let Lx = crate::derivatives::simple_scharr_horizontal(&gaussian);
+    let Ly = crate::derivatives::simple_scharr_vertical(&gaussian);
     let hmax = (1..gaussian.height() - 1)
         .flat_map(|y| (1..gaussian.width() - 1).map(move |x| (x, y)))
         .map(|(x, y)| Lx.get(x, y).powi(2) as f64 + Ly.get(x, y).powi(2) as f64)

akaze/src/derivatives.rs

@@ -1,5 +1,22 @@
-use crate::image::{fill_border, GrayFloatImage};
-use ndarray::{s, Array2, ArrayView2, ArrayViewMut2};
+use crate::image::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.],
+    ))
+}
+
+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.],
+    ))
+}
 
 /// Compute the Scharr derivative horizontally
 ///
@@ -12,18 +29,16 @@ use ndarray::{s, Array2, ArrayView2, ArrayViewMut2};
 /// # Return value
 /// Output image derivative (an image.)
 pub fn scharr_horizontal(image: &GrayFloatImage, sigma_size: u32) -> GrayFloatImage {
-    let img_horizontal = scharr_axis(
-        image,
-        sigma_size,
-        FilterDirection::Horizontal,
-        FilterOrder::Main,
-    );
-    scharr_axis(
-        &img_horizontal,
-        sigma_size,
-        FilterDirection::Vertical,
-        FilterOrder::Off,
-    )
+    if sigma_size == 1 {
+        return simple_scharr_horizontal(image);
+    }
+    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,
+    ))
 }
 
 /// Compute the Scharr derivative vertically
@@ -37,46 +52,16 @@ pub fn scharr_horizontal(image: &GrayFloatImage, sigma_size: u32) -> GrayFloatIm
 /// # Return value
 /// Output image derivative (an image.)
 pub fn scharr_vertical(image: &GrayFloatImage, sigma_size: u32) -> GrayFloatImage {
-    let img_horizontal = scharr_axis(
-        image,
-        sigma_size,
-        FilterDirection::Horizontal,
-        FilterOrder::Off,
-    );
-    scharr_axis(
-        &img_horizontal,
-        sigma_size,
-        FilterDirection::Vertical,
-        FilterOrder::Main,
-    )
-}
-
-/// Multiplies and accumulates
-fn accumulate_mul_offset(
-    mut accumulator: ArrayViewMut2<f32>,
-    source: ArrayView2<f32>,
-    val: f32,
-    border: usize,
-    xoff: usize,
-    yoff: usize,
-) {
-    assert_eq!(source.dim(), accumulator.dim());
-    let dims = source.dim();
-    let mut accumulator =
-        accumulator.slice_mut(s![border..dims.0 - border, border..dims.1 - border]);
-    accumulator.scaled_add(
-        val,
-        &source.slice(s![
-            yoff..dims.0 + yoff - 2 * border,
-            xoff..dims.1 + xoff - 2 * border
-        ]),
-    );
-}
-
-#[derive(Copy, Clone, Debug, PartialEq)]
-enum FilterDirection {
-    Horizontal,
-    Vertical,
+    if sigma_size == 1 {
+        return simple_scharr_vertical(image);
+    }
+    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,
+    ))
 }
 
 #[derive(Copy, Clone, Debug, PartialEq)]
@@ -85,43 +70,26 @@ enum FilterOrder {
     Off,
 }
 
-fn scharr_axis(
-    image: &GrayFloatImage,
-    sigma_size: u32,
-    dir: FilterDirection,
-    order: FilterOrder,
-) -> GrayFloatImage {
-    let mut output = Array2::<f32>::zeros([image.height(), image.width()]);
-    // Get the border size (we wont fill in this border width of the output).
-    let border = sigma_size as usize;
+fn computer_scharr_kernel(sigma_size: u32, order: FilterOrder) -> Vec<f32> {
     // Difference between middle and sides of main axis filter.
     let w = 10.0 / 3.0;
     // Side intensity of filter.
     let norm = (1.0 / (2.0 * f64::from(sigma_size) * (w + 2.0))) as f32;
     // Middle intensity of filter.
     let middle = norm * w as f32;
-
-    let mut offsets = match order {
-        FilterOrder::Main => vec![
-            (norm, [border, 0]),
-            (middle, [border, border]),
-            (norm, [border, 2 * border]),
-        ],
-        FilterOrder::Off => vec![(-1.0, [border, 0]), (1.0, [border, 2 * border])],
-    };
-
-    if dir == FilterDirection::Horizontal {
-        // Swap the offsets if the filter is a horizontal filter.
-        for (_, [x, y]) in &mut offsets {
-            std::mem::swap(x, y);
+    // Size of kernel
+    let ksize = (3 + 2 * (sigma_size - 1)) as usize;
+    let mut kernel = vec![0.0; ksize];
+    match order {
+        FilterOrder::Main => {
+            kernel[0] = -1.0;
+            kernel[ksize - 1] = 1.0;
         }
-    }
-
-    // Accumulate the three components.
-    for (val, [x, y]) in offsets {
-        accumulate_mul_offset(output.view_mut(), image.ref_array2(), val, border, x, y);
-    }
-    let mut output = GrayFloatImage::from_array2(output);
-    fill_border(&mut output, border);
-    output
+        FilterOrder::Off => {
+            kernel[0] = norm;
+            kernel[ksize / 2] = middle;
+            kernel[ksize - 1] = norm;
+        }
+    };
+    kernel
 }

akaze/src/descriptors.rs

@@ -1,4 +1,4 @@
-use crate::{Akaze, EvolutionStep, KeyPoint};
+use crate::{Akaze, Error, EvolutionStep, KeyPoint};
 use bitarray::BitArray;
 
 impl Akaze {
@@ -14,11 +14,16 @@ impl Akaze {
         &self,
         evolutions: &[EvolutionStep],
         keypoints: &[KeyPoint],
-    ) -> Vec<BitArray<64>> {
+    ) -> (Vec<KeyPoint>, Vec<BitArray<64>>) {
         keypoints
             .iter()
-            .map(|keypoint| self.get_mldb_descriptor(keypoint, evolutions))
-            .collect()
+            .filter_map(|&keypoint| {
+                Some((
+                    keypoint,
+                    self.get_mldb_descriptor(&keypoint, evolutions).ok()?,
+                ))
+            })
+            .unzip()
     }
 
     /// Computes the rotation invariant M-LDB binary descriptor (maximum descriptor length)
@@ -33,20 +38,22 @@ impl Akaze {
         &self,
         keypoint: &KeyPoint,
         evolutions: &[EvolutionStep],
-    ) -> BitArray<64> {
+    ) -> Result<BitArray<64>, Error> {
         let mut output = BitArray::zeros();
         let max_channels = 3usize;
         debug_assert!(self.descriptor_channels <= max_channels);
         let mut values: Vec<f32> = vec![0f32; 16 * max_channels];
         let size_mult = [1.0f32, 2.0f32 / 3.0f32, 1.0f32 / 2.0f32];
+
         let ratio = (1u32 << keypoint.octave) as f32;
         let scale = f32::round(0.5f32 * keypoint.size / ratio);
         let xf = keypoint.point.0 / ratio;
         let yf = keypoint.point.1 / ratio;
         let co = f32::cos(keypoint.angle);
         let si = f32::sin(keypoint.angle);
-        let mut dpos = 0usize;
         let pattern_size = self.descriptor_pattern_size as f32;
+
+        let mut dpos = 0usize;
         for (lvl, multiplier) in size_mult.iter().enumerate() {
             let val_count = (lvl + 2usize) * (lvl + 2usize);
             let sample_size = f32::ceil(pattern_size * multiplier) as usize;
@@ -60,7 +67,7 @@ impl Akaze {
                 si,
                 scale,
                 evolutions,
-            );
+            )?;
             mldb_binary_comparisons(
                 &values,
                 output.bytes_mut(),
@@ -69,7 +76,7 @@ impl Akaze {
                 self.descriptor_channels,
             );
         }
-        output
+        Ok(output)
     }
 
     /// Fill the comparison values for the MLDB rotation invariant descriptor
@@ -85,7 +92,7 @@ impl Akaze {
         si: f32,
         scale: f32,
         evolutions: &[EvolutionStep],
-    ) {
+    ) -> Result<(), Error> {
         let pattern_size = self.descriptor_pattern_size;
         let nr_channels = self.descriptor_channels;
         let mut valuepos = 0;
@@ -97,12 +104,22 @@ impl Akaze {
                 let mut nsamples = 0usize;
                 for k in i..(i + (sample_step as i32)) {
                     for l in j..(j + (sample_step as i32)) {
-                        let l = l as f32 + 0.5;
-                        let k = k as f32 + 0.5;
+                        let l = l as f32;
+                        let k = k as f32;
                         let sample_y = yf + (l * co * scale + k * si * scale);
                         let sample_x = xf + (-l * si * scale + k * co * scale);
                         let y1 = f32::round(sample_y) as isize;
                         let x1 = f32::round(sample_x) as isize;
+                        if !(0..evolutions[level].Lt.width() as isize).contains(&x1)
+                            || !(0..evolutions[level].Lt.height() as isize).contains(&y1)
+                        {
+                            return Err(Error::SampleOutOfBounds {
+                                x: x1,
+                                y: y1,
+                                width: evolutions[level].Lt.width(),
+                                height: evolutions[level].Lt.height(),
+                            });
+                        }
                         let y1 = y1 as usize;
                         let x1 = x1 as usize;
                         let ri = evolutions[level].Lt.get(x1, y1);
@@ -138,6 +155,7 @@ impl Akaze {
                 valuepos += nr_channels;
             }
         }
+        Ok(())
     }
 }