From 08738c5247da665805917ac93f2523322d5eab3e Mon Sep 17 00:00:00 2001 From: Bradley Gannon Date: Sun, 10 Mar 2024 09:56:06 -0400 Subject: [PATCH] Implement BRIEF descriptor (#490) --- examples/brief.rs | 111 +++++++++ src/binary_descriptors/brief.rs | 370 ++++++++++++++++++++++++++++ src/binary_descriptors/constants.rs | 2 + src/binary_descriptors/mod.rs | 167 +++++++++++++ src/corners.rs | 211 +++++++++++++++- src/lib.rs | 1 + 6 files changed, 861 insertions(+), 1 deletion(-) create mode 100644 examples/brief.rs create mode 100644 src/binary_descriptors/brief.rs create mode 100644 src/binary_descriptors/constants.rs create mode 100644 src/binary_descriptors/mod.rs diff --git a/examples/brief.rs b/examples/brief.rs new file mode 100644 index 00000000..4a6b8584 --- /dev/null +++ b/examples/brief.rs @@ -0,0 +1,111 @@ +//! An example of a common BRIEF descriptor workflow. First, read the images in. +//! Then, find keypoints and compute descriptors for small patches around them. +//! Match the descriptors and write the output. +//! +//! Use two images of a scene from slightly different positions. Consider +//! downloading images from any of the datasets used in the original BRIEF paper +//! [here][0]. The Wall/Viewpoint dataset seems to work well. +//! +//! [0]: https://www.robots.ox.ac.uk/~vgg/research/affine/ +//! +//! From the root of this crate, run: +//! +//! `cargo run --release --example brief -- ` + +use image::{open, GenericImage, ImageResult, Rgb}; +use imageproc::{ + binary_descriptors::{brief::brief, match_binary_descriptors, BinaryDescriptor}, + corners::{corners_fast9, Corner}, + definitions::Image, + drawing::draw_line_segment_mut, + point::Point, +}; + +use std::{env, path::Path}; + +fn filter_edge_keypoints(keypoint: &Corner, height: u32, width: u32, radius: u32) -> bool { + keypoint.x >= radius + && keypoint.x <= width - radius + && keypoint.y >= radius + && keypoint.y <= height - radius +} + +fn main() -> ImageResult<()> { + if env::args().len() != 4 { + panic!("Please enter two input files and one output file") + } + + let first_image_path_arg = env::args().nth(1).unwrap(); + let second_image_path_arg = env::args().nth(2).unwrap(); + let output_image_path_arg = env::args().nth(3).unwrap(); + + let first_image_path = Path::new(&first_image_path_arg); + let second_image_path = Path::new(&second_image_path_arg); + + if !first_image_path.is_file() { + panic!("First image file does not exist"); + } + if !second_image_path.is_file() { + panic!("Second image file does not exist"); + } + + let first_image = open(first_image_path)?.to_luma8(); + let second_image = open(second_image_path)?.to_luma8(); + + let (first_image_height, first_image_width) = (first_image.height(), first_image.width()); + let (second_image_height, second_image_width) = (second_image.height(), second_image.width()); + + const CORNER_THRESHOLD: u8 = 70; + let first_corners = corners_fast9(&first_image, CORNER_THRESHOLD) + .into_iter() + .filter(|c| filter_edge_keypoints(c, first_image_height, first_image_width, 16)) + .map(|c| c.into()) + .collect::>>(); + println!("Found {} corners in the first image", first_corners.len()); + let second_corners = corners_fast9(&second_image, CORNER_THRESHOLD) + .into_iter() + .filter(|c| filter_edge_keypoints(c, second_image_height, second_image_width, 16)) + .map(|c| c.into()) + .collect::>>(); + println!("Found {} corners in the second image", second_corners.len()); + + let (first_descriptors, test_pairs) = brief(&first_image, &first_corners, 256, None).unwrap(); + let (second_descriptors, _test_pairs) = + brief(&second_image, &second_corners, 256, Some(&test_pairs)).unwrap(); + println!("Computed descriptors"); + + let matches = match_binary_descriptors(&first_descriptors, &second_descriptors, 24, Some(0xc0)); + println!("Matched {} descriptor pairs", matches.len()); + + // now that we've matched descriptors in both images, put them side by side + // and draw lines connecting the descriptors together + let first_image = open(first_image_path)?.to_rgb8(); + let second_image = open(second_image_path)?.to_rgb8(); + + let (output_width, output_height) = ( + first_image_width + second_image_width, + u32::max(first_image_height, second_image_height), + ); + let mut output_image = Image::new(output_width, output_height); + output_image.copy_from(&first_image, 0, 0).unwrap(); + output_image + .copy_from(&second_image, first_image.width(), 0) + .unwrap(); + for keypoint_match in matches.iter() { + let start_point = keypoint_match.0.position(); + let end_point = keypoint_match.1.position(); + draw_line_segment_mut( + &mut output_image, + (start_point.x as f32, start_point.y as f32), + ( + (end_point.x + first_image.width()) as f32, + end_point.y as f32, + ), + Rgb([0, 255, 0]), + ) + } + output_image.save(&output_image_path_arg).unwrap(); + println!("Wrote output image to {}", output_image_path_arg); + + Ok(()) +} diff --git a/src/binary_descriptors/brief.rs b/src/binary_descriptors/brief.rs new file mode 100644 index 00000000..90f20b96 --- /dev/null +++ b/src/binary_descriptors/brief.rs @@ -0,0 +1,370 @@ +//! Structs and functions for finding and computing BRIEF descriptors as +//! described in [Calonder, et. al. (2010)][calonder]. +/// +/// [calonder]: https://www.cs.ubc.ca/~lowe/525/papers/calonder_eccv10.pdf +use image::{GenericImageView, GrayImage, ImageBuffer, Luma}; +use rand_distr::{Distribution, Normal}; + +use crate::{corners::Corner, integral_image::integral_image, point::Point}; + +use super::{ + constants::{BRIEF_PATCH_DIAMETER, BRIEF_PATCH_RADIUS}, + BinaryDescriptor, +}; + +/// BRIEF descriptor as described in [Calonder, et. al. (2010)][calonder]. +/// +/// [calonder]: https://www.cs.ubc.ca/~lowe/525/papers/calonder_eccv10.pdf +#[derive(Clone, PartialEq)] +pub struct BriefDescriptor { + /// Results of the pairwise pixel intensity tests that comprise this BRIEF + /// descriptor. + pub bits: Vec, + /// Pixel location and corner score of the keypoint associated with this + /// BRIEF descriptor. + pub corner: Corner, +} + +impl BinaryDescriptor for BriefDescriptor { + fn get_size(&self) -> u32 { + (self.bits.len() * 128) as u32 + } + + fn hamming_distance(&self, other: &Self) -> u32 { + assert_eq!(self.get_size(), other.get_size()); + self.bits + .iter() + .zip(other.bits.iter()) + .fold(0, |acc, x| acc + (x.0 ^ x.1).count_ones()) + } + + fn get_bit_subset(&self, bits: &[u32]) -> u128 { + assert!( + bits.len() <= 128, + "Can't extract more than 128 bits (found {})", + bits.len() + ); + let mut subset = 0; + for b in bits { + subset <<= 1; + // isolate the bit at index b + subset += (self.bits[(b / 128) as usize] >> (b % 128)) % 2; + } + subset + } + + fn position(&self) -> Point { + self.corner.into() + } +} + +/// Collection of two points that a BRIEF descriptor uses to generate its bits. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct TestPair { + /// The first point in the pair. + pub p0: Point, + /// The second point in the pair. + pub p1: Point, +} + +fn local_pixel_average( + integral_image: &ImageBuffer, Vec>, + x: u32, + y: u32, + radius: u32, +) -> u8 { + if radius == 0 { + return 0; + } + let y_min = if y < radius { 0 } else { y - radius }; + let x_min = if x < radius { 0 } else { x - radius }; + let y_max = u32::min(y + radius + 1, integral_image.height() - 1); + let x_max = u32::min(x + radius + 1, integral_image.width() - 1); + + let pixel_area = (y_max - y_min) * (x_max - x_min); + if pixel_area == 0 { + return 0; + } + + // UNSAFETY JUSTIFICATION + // + // Benefit + // + // Removing the unsafe pixel accesses in this function increases the + // runtimes for bench_brief_fixed_test_pairs_1000_keypoints, + // bench_brief_random_test_pairs_1000_keypoints, and + // bench_rotated_brief_1000_keypoints by about 40%, 30%, and 30%, + // respectively. + // + // Correctness + // + // The values of x_min, x_max, y_min, and y_max are all bounded between zero + // and the integral image dimensions by the checks at the top of this + // function. + let (bottom_right, top_left, top_right, bottom_left) = unsafe { + ( + integral_image.unsafe_get_pixel(x_max, y_max)[0], + integral_image.unsafe_get_pixel(x_min, y_min)[0], + integral_image.unsafe_get_pixel(x_max, y_min)[0], + integral_image.unsafe_get_pixel(x_min, y_max)[0], + ) + }; + let total_intensity = bottom_right + top_left - top_right - bottom_left; + (total_intensity / pixel_area) as u8 +} + +pub(crate) fn brief_impl( + integral_image: &ImageBuffer, Vec>, + keypoints: &[Point], + test_pairs: &[TestPair], + length: usize, +) -> Result, String> { + if length % 128 != 0 { + return Err(format!( + "BRIEF descriptor length must be a multiple of 128 bits (found {})", + length + )); + } + + if length != test_pairs.len() { + return Err(format!( + "BRIEF descriptor length must be equal to the number of test pairs ({} != {})", + length, + test_pairs.len() + )); + } + + let mut descriptors: Vec = Vec::with_capacity(keypoints.len()); + let (width, height) = (integral_image.width(), integral_image.height()); + + for keypoint in keypoints { + // if the keypoint is too close to the edge, return an error + if keypoint.x <= BRIEF_PATCH_RADIUS + || keypoint.x + BRIEF_PATCH_RADIUS >= width + || keypoint.y <= BRIEF_PATCH_RADIUS + || keypoint.y + BRIEF_PATCH_RADIUS >= height + { + return Err(format!( + "Found keypoint within {} px of image edge: ({}, {})", + BRIEF_PATCH_RADIUS + 1, + keypoint.x, + keypoint.y + )); + } + + let patch_top_left = Point { + x: keypoint.x - BRIEF_PATCH_RADIUS, + y: keypoint.y - BRIEF_PATCH_RADIUS, + }; + + let mut descriptor = BriefDescriptor { + bits: Vec::with_capacity(length / 128), + corner: Corner { + x: keypoint.x, + y: keypoint.y, + score: 0., + }, + }; + let mut descriptor_chunk = 0u128; + // for each test pair, compare the pixels within the patch at those points + for (idx, test_pair) in test_pairs.iter().enumerate() { + // if we've entered a new chunk, then save the previous one + if idx != 0 && idx % 128 == 0 { + descriptor.bits.push(descriptor_chunk); + descriptor_chunk = 0; + } + + let p0 = Point { + x: test_pair.p0.x + patch_top_left.x + 1, + y: test_pair.p0.y + patch_top_left.y + 1, + }; + let p1 = Point { + x: test_pair.p1.x + patch_top_left.x + 1, + y: test_pair.p1.y + patch_top_left.y + 1, + }; + + // if p0 < p1, then record true for this test; otherwise, record false + descriptor_chunk += (local_pixel_average(integral_image, p0.x, p0.y, 2) + < local_pixel_average(integral_image, p1.x, p1.y, 2)) + as u128; + descriptor_chunk <<= 1; + } + // save the final chunk too + descriptor.bits.push(descriptor_chunk); + descriptors.push(descriptor); + } + + Ok(descriptors) +} + +/// Generates BRIEF descriptors for small patches around keypoints in an image. +/// +/// Returns a tuple containing a vector of `BriefDescriptor` and a vector of +/// `TestPair`. All returned descriptors are based on the same `TestPair` set. +/// Patches are 31x31 pixels, so keypoints must be at least 17 pixels from any +/// edge. If any keypoints are too close to an edge, returns an error. +/// +/// `length` must be a multiple of 128 bits. Returns an error otherwise. +/// +/// If `override_test_pairs` is `Some`, then those test pairs are used, and none +/// are generated. Use this when you already have test pairs from another run +/// and want to compare the descriptors later. +/// +/// If `override_test_pairs` is `None`, then `TestPair`s are generated according +/// to an isotropic Gaussian. +/// +/// Calonder used Gaussian smoothing to decrease the effects of noise in the +/// patches. This is slow, even with a box filter approximation. For maximum +/// performance, the average intensities of sub-patches of radius 5 around the +/// test points are computed and used instead of the intensities of the test +/// points themselves. This is much faster because the averages come from +/// integral images. Calonder suggests that this approach may be faster, and +/// [Rublee et. al. (2012)][rublee] use this approach to quickly compute ORB +/// descriptors. +/// +/// [rublee]: http://www.gwylab.com/download/ORB_2012.pdf +/// +/// See [Calonder et. al. +/// (2010)][https://www.cs.ubc.ca/~lowe/525/papers/calonder_eccv10.pdf] +pub fn brief( + image: &GrayImage, + keypoints: &[Point], + length: usize, + override_test_pairs: Option<&Vec>, +) -> Result<(Vec, Vec), String> { + // if we have test pairs already, use them; otherwise, generate some + let test_pairs = if let Some(t) = override_test_pairs { + t.clone() + } else { + // generate a set of test pairs within a 31x31 grid with a Gaussian bias (sigma = 6.6) + let test_pair_distribution = Normal::new(BRIEF_PATCH_RADIUS as f32 + 1.0, 6.6).unwrap(); + let mut rng = rand::thread_rng(); + let mut test_pairs: Vec = Vec::with_capacity(length); + while test_pairs.len() < length { + let (x0, y0, x1, y1) = ( + test_pair_distribution.sample(&mut rng) as u32, + test_pair_distribution.sample(&mut rng) as u32, + test_pair_distribution.sample(&mut rng) as u32, + test_pair_distribution.sample(&mut rng) as u32, + ); + if x0 < BRIEF_PATCH_DIAMETER + && y0 < BRIEF_PATCH_DIAMETER + && x1 < BRIEF_PATCH_DIAMETER + && y1 < BRIEF_PATCH_DIAMETER + { + test_pairs.push(TestPair { + p0: Point::new(x0, y0), + p1: Point::new(x1, y1), + }); + } + } + test_pairs.clone() + }; + + let integral_image = integral_image(image); + + let descriptors = brief_impl(&integral_image, keypoints, &test_pairs, length)?; + + // return the descriptors for all the keypoints and the test pairs used + Ok((descriptors, test_pairs)) +} + +#[cfg(test)] +mod tests { + use rand::Rng; + use test::{black_box, Bencher}; + + use crate::utils::gray_bench_image; + + use super::*; + + #[test] + fn test_compute_hamming_distance() { + let d1 = BriefDescriptor { + bits: vec![ + 0xe41749023c71b74df05b57165519a180, + 0x9c06c422620e05d01105618cb3a2dcf1, + ], + corner: Corner::new(0, 0, 0.), + }; + let d2 = BriefDescriptor { + bits: vec![ + 0x8af22bb0596edc267c6f72cf425ebe1a, + 0xc1f6291520e474e8fa114e15420413d1, + ], + corner: Corner::new(0, 0, 0.), + }; + assert_eq!(d1.hamming_distance(&d2), 134); + } + + #[test] + fn test_get_bit_subset() { + let d = BriefDescriptor { + bits: vec![ + 0xdbe3de5bd950adf3d730034f9e4a55f7, + 0xf275f00f6243892a18ffefd0499996ee, + ], + corner: Corner::new(0, 0, 0.), + }; + let bits = vec![ + 226, 38, 212, 210, 60, 205, 68, 184, 47, 105, 152, 169, 11, 39, 76, 217, 183, 113, 189, + 251, 37, 181, 62, 28, 148, 92, 251, 77, 222, 148, 56, 142, + ]; + assert_eq!(d.get_bit_subset(&bits), 0b11001010011100011100011111011110); + } + + #[test] + fn test_local_pixel_average() { + let image = gray_image!( + 186, 106, 86, 22, 191, 10, 204, 217; + 37, 188, 82, 28, 99, 110, 166, 202; + 36, 97, 176, 54, 141, 42, 44, 40; + 248, 163, 218, 204, 117, 121, 151, 135; + 138, 100, 77, 115, 93, 246, 204, 163; + 123, 1, 104, 97, 67, 208, 0, 116; + 5, 237, 254, 171, 172, 165, 50, 39; + 92, 31, 238, 88, 44, 67, 140, 255 + ); + let integral_image: ImageBuffer, Vec> = integral_image(&image); + assert_eq!(local_pixel_average(&integral_image, 3, 3, 2), 117); + } + + #[bench] + #[ignore] + fn bench_brief_random_test_pairs_1000_keypoints(b: &mut Bencher) { + let image = gray_bench_image(640, 480); + let mut rng = rand::thread_rng(); + let keypoints = (0..1000) + .into_iter() + .map(|_| { + Point::new( + rng.gen_range(24, image.width() - 24), + rng.gen_range(24, image.height() - 24), + ) + }) + .collect::>>(); + b.iter(|| { + black_box(brief(&image, &keypoints, 256, None)).unwrap(); + }) + } + + #[bench] + #[ignore] + fn bench_brief_fixed_test_pairs_1000_keypoints(b: &mut Bencher) { + let image = gray_bench_image(640, 480); + let mut rng = rand::thread_rng(); + let keypoints = (0..1000) + .into_iter() + .map(|_| { + Point::new( + rng.gen_range(24, image.width() - 24), + rng.gen_range(24, image.height() - 24), + ) + }) + .collect::>>(); + let (_, test_pairs) = brief(&image, &keypoints, 256, None).unwrap(); + b.iter(|| { + black_box(brief(&image, &keypoints, 256, Some(&test_pairs))).unwrap(); + }) + } +} diff --git a/src/binary_descriptors/constants.rs b/src/binary_descriptors/constants.rs new file mode 100644 index 00000000..5fe5e90e --- /dev/null +++ b/src/binary_descriptors/constants.rs @@ -0,0 +1,2 @@ +pub const BRIEF_PATCH_RADIUS: u32 = 15; +pub const BRIEF_PATCH_DIAMETER: u32 = BRIEF_PATCH_RADIUS * 2 + 1; diff --git a/src/binary_descriptors/mod.rs b/src/binary_descriptors/mod.rs new file mode 100644 index 00000000..158fc113 --- /dev/null +++ b/src/binary_descriptors/mod.rs @@ -0,0 +1,167 @@ +//! Functions for generating and comparing compact binary patch descriptors. + +use rand::{rngs::StdRng, Rng, SeedableRng}; +use std::collections::HashMap; + +use crate::point::Point; + +pub mod brief; +mod constants; + +/// A feature descriptor whose value is given by a string of bits. +pub trait BinaryDescriptor { + /// Returns the length of the descriptor in bits. Typical values are 128, + /// 256, and 512. Will always be a multiple of 128. + fn get_size(&self) -> u32; + /// Returns the number of bits that are different between the two descriptors. + /// + /// Panics if the two descriptors have unequal lengths. The descriptors + /// should have been computed using the same set of test pairs, otherwise + /// comparing them has no meaning. + fn hamming_distance(&self, other: &Self) -> u32; + /// Given a set of bit indices, returns those bits from the descriptor as a + /// single concatenated value. + /// + /// Panics if `bits.len() > 128`. + fn get_bit_subset(&self, bits: &[u32]) -> u128; + /// Returns the pixel location of this binary descriptor in its associated + /// image. + fn position(&self) -> Point; +} + +/// For each descriptor in `d1`, find the descriptor in `d2` with the minimum +/// Hamming distance below `threshold`. If no such descriptor exists in `d2`, +/// the descriptor in `d1` is left unmatched. +/// +/// Descriptors in `d2` may be matched with more than one descriptor in `d1`. +/// +/// Uses [locality-sensitive hashing][lsh] (LSH) for efficient matching. The +/// number of tables is fixed at three, but the hash length is proportional to +/// the log of the size of the largest input array. +/// +/// Returns a vector of references describing the matched pairs. The first +/// reference is to a descriptor in `d1`, and the second reference is to an +/// index into `d2`. +/// +/// [lsh]: +/// https://en.wikipedia.org/wiki/Locality_sensitive_hashing#Bit_sampling_for_Hamming_distance +pub fn match_binary_descriptors<'a, T: BinaryDescriptor>( + d1: &'a [T], + d2: &'a [T], + threshold: u32, + seed: Option, +) -> Vec<(&'a T, &'a T)> { + // early return if either input is empty + if d1.is_empty() || d2.is_empty() { + return Vec::new(); + } + + let mut rng = if let Some(s) = seed { + StdRng::seed_from_u64(s) + } else { + StdRng::from_entropy() + }; + + // locality-sensitive hashing (LSH) + // this algorithm is log(d2.len()) but linear in d1.len(), so swap the inputs if needed + let (queries, database, swapped) = if d1.len() > d2.len() { + (d2, d1, true) + } else { + (d1, d2, false) + }; + + // build l hash tables by selecting k random bits from each descriptor + let l = 3; + // k grows as the log of the database size + // this keeps bucket size roughly constant + let k = (database.len() as f32).log2() as i32; + let mut hash_tables = Vec::with_capacity(l); + for _ in 0..l { + // choose k random bits (not necessarily unique) + let bits = (0..k) + .map(|_| rng.gen_range(0, queries[0].get_size())) + .collect::>(); + + let mut new_hashmap = HashMap::>::with_capacity(database.len()); + + // compute the hash of each descriptor in the database and store its index + // there will be collisions --- we want that to happen + for d in database.iter() { + let hash = d.get_bit_subset(&bits); + if let Some(v) = new_hashmap.get_mut(&hash) { + v.push(d); + } else { + new_hashmap.insert(hash, vec![d]); + } + } + hash_tables.push((bits, new_hashmap)); + } + + // find the hash buckets corresponding to each query descriptor + // then check all bucket members to find the (probable) best match + let mut matches = Vec::with_capacity(queries.len()); + for query in queries.iter() { + // find all buckets for the query descriptor + let mut candidates = Vec::with_capacity(l); + for (bits, table) in hash_tables.iter() { + let query_hash = query.get_bit_subset(bits); + if let Some(m) = table.get(&query_hash) { + for new_candidate in m.clone() { + candidates.push(new_candidate); + } + } + } + // perform linear scan to find the best match + let mut best_score = u32::MAX; + let mut best_candidate = None; + for c in candidates { + let distance = query.hamming_distance(c); + if distance < best_score { + best_score = distance; + best_candidate = Some(c); + } + } + // ignore the match if it's beyond our threshold + if best_score < threshold { + if swapped { + matches.push((best_candidate.unwrap(), query)); + } else { + matches.push((query, best_candidate.unwrap())); + } + } + } + matches +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::{binary_descriptors::brief::brief, utils::gray_bench_image}; + use test::{black_box, Bencher}; + + #[bench] + #[ignore] + fn bench_matcher_1000_keypoints_each(b: &mut Bencher) { + let image = gray_bench_image(640, 480); + let mut rng = rand::thread_rng(); + let keypoints = (0..1000) + .into_iter() + .map(|_| { + Point::new( + rng.gen_range(20, image.width() - 20), + rng.gen_range(20, image.height() - 20), + ) + }) + .collect::>>(); + let (first_descriptors, test_pairs) = brief(&image, &keypoints, 256, None).unwrap(); + let (second_descriptors, _) = brief(&image, &keypoints, 256, Some(&test_pairs)).unwrap(); + b.iter(|| { + black_box(match_binary_descriptors( + &first_descriptors, + &second_descriptors, + 24, + Some(0xc0), + )); + }); + } +} diff --git a/src/corners.rs b/src/corners.rs index 11e436b3..a76d0d28 100644 --- a/src/corners.rs +++ b/src/corners.rs @@ -1,7 +1,12 @@ //! Functions for detecting corners, also known as interest points. -use crate::definitions::{Position, Score}; +use crate::{ + definitions::{Position, Score}, + point::Point, +}; use image::{GenericImageView, GrayImage}; +use rand::{rngs::StdRng, SeedableRng}; +use rand_distr::Distribution; /// A location and score for a detected corner. /// The scores need not be comparable between different @@ -35,6 +40,12 @@ impl Position for Corner { } } +impl From for Point { + fn from(value: Corner) -> Self { + Point::new(value.x, value.y) + } +} + impl Score for Corner { fn score(&self) -> f32 { self.score @@ -90,6 +101,145 @@ pub fn corners_fast9(image: &GrayImage, threshold: u8) -> Vec { corners } +/// A FAST corner with associated orientation as described in [Rublee, et. al. +/// (2012)][rublee]. +/// +/// [rublee]: http://www.gwylab.com/download/ORB_2012.pdf +#[derive(Clone, Copy, PartialEq)] +pub struct OrientedFastCorner { + /// Location and FAST corner score of this corner in its associated image. + pub corner: Corner, + /// Orientation of this FAST corner as determined by computing the intensity + /// centroid of the local patch around the corner. + pub orientation: f32, +} + +fn intensity_centroid(image: &GrayImage, x: u32, y: u32, radius: u32) -> f32 { + let mut y_centroid: i32 = 0; + let mut x_centroid: i32 = 0; + + let (width, height) = image.dimensions(); + let x_min = if x < radius { 0 } else { x - radius }; + let y_min = if y < radius { 0 } else { y - radius }; + let y_max = u32::min(y + radius + 1, height); + let x_max = u32::min(x + radius + 1, width); + + let (mut x_count, mut y_count) = (-(radius as i32), (radius as i32)); + + for y in y_min..y_max { + for x in x_min..x_max { + // UNSAFETY JUSTIFICATION + // + // Benefit + // + // Removing all unsafe pixel accesses in this function increases the + // average runtime for bench_intensity_centroid by about 90%. + // + // Correctness + // + // x will always be greater than or equal to x_min and strictly less + // than x_max due to the range in this for loop. x_min will never be + // less than zero, and x_max will never be greater than the image + // width, both due to the checks earlier in this function. The same + // logic applies to y, y_min, and y_max. + let pixel = unsafe { image.unsafe_get_pixel(x, y).0[0] }; + x_centroid += x_count * (pixel as i32); + x_count += 1; + } + x_count = -(radius as i32); + } + + for x in x_min..x_max { + for y in y_min..y_max { + // See UNSAFETY JUSTIFICATION above. + let pixel = unsafe { image.unsafe_get_pixel(x, y).0[0] }; + y_centroid += y_count * (pixel as i32); + y_count -= 1; + } + y_count = radius as i32; + } + + // Important note: we flip the sign here because there are two coordinate + // systems in play. One is pixel space with the origin in the top left, and + // the other is ordinary Cartesian space with the origin in the bottom left. + // To make the math in later rotation code match the usual convention, we + // hide the coordinate conversion here. + -(y_centroid as f32).atan2(x_centroid as f32) +} + +/// Finds oriented FAST-9 corners as presented in [Rublee et. al. (2012)][rublee]. +/// +/// [rublee]: http://www.gwylab.com/download/ORB_2012.pdf +pub fn oriented_fast( + image: &GrayImage, + threshold: Option, + target_num_corners: usize, + edge_radius: u32, + seed: Option, +) -> Vec { + let (width, height) = image.dimensions(); + let (min_x, max_x) = (edge_radius, width - edge_radius); + let (min_y, max_y) = (edge_radius, height - edge_radius); + let mut corners = vec![]; + + let local_threshold = if let Some(t) = threshold { + t + } else { + // Take a sample of random pixels, compute their FAST scores, and set the + // threshold for the full image accordingly. + const NUM_SAMPLE_POINTS: usize = 1000; + let mut rng = if let Some(s) = seed { + StdRng::seed_from_u64(s) + } else { + StdRng::from_entropy() + }; + let dist_x = rand::distributions::Uniform::new(min_x, max_x); + let dist_y = rand::distributions::Uniform::new(min_y, max_y); + let sample_size = NUM_SAMPLE_POINTS.min((width * height) as usize); + let sample_coords: Vec> = (0..sample_size) + .map(|_| Point::new(dist_x.sample(&mut rng), dist_y.sample(&mut rng))) + .collect(); + let mut fast_scores: Vec = sample_coords + .iter() + .map(|c| fast_corner_score(image, 0, c.x, c.y, Fast::Nine)) + .collect(); + fast_scores.sort(); + let target_corner_fraction = (target_num_corners as f32) / ((width * height) as f32); + let fraction_idx = (NUM_SAMPLE_POINTS as f32 * (1. - target_corner_fraction)) as usize; + + fast_scores[fraction_idx] + }; + + // Iterate over every pixel in the image and find potential corners. + for y in edge_radius..height - edge_radius { + for x in edge_radius..width - edge_radius { + if is_corner_fast9(image, local_threshold, x, y) { + let score = fast_corner_score(image, local_threshold, x, y, Fast::Nine); + corners.push(Corner::new(x, y, score as f32)); + } + } + } + + // Sort descending by Harris corner measure. + corners.sort_by(|a, b| b.score.partial_cmp(&a.score).unwrap()); + + // Keep the top corners and discard the rest. + let top_corners = if corners.len() < target_num_corners { + &corners + } else { + &corners[..target_num_corners] + }; + + // Compute intensity centroids and return oriented FAST corners. + top_corners + .iter() + .map(|c| OrientedFastCorner { + corner: *c, + orientation: intensity_centroid(image, c.x, c.y, 15), + }) + .collect() +} + /// The score of a corner detected using the FAST /// detector is the largest threshold for which this /// pixel is still a corner. We input the threshold at which @@ -487,6 +637,65 @@ mod tests { assert!(is_corner_fast9(&image, 8, 3, 3)); } + #[test] + fn test_intensity_centroid() { + let image = gray_image!( + 00, 00, 10, 10, 10, 00, 00; + 00, 10, 00, 00, 00, 10, 00; + 10, 00, 00, 00, 00, 00, 10; + 10, 00, 00, 00, 00, 00, 10; + 00, 00, 00, 00, 00, 00, 10; + 00, 00, 00, 00, 00, 10, 00; + 00, 00, 00, 10, 10, 00, 00); + + assert_eq!( + intensity_centroid(&image, 3, 3, 3), + -std::f32::consts::FRAC_PI_4 + ); + } + + #[bench] + fn bench_intensity_centroid(b: &mut Bencher) { + let image = gray_image!( + 00, 00, 10, 10, 10, 00, 00; + 00, 10, 00, 00, 00, 10, 00; + 10, 00, 00, 00, 00, 00, 10; + 10, 00, 00, 00, 00, 00, 10; + 00, 00, 00, 00, 00, 00, 10; + 00, 00, 00, 00, 00, 10, 00; + 00, 00, 00, 10, 10, 00, 00); + + b.iter(|| black_box(intensity_centroid(&image, 3, 3, 3))); + } + + #[bench] + fn bench_oriented_fast_corner(b: &mut Bencher) { + let image = gray_image!( + 00, 00, 10, 10, 10, 00, 00; + 00, 10, 00, 00, 00, 10, 00; + 10, 00, 00, 00, 00, 00, 10; + 10, 00, 00, 00, 00, 00, 10; + 00, 00, 00, 00, 00, 00, 10; + 00, 00, 00, 00, 00, 10, 00; + 00, 00, 00, 10, 10, 00, 00); + + b.iter(|| black_box(oriented_fast(&image, Some(0), 1, 0, Some(0xc0)))); + } + + #[bench] + fn bench_oriented_fast_non_corner(b: &mut Bencher) { + let image = gray_image!( + 00, 00, 00, 00, 00, 00, 00; + 00, 00, 00, 00, 00, 00, 00; + 00, 00, 00, 00, 00, 00, 00; + 00, 00, 00, 00, 00, 00, 00; + 00, 00, 00, 00, 00, 00, 00; + 00, 00, 00, 00, 00, 00, 00; + 00, 00, 00, 00, 00, 00, 00); + + b.iter(|| black_box(oriented_fast(&image, Some(255), 0, 0, Some(0xc0)))); + } + #[bench] fn bench_is_corner_fast9_9_contiguous_lighter_pixels(b: &mut Bencher) { let image = black_box(gray_image!( diff --git a/src/lib.rs b/src/lib.rs index b83b75e5..bbe14f8d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -33,6 +33,7 @@ extern crate assert_approx_eq; #[macro_use] pub mod utils; +pub mod binary_descriptors; pub mod contours; pub mod contrast; pub mod corners;