summed_area/
lib.rs

1#![doc = include_str!("../README.md")]
2
3use std::ops::{Add, Sub, RangeBounds};
4pub use imgref::ImgRef;
5
6/// AKA Integral Image. Precomputed sums of subsections of a 2d array.
7///
8/// The implementation is generic and works with anything that can be added and subtracted (`f32`, `u64`, etc.).
9///
10/// The `SumType` should be a larger type for storing the sum, and `InputType` is a type of a slice.
11pub struct SummedArea<SumType> {
12    stride: usize,
13    sums: Vec<SumType>,
14}
15
16impl<SumType> SummedArea<SumType> where SumType: Default + Copy + Add<Output=SumType> + Sub<Output=SumType> {
17    /// Sum area in the slice. Height is implied from slice's length (`len/width`).
18    ///
19    /// Call it like `SummedArea::<u64>::new_slice::<u8>(byte_slice, width)` to specify a larger type to sum into (e.g. sum `u8`s into `u64`, so that they don't overflow).
20    #[inline]
21    #[must_use]
22    pub fn new_slice<InputType>(input: &[InputType], width: usize) -> Self where SumType: From<InputType>, InputType: Copy {
23        let height = input.len()/width;
24        Self::new(ImgRef::new(&input, width, height))
25    }
26
27    /// Sum the 2d rect. See [`ImgRef::new`].
28    #[must_use]
29    pub fn new<InputType>(input: ImgRef<InputType>) -> Self where SumType: From<InputType>, InputType: Copy {
30        let out_width = input.width()+1;
31        let out_height = input.height()+1;
32        let area = out_width.checked_mul(out_height).unwrap();
33
34        let mut sums: Vec<SumType> = Vec::with_capacity(area);
35        sums.resize_with(out_width, SumType::default); // first row of 0s to avoid check on y-1
36
37        let out = &mut sums.spare_capacity_mut()[..area-out_width];
38
39        let mut rows = input.rows().zip(out.chunks_exact_mut(out_width));
40        if let Some((in_row, out_row)) = rows.next() {
41            let (first, out_row) = out_row.split_first_mut().unwrap();
42            debug_assert_eq!(out_row.len(), input.width());
43            first.write(SumType::default()); // avoid check on x-1
44            let mut row_sum = SumType::default();
45            in_row.iter().copied().zip(out_row.iter_mut()).for_each(|(curr, out_col)| {
46                let curr = SumType::from(curr);
47                row_sum = row_sum + curr;
48                out_col.write(row_sum);
49            });
50            let mut prev_out_row = out_row;
51            rows.for_each(|(in_row, out_row)| {
52                let (first, out_row) = out_row.split_first_mut().unwrap();
53                debug_assert_eq!(out_row.len(), input.width());
54                first.write(SumType::default()); // avoid check on x-1
55                let mut up_left = SumType::default();
56                let mut left = SumType::default();
57                let cols = in_row.iter().copied().zip(out_row.iter_mut().zip(prev_out_row.iter_mut()));
58                cols.for_each(|(curr, (out_col, prev_out_px))| {
59                    // safety: previous row has been written to entirely. Use slice_assume_init_ref when stable.
60                    let up = unsafe { prev_out_px.assume_init_read() };
61                    let curr = SumType::from(curr);
62                    let curr_out = curr + up + left - up_left;
63                    out_col.write(curr_out);
64                    left = curr_out;
65                    up_left = up;
66                });
67                prev_out_row = out_row;
68            });
69        }
70        // safety: rows from chunks exact cover this entire area
71        unsafe {
72            sums.set_len(area);
73        }
74        debug_assert_eq!(sums.len()/out_width, out_height);
75        Self {
76            stride: out_width,
77            sums,
78        }
79    }
80
81    /// Width of the rect that has been summed up
82    #[inline]
83    #[must_use]
84    pub fn width(&self) -> usize {
85        self.stride-1
86    }
87
88    /// Height of the rect that has been summed up
89    #[inline]
90    #[must_use]
91    pub fn height(&self) -> usize {
92        self.sums.len() / self.stride - 1
93    }
94
95    /// Sum of all values in a rect at (x,y) being (width,height) elements large.
96    ///
97    /// If x+width or y+height are out of bounds it may return numerically incorrect result or panic.
98    #[inline(always)]
99    #[track_caller]
100    #[must_use]
101    pub fn sum_rect(&self, x: usize, y: usize, width: usize, height: usize) -> SumType {
102        self.try_sum_rect(x, y, width, height).expect("oob")
103    }
104
105    /// Sum of all values in a rect at (x,y) being (width,height) elements large.
106    ///
107    /// If x+width or y+height are out of bounds it may return numerically incorrect result or `None`.
108    #[inline]
109    #[must_use]
110    pub fn try_sum_rect(&self, x1: usize, y1: usize, width: usize, height: usize) -> Option<SumType> {
111        self.try_sum_bounds(x1, y1, x1+width, y1+height)
112    }
113
114    /// Sum of all values in a rect spanning `horizontal` colums and `vertical` columns (indexed from 0).
115    ///
116    /// If ranges are out of bounds it may return numerically incorrect result or panic.
117    #[inline(always)]
118    #[track_caller]
119    #[must_use]
120    pub fn sum_range(&self, horizontal: impl RangeBounds<usize>, vertical: impl RangeBounds<usize>) -> SumType {
121        self.try_sum_range(horizontal, vertical).expect("oob")
122    }
123
124    /// Sum of all values in a rect spanning `horizontal` colums and `vertical` columns (indexed from 0).
125    ///
126    /// If ranges are out of bounds it may return numerically incorrect result or `None`.
127    #[inline]
128    #[must_use]
129    pub fn try_sum_range(&self, horizontal: impl RangeBounds<usize>, vertical: impl RangeBounds<usize>) -> Option<SumType> {
130        let (x1, x2) = bounds(horizontal, self.stride);
131        let (y1, y2) = bounds(vertical, self.height());
132        self.try_sum_bounds(x1, y1, x2, y2)
133    }
134
135    #[inline]
136    fn try_sum_bounds(&self, x1: usize, y1: usize, x2: usize, y2: usize) -> Option<SumType> {
137        let x1y1 = x1 + y1 * self.stride;
138        let x1y2 = x1 + y2 * self.stride;
139        let x2y1 = x2 + y1 * self.stride;
140        let x2y2 = x2 + y2 * self.stride;
141        // there's no way to account for all possible numeric overflows in less than 4 branches
142        let tl = *self.sums.get(x1y1)?;
143        let tr = *self.sums.get(x2y1)?;
144        let bl = *self.sums.get(x1y2)?;
145        let br = *self.sums.get(x2y2)?;
146        Some(tl + br - tr - bl)
147    }
148}
149
150fn bounds(range: impl RangeBounds<usize>, max: usize) -> (usize, usize) {
151    let start = match range.start_bound() {
152        std::ops::Bound::Included(x) => *x,
153        std::ops::Bound::Unbounded => 0,
154        std::ops::Bound::Excluded(x) => *x+1,
155    };
156    let end = match range.end_bound() {
157        std::ops::Bound::Included(x) => *x+1,
158        std::ops::Bound::Excluded(x) => *x,
159        std::ops::Bound::Unbounded => max,
160    };
161    (start, end)
162}
163
164#[test]
165fn bounds_test() {
166    assert_eq!((0, 99), bounds(0.., 99));
167    assert_eq!((0, 0), bounds(0..0, 99));
168    assert_eq!((0, 1), bounds(0..1, 99));
169    assert_eq!((0, 2), bounds(0..=1, 99));
170    assert_eq!((0, 1), bounds(0..=0, 99));
171    assert_eq!((0, 50), bounds(..50, 99));
172    assert_eq!((0, 50), bounds(..=49, 99));
173}
174
175#[test]
176fn wiki() {
177    let s = SummedArea::new_slice(&[
178        31.,2.,4.,33.,5.,36.,
179        12.,26.,9.,10.,29.,25.,
180        13.,17.,21.,22.,20.,18.,
181        24.,23.,15.,16.,14.,19.,
182        30.,8.,28.,27.,11.,7.,
183        1.,35.,34.,3.,32.,6.,
184    ], 6);
185    assert_eq!(s.sums.len(), 7*7);
186    assert_eq!(s.width(), 6);
187    assert_eq!(s.height(), 6);
188    assert_eq!([
189        0., 0.,0.,0.,0.,0.,0.,
190        0., 31.,33.,37.,70.,75.,111.,
191        0., 43.,71.,84.,127.,161.,222.,
192        0., 56.,101.,135.,200.,254.,333.,
193        0., 80.,148.,197.,278.,346.,444.,
194        0., 110.,186.,263.,371.,450.,555.,
195        0., 111.,222.,333.,444.,555.,666.,
196    ], s.sums.as_slice());
197
198    assert_eq!(111., s.sum_rect(2, 3, 3, 2));
199    assert_eq!(111., s.try_sum_range(2..=4, 3..5).unwrap());
200}
201
202#[test]
203fn pixels() {
204    let s = SummedArea::<rgb::RGB<u16>>::new_slice(&[
205        rgb::RGB::<u8>::new(1,2,3),
206        rgb::RGB::<u8>::new(4,5,6),
207    ], 1);
208    assert_eq!(2, s.height());
209    assert_eq!(rgb::RGB::<u16>::new(1,2,3), s.sum_range(0..1, 0..1));
210    assert_eq!(rgb::RGB::new(4,5,6), s.sum_range(0..1, 1..2));
211    assert_eq!(rgb::RGB::new(5,7,9), s.sum_range(0..1, 0..2));
212}
213
214#[test]
215fn ones() {
216    let _ = SummedArea::<i64>::new_slice(&[1i32], 1);
217
218    let s = SummedArea::new_slice(&[
219        1.,1.,1.,
220        1.,1.,1.,
221        1.,1.,1.,
222    ], 3);
223    assert_eq!(s.height(), 3);
224    assert_eq!([
225        0.0, 0.0, 0.0, 0.0,
226        0.0, 1.0, 2.0, 3.0,
227        0.0, 2.0, 4.0, 6.0,
228        0.0, 3.0, 6.0, 9.0,
229    ], s.sums.as_slice());
230
231    assert_eq!(0., s.sum_rect(0,0,0,0));
232    assert_eq!(0., s.sum_rect(2,2,0,0));
233    assert_eq!(0., s.sum_rect(0,0,1,0));
234    assert_eq!(0., s.sum_rect(0,0,0,1));
235    assert_eq!(1., s.sum_rect(0,0,1,1));
236    assert_eq!(1., s.sum_rect(1,1,1,1));
237    assert_eq!(2., s.sum_rect(0,0,1,2));
238    assert_eq!(2., s.sum_rect(0,0,2,1));
239    assert_eq!(4., s.sum_rect(0,0,2,2));
240    assert_eq!(4., s.sum_rect(1,1,2,2));
241    assert_eq!(6., s.sum_rect(0,0,2,3));
242    assert_eq!(6., s.sum_rect(0,0,3,2));
243    assert_eq!(6., s.sum_rect(0,1,3,2));
244    assert_eq!(9., s.sum_rect(0,0,3,3));
245    assert_eq!(1., s.sum_rect(2,2,1,1));
246}
summed_area/lib.rs

summed_area/
lib.rs
