numeris 0.5.13

Pure-Rust numerical algorithms library — high performance with SIMD support while also supporting no-std for embedded and WASM targets.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

use crate::dynmatrix::DynMatrix;
use crate::traits::{FloatScalar, MatrixRef};

/// Compute the summed-area table (SAT) of `src`.
///
/// Returns a `(nrows + 1) × (ncols + 1)` matrix where
/// `sat[i+1, j+1] = Σ_{r ≤ i, c ≤ j} src[r, c]` and the first row and column
/// are zero. The zero-padded layout means rectangle sums over any region
/// `[r0, r1) × [c0, c1)` of the original image reduce to a single
/// four-term expression; see [`integral_rect_sum`].
///
/// Complexity: `O(H·W)` to construct, `O(1)` per rectangle query.
///
/// # Example
///
/// ```
/// use numeris::DynMatrix;
/// use numeris::imageproc::{integral_image, integral_rect_sum};
///
/// let img = DynMatrix::from_rows(
///     3, 3,
///     &[1.0_f64, 2.0, 3.0,
///       4.0,     5.0, 6.0,
///       7.0,     8.0, 9.0],
/// );
/// let sat = integral_image(&img);
/// // Sum over the whole 3×3 image.
/// assert!((integral_rect_sum(&sat, 0, 0, 3, 3) - 45.0).abs() < 1e-12);
/// // Sum over the centre pixel only.
/// assert!((integral_rect_sum(&sat, 1, 1, 2, 2) - 5.0).abs() < 1e-12);
/// ```
pub fn integral_image<T: FloatScalar>(src: &DynMatrix<T>) -> DynMatrix<T> {
    let h = src.nrows();
    let w = src.ncols();
    let mut sat = DynMatrix::<T>::zeros(h + 1, w + 1);

    // Running row sums, then cumulative column sums (two-pass SAT).
    for i in 0..h {
        let mut row_sum = T::zero();
        for j in 0..w {
            row_sum = row_sum + *src.get(i, j);
            // sat[i+1, j+1] = sat[i, j+1] + row_sum_up_to_j
            let above = sat[(i, j + 1)];
            sat[(i + 1, j + 1)] = above + row_sum;
        }
    }
    sat
}

/// Summed-area tables of `src` and of `src²`, built in a single pass.
///
/// Fuses the squaring into the SAT construction so callers (e.g.
/// `local_variance`) do not materialize a full squared copy of the image
/// first. Produces the same tables as `integral_image(src)` and
/// `integral_image(&src²)`, term for term.
pub(crate) fn integral_image_with_squares<T: FloatScalar>(
    src: &DynMatrix<T>,
) -> (DynMatrix<T>, DynMatrix<T>) {
    let h = src.nrows();
    let w = src.ncols();
    let mut sat = DynMatrix::<T>::zeros(h + 1, w + 1);
    let mut sat2 = DynMatrix::<T>::zeros(h + 1, w + 1);

    for i in 0..h {
        let mut row_sum = T::zero();
        let mut row_sum2 = T::zero();
        for j in 0..w {
            let v = *src.get(i, j);
            row_sum = row_sum + v;
            row_sum2 = row_sum2 + v * v;
            sat[(i + 1, j + 1)] = sat[(i, j + 1)] + row_sum;
            sat2[(i + 1, j + 1)] = sat2[(i, j + 1)] + row_sum2;
        }
    }
    (sat, sat2)
}

/// Sum over the half-open rectangle `[r0, r1) × [c0, c1)` of the original
/// image, in O(1), using the SAT produced by [`integral_image`].
///
/// Indices are into the *original* image (not the padded SAT); the `+1`
/// offset is applied internally. Out-of-range indices are clamped to the
/// image extent (so `integral_rect_sum(sat, 0, 0, nrows, ncols)` sums the
/// whole image).
///
/// # Panics
///
/// Panics (debug) if `r0 > r1` or `c0 > c1`.
#[inline]
pub fn integral_rect_sum<T: FloatScalar>(
    sat: &DynMatrix<T>,
    r0: usize,
    c0: usize,
    r1: usize,
    c1: usize,
) -> T {
    debug_assert!(r0 <= r1 && c0 <= c1, "rectangle bounds must be ordered");
    let nrows = sat.nrows().saturating_sub(1);
    let ncols = sat.ncols().saturating_sub(1);
    let r0 = r0.min(nrows);
    let c0 = c0.min(ncols);
    let r1 = r1.min(nrows);
    let c1 = c1.min(ncols);
    // Inclusion–exclusion on the SAT's one-based indices.
    sat[(r1, c1)] + sat[(r0, c0)] - sat[(r1, c0)] - sat[(r0, c1)]
}