chess_corners_core/detect/radon/
response.rs

1//! Whole-image Duda-Frese Radon response computation.
2//!
3//! Computes the dense `R(x, y) = (max_α S_α − min_α S_α)²` response
4//! over four ray angles using summed-area tables for `O(1)` per-pixel
5//! ray sums. Optional 2× bilinear upsampling at the entry boosts
6//! accuracy on small chessboard cells.
7//!
8//! # SAT element type
9//!
10//! Summed-area tables default to `i64`, which is always safe. Enable
11//! the `radon-sat-u32` crate feature to switch to `u32`, which halves
12//! SAT memory and widens SIMD lanes at the cost of a ~16 MP image-size
13//! cap (`255 · W · H ≤ u32::MAX`).
14
15use serde::{Deserialize, Serialize};
16
17#[cfg(feature = "rayon")]
18use rayon::prelude::*;
19
20#[cfg(all(feature = "simd", not(feature = "radon-sat-u32")))]
21use core::simd::Simd;
22
23#[cfg(all(feature = "simd", not(feature = "radon-sat-u32")))]
24use std::simd::cmp::SimdOrd;
25
26use super::primitives::{box_blur_inplace, PeakFitMode};
27use crate::ResponseMap;
28
29/// Number of pixels processed per SIMD iteration in
30/// `compute_response_row_simd`. Eight `i64` lanes is the natural
31/// width on AVX-512 / NEON-pair / SVE machines; smaller-width
32/// architectures fall back to the scalar tail handler in the same
33/// kernel.
34#[cfg(all(feature = "simd", not(feature = "radon-sat-u32")))]
35const RADON_LANES: usize = 8;
36
37/// Summed-area-table element type. Gated by the `radon-sat-u32`
38/// crate feature.
39#[cfg(not(feature = "radon-sat-u32"))]
40pub type SatElem = i64;
41
42/// Summed-area-table element type (feature `radon-sat-u32`).
43#[cfg(feature = "radon-sat-u32")]
44pub type SatElem = u32;
45
46/// Configuration for the whole-image Radon detector.
47#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
48#[serde(default)]
49pub struct RadonDetectorParams {
50    /// Half-length of each ray in **working-resolution** pixels (i.e.
51    /// post-upsample). The ray has `2·ray_radius + 1` samples. Paper
52    /// default at `image_upsample=2` is 4 working pixels ⇒ 2 physical.
53    pub ray_radius: u32,
54    /// Image-level supersampling factor. `1` operates on the input
55    /// pixel grid; `2` bilinearly upsamples first (paper default).
56    /// M1 supports the set `{1, 2}`; values `>= 3` are clamped to `2`
57    /// (see [`MAX_IMAGE_UPSAMPLE`]). Higher factors are future work.
58    pub image_upsample: u32,
59    /// Half-size of the box blur applied to the response map. `0`
60    /// disables blurring; `1` yields a 3×3 box.
61    pub response_blur_radius: u32,
62    /// Peak-fit mode for the 3-point subpixel refinement. Gaussian
63    /// (log-space) is more stable at near-plateau peaks.
64    pub peak_fit: PeakFitMode,
65    /// Relative response threshold as a fraction of the map's max
66    /// value. Used when `threshold_abs` is `None`.
67    pub threshold_rel: f32,
68    /// Absolute response threshold. Overrides `threshold_rel` when set.
69    /// The paper's `(max−min)²` response is always ≥ 0, so the strict
70    /// inequality `R > 0` that ChESS uses is not by itself selective
71    /// enough — use a positive absolute floor in practice.
72    pub threshold_abs: Option<f32>,
73    /// Non-maximum-suppression half-radius (in **working-resolution**
74    /// pixels).
75    pub nms_radius: u32,
76    /// Minimum count of positive-response neighbours in the NMS window
77    /// required to accept a peak. Rejects isolated noise.
78    pub min_cluster_size: u32,
79}
80
81impl Default for RadonDetectorParams {
82    fn default() -> Self {
83        Self {
84            ray_radius: 4,
85            image_upsample: 2,
86            response_blur_radius: 1,
87            peak_fit: PeakFitMode::Gaussian,
88            threshold_rel: 0.01,
89            threshold_abs: None,
90            nms_radius: 4,
91            min_cluster_size: 2,
92        }
93    }
94}
95
96/// Supported image-upsample factors: `{1, 2}`. Anything higher would
97/// need a different upsampler; values `>= 3` are clamped to `2` at the
98/// entry points rather than silently producing mismatched buffer sizes
99/// downstream.
100pub const MAX_IMAGE_UPSAMPLE: u32 = 2;
101
102impl RadonDetectorParams {
103    /// Clamp `image_upsample` into the supported set `{1, 2}`.
104    /// Values outside that range are silently clamped — callers can
105    /// detect truncation by comparing against [`MAX_IMAGE_UPSAMPLE`].
106    #[inline]
107    pub(crate) fn image_upsample_clamped(&self) -> u32 {
108        self.image_upsample.clamp(1, MAX_IMAGE_UPSAMPLE)
109    }
110
111    #[inline]
112    pub(crate) fn ray_radius_clamped(&self) -> u32 {
113        self.ray_radius.max(1)
114    }
115}
116
117/// Reusable scratch for the whole-image Radon detector. Holds the
118/// upsampled image buffer, the four summed-area tables, the response
119/// map, and the box-blur scratch. All buffers grow on demand and are
120/// reused across frames — same pattern as `PyramidBuffers`.
121#[derive(Debug, Default)]
122pub struct RadonBuffers {
123    upsampled: Vec<u8>,
124    working_w: usize,
125    working_h: usize,
126    row_cumsum: Vec<SatElem>,
127    col_cumsum: Vec<SatElem>,
128    diag_pos_cumsum: Vec<SatElem>,
129    diag_neg_cumsum: Vec<SatElem>,
130    response: Vec<f32>,
131    blur_scratch: Vec<f32>,
132}
133
134impl RadonBuffers {
135    /// Create an empty set of buffers. They grow on first use.
136    pub fn new() -> Self {
137        Self::default()
138    }
139
140    fn ensure_capacity(&mut self, input_w: usize, input_h: usize, upsample: u32) {
141        let up = upsample.max(1) as usize;
142        let ww = input_w * up;
143        let wh = input_h * up;
144        let n = ww * wh;
145        self.working_w = ww;
146        self.working_h = wh;
147        if up > 1 {
148            self.upsampled.resize(n, 0);
149        } else {
150            self.upsampled.clear();
151        }
152        self.row_cumsum.resize(n, SatElem::default());
153        self.col_cumsum.resize(n, SatElem::default());
154        self.diag_pos_cumsum.resize(n, SatElem::default());
155        self.diag_neg_cumsum.resize(n, SatElem::default());
156        self.response.resize(n, 0.0);
157        self.blur_scratch.resize(n, 0.0);
158    }
159}
160
161fn upsample_bilinear_2x(src: &[u8], w: usize, h: usize, out: &mut [u8]) {
162    debug_assert_eq!(src.len(), w * h);
163    debug_assert_eq!(out.len(), 4 * w * h);
164    if w == 0 || h == 0 {
165        return;
166    }
167    let ww = 2 * w;
168
169    let row_kernel = |iy: usize, dst: &mut [u8]| {
170        let sy = iy as f32 * 0.5;
171        let y0f = sy.floor();
172        let y0 = (y0f as isize).max(0) as usize;
173        let y1 = (y0 + 1).min(h - 1);
174        let ty = (sy - y0f).clamp(0.0, 1.0);
175        for (ix, out_px) in dst.iter_mut().enumerate() {
176            let sx = ix as f32 * 0.5;
177            let x0f = sx.floor();
178            let x0 = (x0f as isize).max(0) as usize;
179            let x1 = (x0 + 1).min(w - 1);
180            let tx = (sx - x0f).clamp(0.0, 1.0);
181            let i00 = src[y0 * w + x0] as f32;
182            let i10 = src[y0 * w + x1] as f32;
183            let i01 = src[y1 * w + x0] as f32;
184            let i11 = src[y1 * w + x1] as f32;
185            let a = i00 + (i10 - i00) * tx;
186            let b = i01 + (i11 - i01) * tx;
187            let v = a + (b - a) * ty;
188            *out_px = v.round().clamp(0.0, 255.0) as u8;
189        }
190    };
191
192    #[cfg(feature = "rayon")]
193    {
194        out.par_chunks_mut(ww)
195            .enumerate()
196            .for_each(|(iy, row)| row_kernel(iy, row));
197    }
198    #[cfg(not(feature = "rayon"))]
199    {
200        for (iy, row) in out.chunks_mut(ww).enumerate() {
201            row_kernel(iy, row);
202        }
203    }
204}
205
206#[inline]
207fn sat_row(img: &[u8], w: usize, h: usize, row_cumsum: &mut [SatElem]) {
208    debug_assert_eq!(row_cumsum.len(), w * h);
209    for y in 0..h {
210        let mut acc: SatElem = SatElem::default();
211        for x in 0..w {
212            acc += SatElem::from(img[y * w + x]);
213            row_cumsum[y * w + x] = acc;
214        }
215    }
216}
217
218#[inline]
219fn sat_col(img: &[u8], w: usize, h: usize, col_cumsum: &mut [SatElem]) {
220    debug_assert_eq!(col_cumsum.len(), w * h);
221    for x in 0..w {
222        let mut acc: SatElem = SatElem::default();
223        for y in 0..h {
224            acc += SatElem::from(img[y * w + x]);
225            col_cumsum[y * w + x] = acc;
226        }
227    }
228}
229
230#[inline]
231fn sat_diag_pos(img: &[u8], w: usize, h: usize, diag_pos_cumsum: &mut [SatElem]) {
232    debug_assert_eq!(diag_pos_cumsum.len(), w * h);
233    for y in 0..h {
234        for x in 0..w {
235            let prev = if y > 0 && x > 0 {
236                diag_pos_cumsum[(y - 1) * w + (x - 1)]
237            } else {
238                SatElem::default()
239            };
240            diag_pos_cumsum[y * w + x] = prev + SatElem::from(img[y * w + x]);
241        }
242    }
243}
244
245#[inline]
246fn sat_diag_neg(img: &[u8], w: usize, h: usize, diag_neg_cumsum: &mut [SatElem]) {
247    debug_assert_eq!(diag_neg_cumsum.len(), w * h);
248    for y in 0..h {
249        for x in 0..w {
250            let prev = if y > 0 && x + 1 < w {
251                diag_neg_cumsum[(y - 1) * w + (x + 1)]
252            } else {
253                SatElem::default()
254            };
255            diag_neg_cumsum[y * w + x] = prev + SatElem::from(img[y * w + x]);
256        }
257    }
258}
259
260fn build_cumsums(
261    img: &[u8],
262    w: usize,
263    h: usize,
264    row_cumsum: &mut [SatElem],
265    col_cumsum: &mut [SatElem],
266    diag_pos_cumsum: &mut [SatElem],
267    diag_neg_cumsum: &mut [SatElem],
268) {
269    debug_assert_eq!(img.len(), w * h);
270
271    #[cfg(feature = "rayon")]
272    {
273        rayon::join(
274            || {
275                rayon::join(
276                    || sat_row(img, w, h, row_cumsum),
277                    || sat_col(img, w, h, col_cumsum),
278                );
279            },
280            || {
281                rayon::join(
282                    || sat_diag_pos(img, w, h, diag_pos_cumsum),
283                    || sat_diag_neg(img, w, h, diag_neg_cumsum),
284                );
285            },
286        );
287    }
288    #[cfg(not(feature = "rayon"))]
289    {
290        sat_row(img, w, h, row_cumsum);
291        sat_col(img, w, h, col_cumsum);
292        sat_diag_pos(img, w, h, diag_pos_cumsum);
293        sat_diag_neg(img, w, h, diag_neg_cumsum);
294    }
295}
296
297/// Bundle of cumsum tables + geometry, passed to response kernels.
298struct Cumsums<'a> {
299    row: &'a [SatElem],
300    col: &'a [SatElem],
301    diag_pos: &'a [SatElem],
302    diag_neg: &'a [SatElem],
303    w: usize,
304    h: usize,
305}
306
307#[inline(always)]
308fn radon_response_at(cs: &Cumsums<'_>, r: usize, x: usize, y: usize) -> f32 {
309    let w = cs.w;
310    let s_h_hi = cs.row[y * w + (x + r)];
311    let s_h_lo = if x > r {
312        cs.row[y * w + (x - r - 1)]
313    } else {
314        SatElem::default()
315    };
316    let s_h = s_h_hi - s_h_lo;
317
318    let s_v_hi = cs.col[(y + r) * w + x];
319    let s_v_lo = if y > r {
320        cs.col[(y - r - 1) * w + x]
321    } else {
322        SatElem::default()
323    };
324    let s_v = s_v_hi - s_v_lo;
325
326    let s_d1_hi = cs.diag_pos[(y + r) * w + (x + r)];
327    let s_d1_lo = if x > r && y > r {
328        cs.diag_pos[(y - r - 1) * w + (x - r - 1)]
329    } else {
330        SatElem::default()
331    };
332    let s_d1 = s_d1_hi - s_d1_lo;
333
334    let s_d2_hi = cs.diag_neg[(y + r) * w + (x - r)];
335    let s_d2_lo = if y > r && x + r + 1 < w {
336        cs.diag_neg[(y - r - 1) * w + (x + r + 1)]
337    } else {
338        SatElem::default()
339    };
340    let s_d2 = s_d2_hi - s_d2_lo;
341
342    let s = [s_h, s_v, s_d1, s_d2];
343    let (mut mx, mut mn) = (s[0], s[0]);
344    for &v in &s[1..] {
345        if v > mx {
346            mx = v;
347        }
348        if v < mn {
349            mn = v;
350        }
351    }
352    let d = sat_to_f32(mx - mn);
353    d * d
354}
355
356#[inline]
357fn compute_response_row(cs: &Cumsums<'_>, ray_radius: usize, y: usize, row: &mut [f32]) {
358    let w = cs.w;
359    let h = cs.h;
360    let r = ray_radius;
361    if y < r || y + r >= h {
362        for v in row.iter_mut() {
363            *v = 0.0;
364        }
365        return;
366    }
367    for v in row[..r].iter_mut() {
368        *v = 0.0;
369    }
370
371    #[cfg(all(feature = "simd", not(feature = "radon-sat-u32")))]
372    {
373        compute_response_row_simd(cs, r, y, row);
374    }
375    #[cfg(not(all(feature = "simd", not(feature = "radon-sat-u32"))))]
376    {
377        for (x, out_px) in row.iter_mut().enumerate().take(w - r).skip(r) {
378            *out_px = radon_response_at(cs, r, x, y);
379        }
380    }
381
382    for v in row[(w - r)..].iter_mut() {
383        *v = 0.0;
384    }
385}
386
387#[cfg(all(feature = "simd", not(feature = "radon-sat-u32")))]
388#[inline]
389fn compute_response_row_simd(cs: &Cumsums<'_>, r: usize, y: usize, row: &mut [f32]) {
390    type S = Simd<i64, RADON_LANES>;
391
392    let w = cs.w;
393
394    row[r] = radon_response_at(cs, r, r, y);
395
396    let interior_start = r + 1;
397    let interior_end = w - r;
398    let mut x = interior_start;
399
400    if y <= r {
401        for (x, cell) in row
402            .iter_mut()
403            .enumerate()
404            .take(interior_end)
405            .skip(interior_start)
406        {
407            *cell = radon_response_at(cs, r, x, y);
408        }
409        return;
410    }
411
412    let h_row_base = y * w;
413    let v_hi_base = (y + r) * w;
414    let v_lo_base = (y - r - 1) * w;
415    let d1_hi_base = (y + r) * w;
416    let d1_lo_base = (y - r - 1) * w;
417    let d2_hi_base = (y + r) * w;
418    let d2_lo_base = (y - r - 1) * w;
419
420    while x + RADON_LANES < interior_end {
421        let s_h_hi = S::from_slice(&cs.row[h_row_base + x + r..h_row_base + x + r + RADON_LANES]);
422        let s_h_lo =
423            S::from_slice(&cs.row[h_row_base + x - r - 1..h_row_base + x - r - 1 + RADON_LANES]);
424        let s_h = s_h_hi - s_h_lo;
425
426        let s_v_hi = S::from_slice(&cs.col[v_hi_base + x..v_hi_base + x + RADON_LANES]);
427        let s_v_lo = S::from_slice(&cs.col[v_lo_base + x..v_lo_base + x + RADON_LANES]);
428        let s_v = s_v_hi - s_v_lo;
429
430        let s_d1_hi =
431            S::from_slice(&cs.diag_pos[d1_hi_base + x + r..d1_hi_base + x + r + RADON_LANES]);
432        let s_d1_lo = S::from_slice(
433            &cs.diag_pos[d1_lo_base + x - r - 1..d1_lo_base + x - r - 1 + RADON_LANES],
434        );
435        let s_d1 = s_d1_hi - s_d1_lo;
436
437        let s_d2_hi =
438            S::from_slice(&cs.diag_neg[d2_hi_base + x - r..d2_hi_base + x - r + RADON_LANES]);
439        let s_d2_lo = S::from_slice(
440            &cs.diag_neg[d2_lo_base + x + r + 1..d2_lo_base + x + r + 1 + RADON_LANES],
441        );
442        let s_d2 = s_d2_hi - s_d2_lo;
443
444        let mx = s_h.simd_max(s_v).simd_max(s_d1.simd_max(s_d2));
445        let mn = s_h.simd_min(s_v).simd_min(s_d1.simd_min(s_d2));
446        let diff = mx - mn;
447        let arr = diff.to_array();
448        for (lane, v) in arr.iter().enumerate() {
449            let f = *v as f32;
450            row[x + lane] = f * f;
451        }
452
453        x += RADON_LANES;
454    }
455
456    while x < interior_end {
457        row[x] = radon_response_at(cs, r, x, y);
458        x += 1;
459    }
460}
461
462fn compute_response(cs: &Cumsums<'_>, ray_radius: usize, out: &mut [f32]) {
463    let w = cs.w;
464    let h = cs.h;
465    debug_assert_eq!(out.len(), w * h);
466    if w <= 2 * ray_radius || h <= 2 * ray_radius {
467        out.fill(0.0);
468        return;
469    }
470
471    #[cfg(feature = "rayon")]
472    {
473        out.par_chunks_mut(w).enumerate().for_each(|(y, row)| {
474            compute_response_row(cs, ray_radius, y, row);
475        });
476    }
477    #[cfg(not(feature = "rayon"))]
478    {
479        for (y, row) in out.chunks_mut(w).enumerate() {
480            compute_response_row(cs, ray_radius, y, row);
481        }
482    }
483}
484
485#[inline]
486fn sat_to_f32(v: SatElem) -> f32 {
487    v as f32
488}
489
490/// Compute the dense Radon response into `buffers.response` and return
491/// a read-only [`RadonResponseView`] at **working resolution** (i.e.
492/// `input_dim × image_upsample`).
493pub fn radon_response_u8<'a>(
494    img: &[u8],
495    w: usize,
496    h: usize,
497    params: &RadonDetectorParams,
498    buffers: &'a mut RadonBuffers,
499) -> RadonResponseView<'a> {
500    assert_eq!(img.len(), w * h, "img len must equal w*h");
501    let up = params.image_upsample_clamped();
502    buffers.ensure_capacity(w, h, up);
503    let ww = buffers.working_w;
504    let wh = buffers.working_h;
505
506    #[cfg(feature = "radon-sat-u32")]
507    {
508        let pixels = (ww as u64) * (wh as u64);
509        let max_sum = 255u64.checked_mul(pixels);
510        assert!(
511            matches!(max_sum, Some(v) if v <= u32::MAX as u64),
512            "radon-sat-u32: 255*W*H ({}*{}) exceeds u32::MAX; \
513             either rebuild without the radon-sat-u32 feature or \
514             downsample the input",
515            ww,
516            wh,
517        );
518    }
519
520    let working_img: &[u8] = if up > 1 {
521        upsample_bilinear_2x_if_needed(img, w, h, up, &mut buffers.upsampled);
522        &buffers.upsampled
523    } else {
524        img
525    };
526
527    build_cumsums(
528        working_img,
529        ww,
530        wh,
531        &mut buffers.row_cumsum,
532        &mut buffers.col_cumsum,
533        &mut buffers.diag_pos_cumsum,
534        &mut buffers.diag_neg_cumsum,
535    );
536
537    let cs = Cumsums {
538        row: &buffers.row_cumsum,
539        col: &buffers.col_cumsum,
540        diag_pos: &buffers.diag_pos_cumsum,
541        diag_neg: &buffers.diag_neg_cumsum,
542        w: ww,
543        h: wh,
544    };
545    compute_response(
546        &cs,
547        params.ray_radius_clamped() as usize,
548        &mut buffers.response,
549    );
550
551    box_blur_inplace(
552        &mut buffers.response,
553        &mut buffers.blur_scratch,
554        ww,
555        wh,
556        params.response_blur_radius as usize,
557    );
558    RadonResponseView {
559        data: &buffers.response,
560        w: ww,
561        h: wh,
562    }
563}
564
565/// Borrow of the dense working-resolution response map. Cheaply
566/// convertible to a [`ResponseMap`] via [`Self::to_response_map`]
567/// when ownership is required (e.g. for the classic
568/// [`detect_corners_from_response`](crate::detect::detect_corners_from_response)).
569#[derive(Debug)]
570pub struct RadonResponseView<'a> {
571    pub(super) data: &'a [f32],
572    pub(super) w: usize,
573    pub(super) h: usize,
574}
575
576impl<'a> RadonResponseView<'a> {
577    #[inline]
578    pub fn width(&self) -> usize {
579        self.w
580    }
581
582    #[inline]
583    pub fn height(&self) -> usize {
584        self.h
585    }
586
587    #[inline]
588    pub fn data(&self) -> &[f32] {
589        self.data
590    }
591
592    #[inline]
593    pub fn at(&self, x: usize, y: usize) -> f32 {
594        self.data[y * self.w + x]
595    }
596
597    pub fn to_response_map(&self) -> ResponseMap {
598        ResponseMap {
599            w: self.w,
600            h: self.h,
601            data: self.data.to_vec(),
602        }
603    }
604}
605
606#[inline]
607fn upsample_bilinear_2x_if_needed(img: &[u8], w: usize, h: usize, up: u32, out: &mut Vec<u8>) {
608    debug_assert_eq!(up, 2, "image_upsample must be 1 or 2");
609    let _ = up;
610    out.resize(4 * w * h, 0);
611    upsample_bilinear_2x(img, w, h, out);
612}
613
614#[cfg(test)]
615mod tests {
616    use super::super::test_fixtures::synthetic_chessboard_aa;
617    use super::*;
618
619    #[test]
620    fn row_cumsum_matches_naive_sum() {
621        let w = 5usize;
622        let h = 3usize;
623        let img: Vec<u8> = (0..(w * h) as u8).collect();
624        let mut r = vec![SatElem::default(); w * h];
625        let mut c = vec![SatElem::default(); w * h];
626        let mut d1 = vec![SatElem::default(); w * h];
627        let mut d2 = vec![SatElem::default(); w * h];
628        build_cumsums(&img, w, h, &mut r, &mut c, &mut d1, &mut d2);
629        for y in 0..h {
630            let mut expected: SatElem = SatElem::default();
631            for x in 0..w {
632                expected += SatElem::from(img[y * w + x]);
633                assert_eq!(r[y * w + x], expected);
634            }
635        }
636    }
637
638    #[test]
639    fn col_cumsum_matches_naive_sum() {
640        let w = 4usize;
641        let h = 5usize;
642        let img: Vec<u8> = (0..(w * h) as u8).collect();
643        let mut r = vec![SatElem::default(); w * h];
644        let mut c = vec![SatElem::default(); w * h];
645        let mut d1 = vec![SatElem::default(); w * h];
646        let mut d2 = vec![SatElem::default(); w * h];
647        build_cumsums(&img, w, h, &mut r, &mut c, &mut d1, &mut d2);
648        for x in 0..w {
649            let mut expected: SatElem = SatElem::default();
650            for y in 0..h {
651                expected += SatElem::from(img[y * w + x]);
652                assert_eq!(c[y * w + x], expected);
653            }
654        }
655    }
656
657    #[test]
658    fn diag_pos_cumsum_matches_naive_sum() {
659        let w = 4usize;
660        let h = 4usize;
661        let img: Vec<u8> = (0..(w * h) as u8).collect();
662        let mut r = vec![SatElem::default(); w * h];
663        let mut c = vec![SatElem::default(); w * h];
664        let mut d1 = vec![SatElem::default(); w * h];
665        let mut d2 = vec![SatElem::default(); w * h];
666        build_cumsums(&img, w, h, &mut r, &mut c, &mut d1, &mut d2);
667        for y in 0..h {
668            for x in 0..w {
669                let mut expected: SatElem = SatElem::default();
670                let mut xi = x as isize;
671                let mut yi = y as isize;
672                while xi >= 0 && yi >= 0 {
673                    expected += SatElem::from(img[yi as usize * w + xi as usize]);
674                    xi -= 1;
675                    yi -= 1;
676                }
677                assert_eq!(d1[y * w + x], expected, "at ({},{})", x, y);
678            }
679        }
680    }
681
682    #[test]
683    fn diag_neg_cumsum_matches_naive_sum() {
684        let w = 4usize;
685        let h = 4usize;
686        let img: Vec<u8> = (0..(w * h) as u8).collect();
687        let mut r = vec![SatElem::default(); w * h];
688        let mut c = vec![SatElem::default(); w * h];
689        let mut d1 = vec![SatElem::default(); w * h];
690        let mut d2 = vec![SatElem::default(); w * h];
691        build_cumsums(&img, w, h, &mut r, &mut c, &mut d1, &mut d2);
692        for y in 0..h {
693            for x in 0..w {
694                let mut expected: SatElem = SatElem::default();
695                let mut xi = x as isize;
696                let mut yi = y as isize;
697                while xi < w as isize && yi >= 0 {
698                    expected += SatElem::from(img[yi as usize * w + xi as usize]);
699                    xi += 1;
700                    yi -= 1;
701                }
702                assert_eq!(d2[y * w + x], expected, "at ({},{})", x, y);
703            }
704        }
705    }
706
707    #[test]
708    fn ray_sums_via_sat_match_direct_sums() {
709        let w = 15usize;
710        let h = 15usize;
711        let img: Vec<u8> = (0..(w * h)).map(|i| (i % 251) as u8).collect();
712        let mut r = vec![SatElem::default(); w * h];
713        let mut c = vec![SatElem::default(); w * h];
714        let mut d1 = vec![SatElem::default(); w * h];
715        let mut d2 = vec![SatElem::default(); w * h];
716        build_cumsums(&img, w, h, &mut r, &mut c, &mut d1, &mut d2);
717
718        let ray_r = 3usize;
719        for y in ray_r..(h - ray_r) {
720            for x in ray_r..(w - ray_r) {
721                let mut h_sum: SatElem = SatElem::default();
722                for k in 0..=(2 * ray_r) {
723                    let xx = x + k - ray_r;
724                    h_sum += SatElem::from(img[y * w + xx]);
725                }
726                let h_hi = r[y * w + (x + ray_r)];
727                let h_lo = if x > ray_r {
728                    r[y * w + (x - ray_r - 1)]
729                } else {
730                    SatElem::default()
731                };
732                assert_eq!(h_hi - h_lo, h_sum, "horiz at ({},{})", x, y);
733
734                let mut v_sum: SatElem = SatElem::default();
735                for k in 0..=(2 * ray_r) {
736                    let yy = y + k - ray_r;
737                    v_sum += SatElem::from(img[yy * w + x]);
738                }
739                let v_hi = c[(y + ray_r) * w + x];
740                let v_lo = if y > ray_r {
741                    c[(y - ray_r - 1) * w + x]
742                } else {
743                    SatElem::default()
744                };
745                assert_eq!(v_hi - v_lo, v_sum, "vert at ({},{})", x, y);
746
747                let mut d1_sum: SatElem = SatElem::default();
748                for k in 0..=(2 * ray_r) {
749                    let xx = x + k - ray_r;
750                    let yy = y + k - ray_r;
751                    d1_sum += SatElem::from(img[yy * w + xx]);
752                }
753                let d1_hi = d1[(y + ray_r) * w + (x + ray_r)];
754                let d1_lo = if x > ray_r && y > ray_r {
755                    d1[(y - ray_r - 1) * w + (x - ray_r - 1)]
756                } else {
757                    SatElem::default()
758                };
759                assert_eq!(d1_hi - d1_lo, d1_sum, "diag+ at ({},{})", x, y);
760
761                let mut d2_sum: SatElem = SatElem::default();
762                for k in 0..=(2 * ray_r) {
763                    let xx = x + ray_r - k;
764                    let yy = y + k - ray_r;
765                    d2_sum += SatElem::from(img[yy * w + xx]);
766                }
767                let d2_hi = d2[(y + ray_r) * w + (x - ray_r)];
768                let d2_lo = if y > ray_r && x + ray_r + 1 < w {
769                    d2[(y - ray_r - 1) * w + (x + ray_r + 1)]
770                } else {
771                    SatElem::default()
772                };
773                assert_eq!(d2_hi - d2_lo, d2_sum, "diag- at ({},{})", x, y);
774            }
775        }
776    }
777
778    #[test]
779    fn response_map_is_non_negative_everywhere() {
780        const SIZE: usize = 29;
781        const CELL: usize = 6;
782        let img = synthetic_chessboard_aa(SIZE, CELL, (14.2, 14.6), 30, 230);
783        let params = RadonDetectorParams {
784            image_upsample: 1,
785            ..RadonDetectorParams::default()
786        };
787        let mut buffers = RadonBuffers::new();
788        let resp = radon_response_u8(&img, SIZE, SIZE, &params, &mut buffers);
789        for &v in resp.data() {
790            assert!(v >= 0.0, "negative response value: {v}");
791        }
792    }
793
794    #[test]
795    fn image_upsample_above_cap_is_clamped_not_panicked() {
796        const SIZE: usize = 29;
797        const CELL: usize = 6;
798        let img = synthetic_chessboard_aa(SIZE, CELL, (14.2, 14.6), 30, 230);
799        let params = RadonDetectorParams {
800            image_upsample: 5,
801            ..RadonDetectorParams::default()
802        };
803        assert_eq!(params.image_upsample_clamped(), MAX_IMAGE_UPSAMPLE);
804
805        let mut buffers = RadonBuffers::new();
806        let resp = radon_response_u8(&img, SIZE, SIZE, &params, &mut buffers);
807        assert_eq!(resp.width(), SIZE * MAX_IMAGE_UPSAMPLE as usize);
808        assert_eq!(resp.height(), SIZE * MAX_IMAGE_UPSAMPLE as usize);
809    }
810
811    #[test]
812    fn image_upsample_zero_is_clamped_to_one() {
813        const SIZE: usize = 21;
814        const CELL: usize = 5;
815        let img = synthetic_chessboard_aa(SIZE, CELL, (10.1, 10.4), 30, 230);
816        let params = RadonDetectorParams {
817            image_upsample: 0,
818            ..RadonDetectorParams::default()
819        };
820        assert_eq!(params.image_upsample_clamped(), 1);
821
822        let mut buffers = RadonBuffers::new();
823        let resp = radon_response_u8(&img, SIZE, SIZE, &params, &mut buffers);
824        assert_eq!(resp.width(), SIZE);
825        assert_eq!(resp.height(), SIZE);
826    }
827
828    #[test]
829    fn radon_response_u8_handles_zero_extent_image() {
830        let img: Vec<u8> = Vec::new();
831        let params = RadonDetectorParams::default();
832        let mut buffers = RadonBuffers::new();
833        let resp = radon_response_u8(&img, 0, 0, &params, &mut buffers);
834        assert_eq!(resp.width(), 0);
835        assert_eq!(resp.height(), 0);
836        assert!(resp.data().is_empty());
837
838        let params_no_upsample = RadonDetectorParams {
839            image_upsample: 1,
840            ..RadonDetectorParams::default()
841        };
842        let resp = radon_response_u8(&img, 0, 0, &params_no_upsample, &mut buffers);
843        assert_eq!(resp.width(), 0);
844        assert_eq!(resp.height(), 0);
845    }
846
847    #[cfg(all(feature = "simd", not(feature = "radon-sat-u32")))]
848    #[test]
849    fn simd_kernel_matches_scalar_at_every_interior_pixel() {
850        for w in [16usize, 17, 18, 23, 24, 25, 32, 33] {
851            let h = 24usize;
852            let mut img = vec![0u8; w * h];
853            for (i, p) in img.iter_mut().enumerate() {
854                *p = ((i.wrapping_mul(37) ^ (i >> 3)) & 0xff) as u8;
855            }
856            let params = RadonDetectorParams {
857                image_upsample: 1,
858                response_blur_radius: 0,
859                threshold_abs: Some(0.0),
860                ..RadonDetectorParams::default()
861            };
862            let mut buffers = RadonBuffers::new();
863            let response_snapshot: Vec<f32> = {
864                let resp = radon_response_u8(&img, w, h, &params, &mut buffers);
865                resp.data().to_vec()
866            };
867
868            let r = params.ray_radius_clamped() as usize;
869            let cs = Cumsums {
870                row: &buffers.row_cumsum,
871                col: &buffers.col_cumsum,
872                diag_pos: &buffers.diag_pos_cumsum,
873                diag_neg: &buffers.diag_neg_cumsum,
874                w,
875                h,
876            };
877            for y in r..(h - r) {
878                for x in r..(w - r) {
879                    let expected = radon_response_at(&cs, r, x, y);
880                    let got = response_snapshot[y * w + x];
881                    assert!(
882                        (expected - got).abs() < 1e-3,
883                        "mismatch at w={w}, (x={x}, y={y}): scalar={expected}, simd={got}",
884                    );
885                }
886            }
887        }
888    }
889}
chess_corners_core/detect/radon/response.rs

chess_corners_core/detect/radon/
response.rs
