spatialrust_vision/
corners.rs

1//! Classical corner and accelerated segment-test detectors.
2
3use std::cmp::Ordering;
4
5use spatialrust_image::{Image, ImageView};
6
7use crate::border::fetch;
8use crate::{sobel, BorderMode, Keypoint2, VisionError, VisionResult};
9
10/// Shared response thresholding and spatial suppression options.
11#[derive(Clone, Copy, Debug, PartialEq)]
12pub struct CornerSelectionOptions {
13    /// Maximum retained corners; zero keeps every accepted corner.
14    pub max_corners: usize,
15    /// Fraction of the strongest positive response in `(0, 1]`.
16    pub quality_level: f32,
17    /// Minimum Euclidean distance between retained corners in pixels.
18    pub min_distance: f32,
19    /// Odd structure-tensor summation window size.
20    pub block_size: usize,
21    /// Sobel aperture size: 3, 5, or 7.
22    pub gradient_size: usize,
23    /// Source-image border extrapolation used by derivatives.
24    pub border: BorderMode<u8, 1>,
25}
26
27impl Default for CornerSelectionOptions {
28    fn default() -> Self {
29        Self {
30            max_corners: 0,
31            quality_level: 0.01,
32            min_distance: 1.0,
33            block_size: 3,
34            gradient_size: 3,
35            border: BorderMode::Reflect101,
36        }
37    }
38}
39
40impl CornerSelectionOptions {
41    fn validate(self) -> VisionResult<Self> {
42        if !self.quality_level.is_finite() || self.quality_level <= 0.0 || self.quality_level > 1.0
43        {
44            return Err(VisionError::InvalidParameter(
45                "corner quality_level must be finite and in (0, 1]".into(),
46            ));
47        }
48        if !self.min_distance.is_finite() || self.min_distance < 0.0 {
49            return Err(VisionError::InvalidParameter(
50                "corner min_distance must be finite and non-negative".into(),
51            ));
52        }
53        if self.block_size == 0 || self.block_size % 2 == 0 {
54            return Err(VisionError::InvalidParameter(
55                "corner block_size must be positive and odd".into(),
56            ));
57        }
58        if !matches!(self.gradient_size, 3 | 5 | 7) {
59            return Err(VisionError::InvalidParameter(
60                "corner gradient_size must be 3, 5, or 7".into(),
61            ));
62        }
63        Ok(self)
64    }
65}
66
67/// Harris detector configuration.
68#[derive(Clone, Copy, Debug, PartialEq)]
69pub struct HarrisOptions {
70    /// Shared candidate selection options.
71    pub selection: CornerSelectionOptions,
72    /// Harris trace penalty, normally near `0.04`.
73    pub k: f32,
74}
75
76impl Default for HarrisOptions {
77    fn default() -> Self {
78        Self { selection: CornerSelectionOptions::default(), k: 0.04 }
79    }
80}
81
82/// Shi–Tomasi minimum-eigenvalue detector configuration.
83#[derive(Clone, Copy, Debug, Default, PartialEq)]
84pub struct ShiTomasiOptions {
85    /// Shared candidate selection options.
86    pub selection: CornerSelectionOptions,
87}
88
89/// Computes an unnormalized Harris structure-tensor response image.
90pub fn harris_response(
91    input: ImageView<'_, u8, 1>,
92    block_size: usize,
93    gradient_size: usize,
94    k: f32,
95    border: BorderMode<u8, 1>,
96) -> VisionResult<Image<f32, 1>> {
97    if !k.is_finite() || k < 0.0 {
98        return Err(VisionError::InvalidParameter(
99            "Harris k must be finite and non-negative".into(),
100        ));
101    }
102    corner_response(input, block_size, gradient_size, border, |a, b, c| {
103        a.mul_add(c, -(b * b)) - k * (a + c) * (a + c)
104    })
105}
106
107/// Computes an unnormalized Shi–Tomasi minimum-eigenvalue response image.
108pub fn shi_tomasi_response(
109    input: ImageView<'_, u8, 1>,
110    block_size: usize,
111    gradient_size: usize,
112    border: BorderMode<u8, 1>,
113) -> VisionResult<Image<f32, 1>> {
114    corner_response(input, block_size, gradient_size, border, |a, b, c| {
115        0.5 * (a + c - ((a - c) * (a - c) + 4.0 * b * b).sqrt())
116    })
117}
118
119/// Detects Harris corners and returns strongest-first keypoints.
120pub fn detect_harris(
121    input: ImageView<'_, u8, 1>,
122    options: HarrisOptions,
123) -> VisionResult<Vec<Keypoint2>> {
124    let selection = options.selection.validate()?;
125    let response = harris_response(
126        input,
127        selection.block_size,
128        selection.gradient_size,
129        options.k,
130        selection.border,
131    )?;
132    select_corner_responses(response.view(), selection)
133}
134
135/// Detects Shi–Tomasi corners and returns strongest-first keypoints.
136pub fn detect_shi_tomasi(
137    input: ImageView<'_, u8, 1>,
138    options: ShiTomasiOptions,
139) -> VisionResult<Vec<Keypoint2>> {
140    let selection = options.selection.validate()?;
141    let response = shi_tomasi_response(
142        input,
143        selection.block_size,
144        selection.gradient_size,
145        selection.border,
146    )?;
147    select_corner_responses(response.view(), selection)
148}
149
150fn corner_response(
151    input: ImageView<'_, u8, 1>,
152    block_size: usize,
153    gradient_size: usize,
154    border: BorderMode<u8, 1>,
155    score: impl Fn(f32, f32, f32) -> f32,
156) -> VisionResult<Image<f32, 1>> {
157    let validation = CornerSelectionOptions {
158        block_size,
159        gradient_size,
160        border,
161        ..CornerSelectionOptions::default()
162    }
163    .validate()?;
164    let gradient_x = sobel(input, 1, 0, validation.gradient_size, 1.0, 0.0, border)?;
165    let gradient_y = sobel(input, 0, 1, validation.gradient_size, 1.0, 0.0, border)?;
166    let products_xx = Image::try_new(
167        input.width(),
168        input.height(),
169        gradient_x.as_slice().iter().map(|value| value * value).collect(),
170    )?;
171    let products_xy = Image::try_new(
172        input.width(),
173        input.height(),
174        gradient_x.as_slice().iter().zip(gradient_y.as_slice()).map(|(x, y)| x * y).collect(),
175    )?;
176    let products_yy = Image::try_new(
177        input.width(),
178        input.height(),
179        gradient_y.as_slice().iter().map(|value| value * value).collect(),
180    )?;
181    let product_border = product_border(border);
182    let radius = (block_size / 2) as isize;
183    let mut output = Vec::with_capacity(input.width() * input.height());
184    for y in 0..input.height() {
185        for x in 0..input.width() {
186            let mut a = 0.0_f32;
187            let mut b = 0.0_f32;
188            let mut c = 0.0_f32;
189            for dy in -radius..=radius {
190                for dx in -radius..=radius {
191                    a +=
192                        fetch(products_xx.view(), x as isize + dx, y as isize + dy, product_border)
193                            [0];
194                    b +=
195                        fetch(products_xy.view(), x as isize + dx, y as isize + dy, product_border)
196                            [0];
197                    c +=
198                        fetch(products_yy.view(), x as isize + dx, y as isize + dy, product_border)
199                            [0];
200                }
201            }
202            output.push(score(a, b, c));
203        }
204    }
205    Ok(Image::try_new(input.width(), input.height(), output)?)
206}
207
208fn product_border(border: BorderMode<u8, 1>) -> BorderMode<f32, 1> {
209    match border {
210        BorderMode::Constant(_) => BorderMode::Constant([0.0]),
211        BorderMode::Replicate => BorderMode::Replicate,
212        BorderMode::Reflect => BorderMode::Reflect,
213        BorderMode::Reflect101 => BorderMode::Reflect101,
214        BorderMode::Wrap => BorderMode::Wrap,
215    }
216}
217
218fn select_corner_responses(
219    response: ImageView<'_, f32, 1>,
220    options: CornerSelectionOptions,
221) -> VisionResult<Vec<Keypoint2>> {
222    let maximum = (0..response.height())
223        .flat_map(|y| response.row(y).expect("coordinate in bounds").iter().copied())
224        .filter(|value| value.is_finite())
225        .fold(f32::NEG_INFINITY, f32::max);
226    if maximum <= 0.0 || !maximum.is_finite() {
227        return Ok(Vec::new());
228    }
229    if response.width() < 3 || response.height() < 3 {
230        return Ok(Vec::new());
231    }
232    let threshold = maximum * options.quality_level;
233    let mut candidates = Vec::new();
234    // OpenCV goodFeaturesToTrack applies 3x3 dilation/NMS only to pixels with
235    // a complete immediate neighborhood, excluding the outermost image ring.
236    for y in 1..response.height() - 1 {
237        for x in 1..response.width() - 1 {
238            let value = response.get(x, y).expect("coordinate in bounds")[0];
239            if value <= threshold || !is_local_maximum(response, x, y, value) {
240                continue;
241            }
242            candidates.push((x, y, value));
243        }
244    }
245    candidates.sort_by(|left, right| {
246        right
247            .2
248            .partial_cmp(&left.2)
249            .unwrap_or(Ordering::Equal)
250            .then_with(|| left.1.cmp(&right.1))
251            .then_with(|| left.0.cmp(&right.0))
252    });
253    let minimum_squared = options.min_distance * options.min_distance;
254    let mut selected = Vec::<Keypoint2>::new();
255    for (x, y, value) in candidates {
256        if options.min_distance > 0.0
257            && selected.iter().any(|keypoint| {
258                let dx = keypoint.x() - x as f32;
259                let dy = keypoint.y() - y as f32;
260                dx.mul_add(dx, dy * dy) < minimum_squared
261            })
262        {
263            continue;
264        }
265        selected.push(
266            Keypoint2::try_new(x as f32, y as f32, value)?.with_size(options.block_size as f32)?,
267        );
268        if options.max_corners != 0 && selected.len() == options.max_corners {
269            break;
270        }
271    }
272    Ok(selected)
273}
274
275fn is_local_maximum(response: ImageView<'_, f32, 1>, x: usize, y: usize, value: f32) -> bool {
276    for dy in -1_isize..=1 {
277        for dx in -1_isize..=1 {
278            if dx == 0 && dy == 0 {
279                continue;
280            }
281            let nx = x as isize + dx;
282            let ny = y as isize + dy;
283            if nx >= 0
284                && ny >= 0
285                && (nx as usize) < response.width()
286                && (ny as usize) < response.height()
287                && response.get(nx as usize, ny as usize).expect("checked coordinate")[0] > value
288            {
289                return false;
290            }
291        }
292    }
293    true
294}
295
296/// FAST detector configuration for the standard 9-of-16 radius-three circle.
297#[derive(Clone, Copy, Debug, PartialEq, Eq)]
298pub struct FastOptions {
299    /// Strict center-to-circle intensity threshold.
300    pub threshold: u8,
301    /// Retain only strict 3x3 maxima of the FAST score.
302    pub nonmax_suppression: bool,
303}
304
305impl Default for FastOptions {
306    fn default() -> Self {
307        Self { threshold: 10, nonmax_suppression: true }
308    }
309}
310
311const FAST_CIRCLE: [(isize, isize); 16] = [
312    (0, -3),
313    (1, -3),
314    (2, -2),
315    (3, -1),
316    (3, 0),
317    (3, 1),
318    (2, 2),
319    (1, 3),
320    (0, 3),
321    (-1, 3),
322    (-2, 2),
323    (-3, 1),
324    (-3, 0),
325    (-3, -1),
326    (-2, -2),
327    (-1, -3),
328];
329
330/// Detects FAST-9/16 corners in scan order, matching OpenCV's radius-three ROI.
331pub fn detect_fast(
332    input: ImageView<'_, u8, 1>,
333    options: FastOptions,
334) -> VisionResult<Vec<Keypoint2>> {
335    if input.width() < 7 || input.height() < 7 {
336        return Ok(Vec::new());
337    }
338    let mut scores = vec![0_u8; input.width() * input.height()];
339    let mut candidates = Vec::new();
340    for y in 3..input.height() - 3 {
341        for x in 3..input.width() - 3 {
342            let score = fast_score(input, x, y);
343            if score >= options.threshold {
344                scores[y * input.width() + x] = score;
345                candidates.push((x, y));
346            }
347        }
348    }
349    let mut keypoints = Vec::new();
350    for (x, y) in candidates {
351        let score = scores[y * input.width() + x];
352        if options.nonmax_suppression {
353            let strict_maximum = (-1_isize..=1).all(|dy| {
354                (-1_isize..=1).all(|dx| {
355                    (dx == 0 && dy == 0)
356                        || score
357                            > scores[(y as isize + dy) as usize * input.width()
358                                + (x as isize + dx) as usize]
359                })
360            });
361            if !strict_maximum {
362                continue;
363            }
364        }
365        keypoints.push(
366            Keypoint2::try_new(
367                x as f32,
368                y as f32,
369                if options.nonmax_suppression { f32::from(score) } else { 0.0 },
370            )?
371            .with_size(7.0)?,
372        );
373    }
374    Ok(keypoints)
375}
376
377fn fast_score(input: ImageView<'_, u8, 1>, x: usize, y: usize) -> u8 {
378    let center = i16::from(input.get(x, y).expect("coordinate in bounds")[0]);
379    let differences = std::array::from_fn::<_, 16, _>(|index| {
380        let (dx, dy) = FAST_CIRCLE[index];
381        i16::from(
382            input
383                .get((x as isize + dx) as usize, (y as isize + dy) as usize)
384                .expect("FAST radius checked")[0],
385        ) - center
386    });
387    let mut best = 0_i16;
388    for start in 0..16 {
389        let mut bright = i16::MAX;
390        let mut dark = i16::MAX;
391        for offset in 0..9 {
392            let difference = differences[(start + offset) % 16];
393            bright = bright.min(difference);
394            dark = dark.min(-difference);
395        }
396        best = best.max(bright).max(dark);
397    }
398    best.saturating_sub(1).clamp(0, 255) as u8
399}
400
401#[cfg(test)]
402mod tests {
403    use super::{
404        detect_fast, detect_harris, detect_shi_tomasi, FastOptions, HarrisOptions, ShiTomasiOptions,
405    };
406    use spatialrust_image::{Image, ImageRegion};
407
408    fn square() -> Image<u8, 1> {
409        let mut image = Image::try_new(17, 15, vec![0; 17 * 15]).unwrap();
410        for y in 4..11 {
411            for x in 5..13 {
412                image.get_mut(x, y).unwrap()[0] = 255;
413            }
414        }
415        image
416    }
417
418    #[test]
419    fn harris_and_shi_tomasi_find_square_corners() {
420        let image = square();
421        let harris = detect_harris(image.view(), HarrisOptions::default()).unwrap();
422        let shi = detect_shi_tomasi(image.view(), ShiTomasiOptions::default()).unwrap();
423        for corners in [&harris, &shi] {
424            assert!(corners.iter().any(|point| point.x() <= 6.0 && point.y() <= 5.0));
425            assert!(corners.iter().any(|point| point.x() >= 11.0 && point.y() >= 9.0));
426        }
427    }
428
429    #[test]
430    fn corner_detectors_accept_strided_roi() {
431        let image = square();
432        let roi = image.view().subview(ImageRegion::new(2, 2, 13, 11)).unwrap();
433        assert!(!detect_harris(roi, HarrisOptions::default()).unwrap().is_empty());
434    }
435
436    #[test]
437    fn fast_detects_contrast_corner_and_handles_tiny_images() {
438        let mut image = Image::try_new(9, 9, vec![0; 81]).unwrap();
439        image.get_mut(4, 4).unwrap()[0] = 255;
440        let corners =
441            detect_fast(image.view(), FastOptions { threshold: 20, nonmax_suppression: true })
442                .unwrap();
443        assert_eq!(corners.len(), 1);
444        assert_eq!((corners[0].x(), corners[0].y()), (4.0, 4.0));
445        assert!(corners.iter().all(|point| point.size() == 7.0));
446        let tiny = Image::try_new(6, 6, vec![0; 36]).unwrap();
447        assert!(detect_fast(tiny.view(), FastOptions::default()).unwrap().is_empty());
448    }
449}
spatialrust_vision/corners.rs

spatialrust_vision/
corners.rs
