eventcv_core/

features.rs

//! Event feature detection (OpenCV `features2d` analogue). Both detectors maintain a per-polarity
//! Surface of Active Events (SAE — the latest timestamp seen at each pixel) and return a **new**
//! [`EventStream`] holding only the events that sit on a moving corner, so detection **chains**
//! like a denoising filter. Both assume events arrive in **ascending time order** (what the
//! readers produce); call [`EventStream::sort_by_time`] first if a stream might be unordered.

use crate::{EventStream, EventStreamBuilder};

/// Bresenham circle of radius 3 (16 pixels), in contiguous clockwise order — the FAST ring.
const INNER_CIRCLE: [(i32, i32); 16] = [
    (0, 3),
    (1, 3),
    (2, 2),
    (3, 1),
    (3, 0),
    (3, -1),
    (2, -2),
    (1, -3),
    (0, -3),
    (-1, -3),
    (-2, -2),
    (-3, -1),
    (-3, 0),
    (-3, 1),
    (-2, 2),
    (-1, 3),
];

/// Bresenham circle of radius 4 (20 pixels), contiguous clockwise — eFAST's outer ring.
const OUTER_CIRCLE: [(i32, i32); 20] = [
    (0, 4),
    (1, 4),
    (2, 3),
    (3, 2),
    (4, 1),
    (4, 0),
    (4, -1),
    (3, -2),
    (2, -3),
    (1, -4),
    (0, -4),
    (-1, -4),
    (-2, -3),
    (-3, -2),
    (-4, -1),
    (-4, 0),
    (-4, 1),
    (-3, 2),
    (-2, 3),
    (-1, 4),
];

/// eFAST arc-length bounds: a corner's most-recent pixels form a contiguous arc this many pixels
/// long. A straight edge fills ~half the ring (too long); noise fills too little.
const INNER_ARC: (usize, usize) = (3, 6);
const OUTER_ARC: (usize, usize) = (4, 8);

/// Harris uses a fixed 9×9 window (radius 4) of the time surface around each event.
const HARRIS_RADIUS: usize = 4;
/// Harris' empirical sensitivity constant `k` in `det - k·trace²`.
const HARRIS_K: f64 = 0.04;

impl EventStream {
    /// **eFAST** event corner detector (Mueggler et al., *Fast Event-based Corner Detection*,
    /// BMVC 2017). For each event it updates its polarity's SAE, then tests two Bresenham rings
    /// (radius 3 and 4) around the pixel: an event is a corner when, on **both** rings, the most
    /// recent timestamps form a contiguous arc within the `INNER_ARC`/`OUTER_ARC` bounds — the
    /// signature of a moving corner rather than a straight edge. Events too close to the border to
    /// evaluate the outer ring are dropped. Returns the corner events as a new stream.
    pub fn efast(&self) -> EventStream {
        let (width, height) = self.sensor_size();
        let (xs, ys, ts, ps) = (self.xs(), self.ys(), self.ts(), self.ps());
        let mut builder =
            EventStreamBuilder::with_capacity(width, height, self.timestamp_scale_ms(), self.len());
        if width == 0 || height == 0 {
            return builder.build();
        }
        // Separate surfaces per polarity, as in the paper.
        let mut sae_on = vec![i64::MIN; width * height];
        let mut sae_off = vec![i64::MIN; width * height];

        for index in 0..self.len() {
            let (x, y, t, p) = (xs[index] as usize, ys[index] as usize, ts[index], ps[index]);
            let sae = if p { &mut sae_on } else { &mut sae_off };
            sae[y * width + x] = t; // the current pixel is the newest by construction

            // Need radius 4 clear of every border to sample the outer ring.
            if x < 4 || y < 4 || x + 4 >= width || y + 4 >= height {
                continue;
            }
            if ring_is_corner(sae, x, y, width, &INNER_CIRCLE, INNER_ARC)
                && ring_is_corner(sae, x, y, width, &OUTER_CIRCLE, OUTER_ARC)
            {
                builder.push(xs[index], ys[index], t, p);
            }
        }
        builder.build()
    }

    /// **Harris corner score on the Surface of Active Events.** For each event it updates a
    /// merged SAE of raw latest timestamps, then computes the Harris response
    /// `det(M) - k·trace(M)²` of the structure tensor `M = Σ ∇T ∇Tᵀ` of the SAE's spatial
    /// gradient over a 9×9 window. Because the SAE is a local time *ramp*, a straight moving edge
    /// has a constant gradient direction (rank-1 `M`, `R < 0`) while a corner mixes gradient
    /// directions (rank-2 `M`, `R > 0`) — so the default `threshold = 0` keeps corners and rejects
    /// edges. Raise `threshold` to be stricter. Returns the corner events as a new stream; a
    /// score-based complement to [`Self::efast`].
    pub fn harris_corners(&self, threshold: f64) -> EventStream {
        let (width, height) = self.sensor_size();
        let (xs, ys, ts, ps) = (self.xs(), self.ys(), self.ts(), self.ps());
        let scale = self.timestamp_scale_ms();
        let mut builder =
            EventStreamBuilder::with_capacity(width, height, self.timestamp_scale_ms(), self.len());
        let margin = HARRIS_RADIUS + 1; // central differences need one pixel beyond the window
        if width < 2 * margin + 1 || height < 2 * margin + 1 {
            return builder.build();
        }
        // One merged ramp — flow/corner geometry is independent of contrast polarity.
        let mut sae = vec![i64::MIN; width * height];

        for index in 0..self.len() {
            let (x, y, t, p) = (xs[index] as usize, ys[index] as usize, ts[index], ps[index]);
            sae[y * width + x] = t;

            if x < margin || y < margin || x + margin >= width || y + margin >= height {
                continue;
            }
            if harris_response(&sae, x, y, width, scale) > threshold {
                builder.push(xs[index], ys[index], t, p);
            }
        }
        builder.build()
    }
}

/// Tests whether the ring's most-recent timestamps form a contiguous arc whose length lies in
/// `[min_arc, max_arc]` — i.e. some rotation of the ring has a run of pixels all strictly newer
/// than every pixel outside the run. `circle` offsets are contiguous around the ring so wrap-around
/// runs are handled by indexing modulo the ring length.
fn ring_is_corner(
    sae: &[i64],
    cx: usize,
    cy: usize,
    width: usize,
    circle: &[(i32, i32)],
    (min_arc, max_arc): (usize, usize),
) -> bool {
    let n = circle.len();
    let mut times = [i64::MIN; 20]; // 20 = the largest ring
    for (slot, &(dx, dy)) in times.iter_mut().zip(circle) {
        let px = (cx as i32 + dx) as usize;
        let py = (cy as i32 + dy) as usize;
        *slot = sae[py * width + px];
    }
    let times = &times[..n];

    for length in min_arc..=max_arc {
        for start in 0..n {
            let arc_min = (0..length)
                .map(|k| times[(start + k) % n])
                .min()
                .expect("arc length is at least 1");
            let rest_max = (length..n)
                .map(|k| times[(start + k) % n])
                .max()
                .expect("ring is longer than the arc");
            if arc_min > rest_max {
                return true;
            }
        }
    }
    false
}

/// Minimum number of valid gradient samples in the window for a Harris score to be meaningful.
const HARRIS_MIN_SAMPLES: usize = 3;

/// Harris response of the raw SAE ramp over the `(2R+1)²` window around `(cx, cy)`. Reads the SAE
/// as a time surface `T` (in ms), takes central-difference gradients wherever both neighbours have
/// fired, and returns `det(M) - k·trace(M)²` for the structure tensor `M = Σ ∇T ∇Tᵀ`. Unfired
/// pixels are skipped (their gradient is undefined), so a window with too little structure returns
/// `-∞` (never a corner). The caller guarantees the window lies `HARRIS_RADIUS + 1` inside every
/// border, so the `x ± 1` / `y ± 1` reads are in bounds.
fn harris_response(sae: &[i64], cx: usize, cy: usize, width: usize, scale: f64) -> f64 {
    let t_at = |x: usize, y: usize| -> Option<f64> {
        let t = sae[y * width + x];
        (t != i64::MIN).then_some(t as f64 * scale)
    };
    // Gradient at a fired pixel from central differences, needing both neighbours on each axis.
    let gradient = |x: usize, y: usize| -> Option<(f64, f64)> {
        t_at(x, y)?; // the pixel itself must have fired
        let gx = (t_at(x + 1, y)? - t_at(x - 1, y)?) / 2.0;
        let gy = (t_at(x, y + 1)? - t_at(x, y - 1)?) / 2.0;
        Some((gx, gy))
    };

    let (mut sxx, mut syy, mut sxy) = (0.0, 0.0, 0.0);
    let mut samples = 0;
    for y in cy - HARRIS_RADIUS..=cy + HARRIS_RADIUS {
        for x in cx - HARRIS_RADIUS..=cx + HARRIS_RADIUS {
            if let Some((ix, iy)) = gradient(x, y) {
                sxx += ix * ix;
                syy += iy * iy;
                sxy += ix * iy;
                samples += 1;
            }
        }
    }
    if samples < HARRIS_MIN_SAMPLES {
        return f64::NEG_INFINITY;
    }
    let det = sxx * syy - sxy * sxy;
    let trace = sxx + syy;
    det - HARRIS_K * trace * trace
}

#[cfg(test)]
mod tests {
    use ndarray::{array, Array2};

    use super::{ring_is_corner, INNER_ARC, INNER_CIRCLE};
    use crate::EventStream;

    fn empty(width: usize, height: usize) -> EventStream {
        EventStream::from_array2(Array2::zeros((0, 4)), width, height, 0.001)
    }

    // A tiny SAE laid out so a 4-pixel recent arc appears on the 16-pixel inner ring.
    fn sae_with_arc(recent: &[usize]) -> Vec<i64> {
        // 7×7 grid, centre at (3,3); ring pixels default old, `recent` indices set newest.
        let width = 7;
        let mut sae = vec![0_i64; width * width];
        for (k, &(dx, dy)) in INNER_CIRCLE.iter().enumerate() {
            let x = (3 + dx) as usize;
            let y = (3 + dy) as usize;
            sae[y * width + x] = if recent.contains(&k) { 100 } else { 1 };
        }
        sae
    }

    #[test]
    fn ring_detects_contiguous_arc_and_rejects_scattered_or_long() {
        let width = 7;
        // All ring pixels equal → flat, no corner.
        let flat = sae_with_arc(&[]);
        assert!(!ring_is_corner(
            &flat,
            3,
            3,
            width,
            &INNER_CIRCLE,
            INNER_ARC
        ));
        // 4 contiguous newest pixels → corner.
        let corner = sae_with_arc(&[2, 3, 4, 5]);
        assert!(ring_is_corner(
            &corner,
            3,
            3,
            width,
            &INNER_CIRCLE,
            INNER_ARC
        ));
        // Half the ring newest (8 pixels) → straight edge, exceeds max arc.
        let edge = sae_with_arc(&[0, 1, 2, 3, 4, 5, 6, 7]);
        assert!(!ring_is_corner(
            &edge,
            3,
            3,
            width,
            &INNER_CIRCLE,
            INNER_ARC
        ));
        // Two isolated newest pixels on opposite sides → not contiguous.
        let scattered = sae_with_arc(&[0, 8]);
        assert!(!ring_is_corner(
            &scattered,
            3,
            3,
            width,
            &INNER_CIRCLE,
            INNER_ARC
        ));
    }

    #[test]
    fn efast_empty_stream_is_empty() {
        assert_eq!(empty(20, 20).efast().len(), 0);
        assert_eq!(empty(0, 0).efast().len(), 0);
    }

    #[test]
    fn efast_and_harris_return_a_subset() {
        // A moving L-corner: a horizontal then vertical arm sweeping across time.
        let mut rows = Vec::new();
        let mut t = 0_u64;
        for x in 0..20u64 {
            rows.push([x, 10, t, 1]);
            t += 10;
        }
        for y in 0..20u64 {
            rows.push([10, y, t, 1]);
            t += 10;
        }
        let events = Array2::from_shape_vec((rows.len(), 4), rows.concat()).unwrap();
        let stream = EventStream::from_array2(events, 20, 20, 0.001);

        let corners = stream.efast();
        assert!(corners.len() <= stream.len());
        let harris = stream.harris_corners(0.0);
        assert!(harris.len() <= stream.len());
    }

    #[test]
    fn harris_empty_and_tiny_sensor_are_empty() {
        assert_eq!(empty(20, 20).harris_corners(0.0).len(), 0);
        // Sensor smaller than the Harris window → nothing evaluable.
        assert_eq!(
            EventStream::from_array2(array![[1, 1, 5, 1]], 4, 4, 0.001)
                .harris_corners(0.0)
                .len(),
            0
        );
    }

    /// A straight edge sweeping across the sensor is rank-1 on the SAE ramp (`R < 0`), so the
    /// `threshold = 0` Harris keeps almost nothing — the property that makes it a corner detector
    /// rather than an edge detector.
    #[test]
    fn harris_rejects_a_straight_moving_edge() {
        let mut rows = Vec::new();
        let mut t = 0_u64;
        for step in 0..30u64 {
            for y in 2..28u64 {
                rows.push([5 + step, y, t, 1]); // vertical edge column advancing in +x
            }
            t += 100;
        }
        let events = Array2::from_shape_vec((rows.len(), 4), rows.concat()).unwrap();
        let stream = EventStream::from_array2(events, 40, 30, 0.001);
        let kept = stream.harris_corners(0.0).len();
        // Far below 5% of the edge events survive (no true corner is present).
        assert!(
            kept * 20 < stream.len(),
            "straight edge should yield ~no corners, got {kept}/{}",
            stream.len()
        );
    }
}