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

chess_corners_core/detect/radon/
response.rs
