diff --git a/.gitignore b/.gitignore index 90907854..7af95127 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ Cargo.lock # MacOS .DS_Store +.idea/ diff --git a/Cargo.toml b/Cargo.toml index cef8ef02..38bfc751 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,6 +7,7 @@ members = [ "cv-pinhole", "cv-optimize", "cv-sfm", + "cv-sift", "akaze", "eight-point", "lambda-twist", diff --git a/cv-sift/.cargo/config.toml b/cv-sift/.cargo/config.toml new file mode 100644 index 00000000..1fbbe44f --- /dev/null +++ b/cv-sift/.cargo/config.toml @@ -0,0 +1,2 @@ +[env] +RUST_LOG="trace" \ No newline at end of file diff --git a/cv-sift/.gitignore b/cv-sift/.gitignore new file mode 100644 index 00000000..830f5817 --- /dev/null +++ b/cv-sift/.gitignore @@ -0,0 +1 @@ +testdata \ No newline at end of file diff --git a/cv-sift/Cargo.toml b/cv-sift/Cargo.toml new file mode 100644 index 00000000..6f804c1a --- /dev/null +++ b/cv-sift/Cargo.toml @@ -0,0 +1,27 @@ +[package] +name = "cv-sift" +version = "0.1.0" +edition = "2018" +authors = ["Aalekh Patel "] +description = "SIFT feature extraction algorithm for computer vision" +keywords = ["keypoint", "descriptor", "vision", "sfm", "slam", "sift"] +categories = ["computer-vision", "science::robotics"] +repository = "https://github.com/rust-cv/cv" +documentation = "https://docs.rs/cv-sift/" +license = "MIT" +readme = "README.md" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html + +[dependencies] +# cv = { version = "0.6.0", path="../cv" } +image = { version = "0.24.6" } +nalgebra = "0.32.2" +num = "0.4" +itertools = "0.11.0" +tracing = "0.1.37" +thiserror = "1.0.40" + +[dev-dependencies] +test-case = "3.1.0" +tracing-subscriber = { version = "0.3.17", features = ["env-filter"] } diff --git a/cv-sift/LICENSE b/cv-sift/LICENSE new file mode 100644 index 00000000..fad12297 --- /dev/null +++ b/cv-sift/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2021 Aalekh Patel + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. diff --git a/cv-sift/README.md b/cv-sift/README.md new file mode 100644 index 00000000..e3eedd08 --- /dev/null +++ b/cv-sift/README.md @@ -0,0 +1,5 @@ +# SIFT + +Implementation of [SIFT](https://docs.opencv.org/master/da/df5/tutorial_py_sift_intro.html), one of the first breakthroughs in Computer Vision that introduced feature descriptors and keypoints to us common folks. + +## WIP diff --git a/cv-sift/src/config.rs b/cv-sift/src/config.rs new file mode 100644 index 00000000..05c8f2e8 --- /dev/null +++ b/cv-sift/src/config.rs @@ -0,0 +1,27 @@ + + +#[derive(Debug, Clone, Copy)] +pub struct SIFTConfig { + pub sigma: f64, + pub num_intervals: usize, + pub assumed_blur: f64, + pub image_border_width: usize +} + + +impl Default for SIFTConfig { + fn default() -> Self { + SIFTConfig { + sigma: 1.6, + num_intervals: 3, + assumed_blur: 0.5, + image_border_width: 5 + } + } +} + +impl SIFTConfig { + pub fn new() -> Self { + SIFTConfig::default() + } +} \ No newline at end of file diff --git a/cv-sift/src/conversion.rs b/cv-sift/src/conversion.rs new file mode 100644 index 00000000..25714984 --- /dev/null +++ b/cv-sift/src/conversion.rs @@ -0,0 +1,68 @@ +use image::{DynamicImage}; +use crate::ImageRgb32F; +use crate::errors::{Result, SIFTError}; + + +pub fn try_get_rgb_32f(img: &DynamicImage) -> Result { + match img { + + DynamicImage::ImageRgb8(p) + => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| pixel.0.iter().cloned().map(|v| v as f32)).collect::>() + ) + .unwrap()) + }, + DynamicImage::ImageRgb32F(p) + => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| pixel.0.iter().cloned()).collect::>() + ) + .unwrap()) + }, + DynamicImage::ImageRgba32F(p) + => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| pixel.0.iter().cloned()).collect::>() + ) + .unwrap()) + }, + DynamicImage::ImageLuma16(p) => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| [pixel.0[0] as f32; 3]).collect::>() + ).unwrap()) + } + DynamicImage::ImageLumaA16(p) => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| [pixel.0[0] as f32; 3]).collect::>() + ).unwrap()) + } + DynamicImage::ImageLuma8(p) => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| [pixel.0[0] as f32; 3]).collect::>() + ).unwrap()) + } + DynamicImage::ImageLumaA8(p) => { + Ok(ImageRgb32F::from_raw( + p.width(), + p.height(), + p.pixels().into_iter().flat_map(|pixel| [pixel.0[0] as f32; 3]).collect::>() + ).unwrap()) + } + _ => { + Err(SIFTError::Unsupported("ImageRgb32F image can only be one of ImageLuma8, ImageLumaA8, ImageLuma16, ImageLumaA16, ImageImageRgb32F, ImageRgba32F.".to_string())) + } + } +} diff --git a/cv-sift/src/errors.rs b/cv-sift/src/errors.rs new file mode 100644 index 00000000..7457f473 --- /dev/null +++ b/cv-sift/src/errors.rs @@ -0,0 +1,9 @@ +use thiserror::Error; + +#[derive(Debug, Error)] +pub enum SIFTError { + #[error("{0}")] + Unsupported(String), +} + +pub type Result = std::result::Result; \ No newline at end of file diff --git a/cv-sift/src/ext.rs b/cv-sift/src/ext.rs new file mode 100644 index 00000000..79ed8cb8 --- /dev/null +++ b/cv-sift/src/ext.rs @@ -0,0 +1,38 @@ +use crate::ImageRgb32F; +use image::Pixel; + +/// How close do we want two floats to be +/// before they're considered equal. +pub const TOLERANCE: f32 = 1e-8; + +pub trait ImageExt { + fn is_zero(&self) -> bool; + fn is_same_as(&self, other: &Self) -> bool; +} + +impl ImageExt for ImageRgb32F { + fn is_zero(&self) -> bool { + self + .pixels() + .all( + |pixel| + pixel + .channels() + .iter() + .find(|p| (**p).abs() > TOLERANCE) + .is_none() + ) + } + + fn is_same_as(&self, other: &Self) -> bool { + + for (own_pixel, other_pixel) in self.pixels().zip(other.pixels()) { + for (self_p, other_p) in own_pixel.channels().iter().zip(other_pixel.channels().iter()) { + if (*self_p - *other_p).abs() > TOLERANCE { + return false; + } + } + } + true + } +} diff --git a/cv-sift/src/lib.rs b/cv-sift/src/lib.rs new file mode 100644 index 00000000..02019159 --- /dev/null +++ b/cv-sift/src/lib.rs @@ -0,0 +1,12 @@ +pub mod config; +pub mod pyramid; +// Expose all utils. +pub mod utils; +pub mod conversion; +// #[cfg(test)] +pub mod ext; +mod errors; + +pub type ImageRgb32F = image::ImageBuffer, Vec>; + +pub use errors::*; \ No newline at end of file diff --git a/cv-sift/src/pyramid/base_image.rs b/cv-sift/src/pyramid/base_image.rs new file mode 100644 index 00000000..5ffc1049 --- /dev/null +++ b/cv-sift/src/pyramid/base_image.rs @@ -0,0 +1,107 @@ +use std::convert::TryInto; +use image::{DynamicImage, Rgb}; +use tracing::{ + trace, + debug +}; +use image::imageops; +use crate::ImageRgb32F; +use crate::errors::Result; +use image::Pixel; +use crate::conversion::try_get_rgb_32f; + + + +/// Scale the image by a factor of 2 and apply some blur. +/// # Examples +/// ``` +/// use cv_sift::pyramid::generate_base_image; +/// use image::{DynamicImage}; +/// +/// let img = image::open("tests/fixtures/box.png").unwrap(); +/// assert_eq!(img.height(), 223); +/// assert_eq!(img.width(), 324); +/// let base_img = generate_base_image(&img, 1.6, 0.5).unwrap(); +/// assert_eq!(base_img.height(), 446); +/// assert_eq!(base_img.width(), 648); +/// +/// ``` +pub fn generate_base_image( + img: &DynamicImage, + sigma: f64, + assumed_blur: f64 +) -> Result { + let (height, width) = (img.height(), img.width()); + + trace!( + height, + width, + sigma, + assumed_blur, + "Computing base image" + ); + + let scaled = img.resize( + width * 2, + height * 2, + imageops::FilterType::Triangle + ); + + trace!( + height = height * 2, + width = width * 2, + interpolation = "linear", + "Scaled image." + ); + + let final_sigma = { + let sigma_val = sigma * sigma - 4.0 * assumed_blur * assumed_blur; + if sigma_val > 0.01 { + sigma_val.sqrt() + } else { + 0.01_f64.sqrt() + } + }; + + debug!(final_sigma, "Computed final_sigma for blurring."); + + let blurred = scaled.blur(final_sigma as f32); + try_get_rgb_32f(&blurred) +} + + +/// Returns a difference of two images. +pub fn subtract(minuend: &ImageRgb32F, subtrahend: &ImageRgb32F) -> Result { + assert_eq!(minuend.height(), subtrahend.height()); + assert_eq!(minuend.width(), subtrahend.width()); + + let (width, height) = (minuend.width() as u32, minuend.height() as u32); + + let mut result_mat = ImageRgb32F::new(width, height); + let mut result_mat_pixels = result_mat.pixels_mut(); + + for (minuend_pixel, subtrahend_pixel) in minuend.pixels().zip(subtrahend.pixels()){ + + let output_pixel: [f32; 3] = + minuend_pixel + .channels() + .iter() + .zip( + subtrahend_pixel + .channels() + .iter() + ) + .map(|(minuend_p, subtrahend_p)| { + *minuend_p - *subtrahend_p + }) + .collect::>() + .try_into() + .unwrap(); + + let next_pixel = result_mat_pixels.next().unwrap(); + *next_pixel = Rgb(output_pixel); + } + Ok(result_mat) +} + + diff --git a/cv-sift/src/pyramid/kernels.rs b/cv-sift/src/pyramid/kernels.rs new file mode 100644 index 00000000..6f81539c --- /dev/null +++ b/cv-sift/src/pyramid/kernels.rs @@ -0,0 +1,36 @@ +/// List of gaussian kernels at which to blur the input image. +/// # Examples +/// ``` +/// use cv_sift::pyramid::gaussian_kernels; +/// use cv_sift::utils::assert_similar; +/// +/// let kernels = gaussian_kernels(1.6, 3); +/// let expected: [f64; 6] = [ +/// 1.6, +/// 1.2262735, +/// 1.54500779, +/// 1.94658784, +/// 2.452547, +/// 3.09001559 +/// ]; +/// assert_similar(&kernels, &expected); +/// ``` +pub fn gaussian_kernels( + sigma: f64, + num_intervals: usize +) -> Vec { + + let images_per_octave = num_intervals + 3; + let k: f64 = (2.0_f64).powf(1.0 / num_intervals as f64); + let mut kernels: Vec = vec![0.0; images_per_octave]; + + kernels[0] = sigma; + for (idx, item) in kernels.iter_mut().enumerate().take(images_per_octave).skip(1) { + let sigma_previous = (k.powf(idx as f64 - 1.0)) * sigma; + let sigma_total = k * sigma_previous; + *item = (sigma_total.powf(2.0) - sigma_previous.powf(2.0)).sqrt(); + + } + + kernels +} diff --git a/cv-sift/src/pyramid/mod.rs b/cv-sift/src/pyramid/mod.rs new file mode 100644 index 00000000..b09ed0c8 --- /dev/null +++ b/cv-sift/src/pyramid/mod.rs @@ -0,0 +1,8 @@ +mod base_image; +pub use base_image::*; + +mod kernels; +pub use kernels::*; + +mod octaves; +pub use octaves::*; \ No newline at end of file diff --git a/cv-sift/src/pyramid/octaves.rs b/cv-sift/src/pyramid/octaves.rs new file mode 100644 index 00000000..3a772713 --- /dev/null +++ b/cv-sift/src/pyramid/octaves.rs @@ -0,0 +1,15 @@ +/// Compute the number of octaves in the image pyramid as a function of height and width of the image. +/// # Examples +/// ``` +/// use cv_sift::pyramid::number_of_octaves; +/// +/// let num_octaves = number_of_octaves(223, 324); +/// assert_eq!(num_octaves, 7); +/// ``` +pub fn number_of_octaves(height: u32, width: u32) -> u32 { + if height < width { + ((height as f64).log2() - 1.0).round() as u32 + } else { + ((width as f64).log2() - 1.0).round() as u32 + } +} diff --git a/cv-sift/src/utils.rs b/cv-sift/src/utils.rs new file mode 100644 index 00000000..1c371428 --- /dev/null +++ b/cv-sift/src/utils.rs @@ -0,0 +1,53 @@ +/// Assert two given slices of floats are equal. +/// # Examples: +/// +/// ### Two exactly equal slices must be equal. +/// ``` +/// use cv_sift::utils::assert_similar; +/// +/// let a = vec![1.0, 2.0, 3.0]; +/// let b = vec![1.0, 2.0, 3.0]; +/// assert_similar(&a, &b); +/// ``` +/// +/// ### Elements differ no more than 1e-8. +/// ``` +/// use cv_sift::utils::assert_similar; +/// +/// let a = vec![1.0, 2.0 - 1e-8, 3.0]; +/// let b = vec![1.0, 2.0, 3.0 - 1e-9]; +/// +/// assert_similar(&a, &b); +/// ``` +pub fn assert_similar(v1: &[f64], v2: &[f64]) { + assert_eq!(v1.len(), v2.len()); + let result = v1 + .iter() + .zip(v2.iter()) + .all(|(&f1, &f2)| (f1 - f2).abs() <= 1e-8); + assert!(result); +} + +/// Assert two given slices of floats are not equal. +/// ### Some elements differ more than 1e-8. +/// ``` +/// use cv_sift::utils::assert_not_similar; +/// +/// let a = vec![1.0, 2.0 - 1e-8, 3.0]; +/// let b = vec![1.0 - 1e-5, 2.0, 3.0 - 1e-9]; +/// +/// assert_not_similar(&a, &b); +/// ``` +pub fn assert_not_similar(v1: &[f64], v2: &[f64]) { + assert_eq!(v1.len(), v2.len()); + let result = v1 + .iter() + .zip(v2.iter()) + .any(|(&f1, &f2)| (f1 - f2).abs() > 1e-8); + assert!(result); +} + +pub fn open>(p: P) -> crate::errors::Result { + image::open(p) + .map_err(|err| crate::errors::SIFTError::Unsupported(err.to_string())) +} diff --git a/cv-sift/tests/fixtures/box.png b/cv-sift/tests/fixtures/box.png new file mode 100644 index 00000000..6f01082f Binary files /dev/null and b/cv-sift/tests/fixtures/box.png differ diff --git a/cv-sift/tests/fixtures/sus.png b/cv-sift/tests/fixtures/sus.png new file mode 100644 index 00000000..ccf6d98e Binary files /dev/null and b/cv-sift/tests/fixtures/sus.png differ diff --git a/cv-sift/tests/test_base_image.rs b/cv-sift/tests/test_base_image.rs new file mode 100644 index 00000000..a162139f --- /dev/null +++ b/cv-sift/tests/test_base_image.rs @@ -0,0 +1,19 @@ +use test_case::test_case; + + +#[test_case("testdata/box.png"; "box_")] +fn test_base_image(in_path: &str) -> Result<(), Box> { + tracing_subscriber::fmt::init(); + + let img = cv_sift::utils::open(in_path.clone()).unwrap(); + eprintln!("height: {}, width: {}, name: {}", img.height(), img.width(), in_path); + + let base_img = cv_sift::pyramid::generate_base_image(&img, 1.6, 0.5)?; + + assert_eq!(img.height() * 2, base_img.height()); + assert_eq!(img.width() * 2, base_img.width()); + + eprintln!("height: {}, width: {}, name: {}", base_img.height(), base_img.width(), "base_img"); + + Ok(()) +} diff --git a/cv-sift/tests/test_gaussian_kernels.rs b/cv-sift/tests/test_gaussian_kernels.rs new file mode 100644 index 00000000..a86bb1f1 --- /dev/null +++ b/cv-sift/tests/test_gaussian_kernels.rs @@ -0,0 +1,11 @@ +use cv_sift::pyramid::gaussian_kernels; +use cv_sift::utils::assert_similar; +use test_case::test_case; + +#[test_case(0.0, 1, vec![0.; 4]; "0 point 0 and 1")] +#[test_case(2., 5, vec![2.0, 1.13050062, 1.2986042, 1.49170451, 1.71351851, 1.9683159, 2.26100123, 2.5972084]; "2 point 0 and 5")] +#[test_case(1.6, 3, vec![1.6, 1.2262735, 1.54500779, 1.94658784, 2.452547, 3.09001559]; "default params")] +fn test_gaussian_kernels(sigma: f64, num_intervals: usize, expected: Vec) { + let observed = gaussian_kernels(sigma, num_intervals); + assert_similar(&observed, &expected); +} \ No newline at end of file diff --git a/cv-sift/tests/test_number_of_octaves.rs b/cv-sift/tests/test_number_of_octaves.rs new file mode 100644 index 00000000..fcb0f449 --- /dev/null +++ b/cv-sift/tests/test_number_of_octaves.rs @@ -0,0 +1,8 @@ +use cv_sift::pyramid::number_of_octaves; +use test_case::test_case; + +#[test_case(223, 324, 7; "223x324")] +#[test_case(100, 200, 6; "100x200")] +fn test_number_of_octaves(height: u32, width: u32, expected: u32) { + assert_eq!(number_of_octaves(height, width), expected); +} diff --git a/cv-sift/tests/test_subtraction.rs b/cv-sift/tests/test_subtraction.rs new file mode 100644 index 00000000..8771ed56 --- /dev/null +++ b/cv-sift/tests/test_subtraction.rs @@ -0,0 +1,27 @@ +use cv_sift::conversion::try_get_rgb_32f; +use cv_sift::ext::ImageExt; +use cv_sift::ImageRgb32F; +use cv_sift::pyramid::subtract; + + +#[test] +fn subtraction_reflexivity() { + let ima = image::open("tests/fixtures/sus.png").unwrap(); + let minuend = try_get_rgb_32f(&ima).unwrap(); + let subtrahend = try_get_rgb_32f(&ima).unwrap(); + + let result = subtract(&minuend, &subtrahend).unwrap(); + assert!(result.is_zero()); +} + +#[test] +fn subtraction_of_zero_matrix() { + let ima = image::open("tests/fixtures/sus.png").unwrap(); + let minuend = try_get_rgb_32f(&ima).unwrap(); + + let subtrahend = ImageRgb32F::new(minuend.width(), minuend.height()); + assert!(subtrahend.is_zero()); + + let result = subtract(&minuend, &subtrahend).unwrap(); + assert!(result.is_same_as(&minuend)); +}