implemented gftt_response

crates/kornia-imgproc/Cargo.toml

@@ -14,6 +14,7 @@ version.workspace = true
 [dependencies]
 fast_image_resize = "5.1.0"
 kornia-tensor = { workspace = true }
+kornia-tensor-ops = { workspace = true }
 kornia-image = { workspace = true }
 num-traits = { workspace = true }
 rayon = "1.10"

crates/kornia-imgproc/benches/bench_features.rs

@@ -9,7 +9,7 @@ use rand::Rng;
 fn bench_fast_corner_detect(c: &mut Criterion) {
     let mut group = c.benchmark_group("FastCornerDetect");
 
-    let img_rgb8 = io::read_image_any_rgb8("/home/edgar/Downloads/kodim08_grayscale.png").unwrap();
+    let img_rgb8 = io::read_image_any_rgb8("../../tests/data/dog-rgb8.png").unwrap();
 
     let new_size = [1920, 1080].into();
     let mut img_resized = Image::from_size_val(new_size, 0).unwrap();
@@ -32,6 +32,38 @@ fn bench_fast_corner_detect(c: &mut Criterion) {
     );
 }
 
+fn bench_gftt_response(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Features");
+    let mut rng = rand::thread_rng();
+
+    for (width, height) in [(224, 224), (512, 512), (1920, 1080)].iter() {
+        group.throughput(criterion::Throughput::Elements((*width * *height) as u64));
+
+        let parameter_string = format!("{}x{}", width, height);
+
+        // input image
+        let image_data: Vec<f32> = (0..(*width * *height))
+            .map(|_| rng.gen_range(0.0..1.0))
+            .collect();
+        let image_size = [*width, *height].into();
+
+        let image_f32: Image<f32, 1> = Image::new(image_size, image_data).unwrap();
+
+        // output image
+        let response_f32: Image<f32, 1> = Image::from_size_val(image_size, 0.0).unwrap();
+
+        group.bench_with_input(
+            BenchmarkId::new("gftt", &parameter_string),
+            &(&image_f32, &response_f32),
+            |b, i| {
+                let (src, mut dst) = (i.0, i.1.clone());
+                b.iter(|| black_box(gftt_response(src, &mut dst)))
+            },
+        );
+    }
+    group.finish();
+}
+
 fn bench_harris_response(c: &mut Criterion) {
     let mut group = c.benchmark_group("Features");
     let mut rng = rand::thread_rng();
@@ -100,7 +132,7 @@ fn bench_dog_response(c: &mut Criterion) {
 criterion_group!(
     name = benches;
     config = Criterion::default().warm_up_time(std::time::Duration::new(10, 0));
-    targets = bench_harris_response, bench_dog_response, bench_fast_corner_detect
+    targets = bench_harris_response, bench_dog_response, bench_fast_corner_detect, bench_gftt_response
 );
 criterion_main!(benches);
 

crates/kornia-imgproc/src/features/responses.rs

@@ -1,5 +1,6 @@
-use crate::filter::gaussian_blur;
+use crate::filter::{gaussian_blur, spatial_gradient_float_parallel_row};
 use kornia_image::{Image, ImageError, ImageSize};
+use kornia_tensor_ops::TensorOps;
 use rayon::prelude::*;
 
 /// Method to calculate gradient for feature response
@@ -24,6 +25,64 @@ fn _get_kernel_size(sigma: f32) -> usize {
 
     ksize
 }
+///Compute the Shi-Tomasi cornerness function.
+///
+/// The Shi-Tomasi cornerness function is computed as the minimum eigenvalue of the gradient matrix.
+///
+/// Args:
+///    src: The source image.
+///   dst: The destination image.
+pub fn gftt_response(src: &Image<f32, 1>, dst: &mut Image<f32, 1>) -> Result<(), ImageError> {
+    if src.size() != dst.size() {
+        return Err(ImageError::InvalidImageSize(
+            src.cols(),
+            src.rows(),
+            dst.cols(),
+            dst.rows(),
+        ));
+    }
+    let size = src.size();
+    let mut dx = Image::from_size_val(size, 0.0)?;
+    let mut dy = Image::from_size_val(size, 0.0)?;
+
+    let mut dx2_g: Image<f32, 1> = Image::from_size_val(size, 0.0)?;
+    let mut dy2_g: Image<f32, 1> = Image::from_size_val(size, 0.0)?;
+    let mut dxy_g: Image<f32, 1> = Image::from_size_val(size, 0.0)?;
+    spatial_gradient_float_parallel_row(src, &mut dx, &mut dy)?;
+
+    let dx2 = Image(dx.mul(&dx).unwrap());
+    let dy2 = Image(dy.mul(&dy).unwrap());
+    let dxy = Image(dx.mul(&dy).unwrap());
+
+    gaussian_blur(&dx2, &mut dx2_g, (7, 7), (1.0, 1.0))?;
+    gaussian_blur(&dy2, &mut dy2_g, (7, 7), (1.0, 1.0))?;
+    gaussian_blur(&dxy, &mut dxy_g, (7, 7), (1.0, 1.0))?;
+    unsafe {
+        let det = dx2_g
+            .mul(&dy2_g)
+            .unwrap_unchecked()
+            .sub(&dxy_g.mul(&dxy_g).unwrap_unchecked())
+            .unwrap_unchecked();
+        let trace = dx2_g.add(&dy2_g).unwrap_unchecked();
+
+        let trace_sq = trace.mul(&trace).unwrap_unchecked();
+        let four_det = det.mul_scalar(4.0);
+        let discr = trace_sq.sub(&four_det).unwrap_unchecked().abs().powf(0.5);
+        let half = 0.5;
+        let e1 = trace.add(&discr).unwrap_unchecked().mul_scalar(half);
+        let e2 = trace.sub(&discr).unwrap_unchecked().mul_scalar(half);
+
+        let score = e1.min(&e2).unwrap_unchecked();
+
+        let sd = score.as_slice();
+        let dd = dst.as_slice_mut();
+        for i in 0..sd.len() {
+            dd[i] = sd[i];
+        }
+    }
+
+    Ok(())
+}
 
 /// Compute the Hessian response of an image.
 ///
@@ -325,7 +384,44 @@ pub fn dog_response(
 #[cfg(test)]
 mod tests {
     use super::*;
+    #[test]
+    fn test_gftt_response() -> Result<(), ImageError> {
+        #[rustfmt::skip]
+        let src = Image::from_size_slice(
+            [9, 9].into(),
+            &[
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0,
+                0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0,
+                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+                ],
+        )?;
+
+        let mut dst = Image::from_size_val([9, 9].into(), 0.0)?;
+        gftt_response(&src, &mut dst)?;
 
+        #[rustfmt::skip]
+        let expected_center_value = 0.1274;
+        assert!(
+            (dst.as_slice()[4 * 9 + 4] - expected_center_value).abs() < 1e-4,
+            "Center value should be close to expected value"
+        );
+        let max = dst
+            .as_slice()
+            .iter()
+            .max_by(|a, b| a.partial_cmp(b).unwrap())
+            .unwrap();
+        assert!(
+            (*max - expected_center_value).abs() < 1e-4,
+            "Max value should be close to centre value"
+        );
+        Ok(())
+    }
     #[test]
     fn test_hessian_response() -> Result<(), ImageError> {
         #[rustfmt::skip]

crates/kornia-imgproc/src/filter/ops.rs

@@ -220,37 +220,38 @@ pub fn spatial_gradient_float_parallel_row<const C: usize>(
     }
 
     let (sobel_x, sobel_y) = kernels::normalized_sobel_kernel3();
+    let rows = src.rows();
     let cols = src.cols();
 
     let src_data = src.as_slice();
 
-    dx.as_slice_mut()
+    let dx_data = dx.as_slice_mut();
+    let dy_data = dy.as_slice_mut();
+
+    dx_data
         .par_chunks_mut(cols * C)
-        .zip(dy.as_slice_mut().par_chunks_mut(cols * C))
+        .zip(dy_data.par_chunks_mut(cols * C))
         .enumerate()
         .for_each(|(r, (dx_row, dy_row))| {
-            dx_row
-                .chunks_mut(C)
-                .zip(dy_row.chunks_mut(C))
-                .enumerate()
-                .for_each(|(c, (dx_c, dy_c))| {
-                    let mut sum_x = [0.0; C];
-                    let mut sum_y = [0.0; C];
-                    for dy in 0..3 {
-                        for dx in 0..3 {
-                            let row = (r + dy).min(src.rows()).max(1) - 1;
-                            let col = (c + dx).min(src.cols()).max(1) - 1;
-                            for ch in 0..C {
-                                let src_pix_offset = (row * src.cols() + col) * C + ch;
-                                let val = unsafe { src_data.get_unchecked(src_pix_offset) };
-                                sum_x[ch] += val * sobel_x[dy][dx];
-                                sum_y[ch] += val * sobel_y[dy][dx];
-                            }
+            for c in 0..cols {
+                let mut sum_x = [0.0f32; C];
+                let mut sum_y = [0.0f32; C];
+                for ky in 0..3 {
+                    let y = (r + ky).saturating_sub(1).min(rows - 1);
+                    for kx in 0..3 {
+                        let x = (c + kx).saturating_sub(1).min(cols - 1);
+                        let src_idx = (y * cols + x) * C;
+                        for ch in 0..C {
+                            let val = unsafe { *src_data.get_unchecked(src_idx + ch) };
+                            sum_x[ch] += val * sobel_x[ky][kx];
+                            sum_y[ch] += val * sobel_y[ky][kx];
                         }
                     }
-                    dx_c.copy_from_slice(&sum_x);
-                    dy_c.copy_from_slice(&sum_y);
-                });
+                }
+                let out_idx = c * C;
+                dx_row[out_idx..out_idx + C].copy_from_slice(&sum_x);
+                dy_row[out_idx..out_idx + C].copy_from_slice(&sum_y);
+            }
         });
 
     Ok(())