ndarray_ndimage/filters/
min_max.rs

1use std::collections::VecDeque;
2
3use ndarray::{Array, Array1, ArrayBase, Axis, Data, Dimension, ScalarOperand, Zip};
4use num_traits::{FromPrimitive, Num};
5
6use crate::{array_like, filters::origin_check, pad_to, BorderMode};
7
8/// Calculate a 1-D maximum filter along the given axis.
9///
10/// The lines of the array along the given axis are filtered with a maximum filter of given size.
11///
12/// * `data` - The input N-D data.
13/// * `size` - Length along which to calculate 1D maximum.
14/// * `axis` - The axis of input along which to calculate.
15/// * `mode` - Method that will be used to select the padded values. See the
16///   [`CorrelateMode`](crate::CorrelateMode) enum for more information.
17/// * `origin` - Controls the placement of the filter on the input array’s pixels. A value of 0
18///   centers the filter over the pixel, with positive values shifting the filter to the left, and
19///   negative ones to the right.
20pub fn maximum_filter1d<S, A, D>(
21    data: &ArrayBase<S, D>,
22    size: usize,
23    axis: Axis,
24    mode: BorderMode<A>,
25    origin: isize,
26) -> Array<A, D>
27where
28    S: Data<Elem = A>,
29    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
30    D: Dimension,
31{
32    let mut output = data.to_owned();
33    maximum_filter1d_to(data, size, axis, mode, origin, &mut output);
34    output
35}
36
37/// Calculate a multidimensional maximum filter.
38///
39/// * `data` - The input N-D data.
40/// * `size` - Length along which to calculate 1D maximum.
41/// * `mode` - Method that will be used to select the padded values. See the
42///   [`CorrelateMode`](crate::CorrelateMode) enum for more information.
43/// * `origin` - Controls the placement of the filter on the input array’s pixels. A value of 0
44///   centers the filter over the pixel, with positive values shifting the filter to the left, and
45///   negative ones to the right.
46pub fn maximum_filter<S, A, D>(
47    data: &ArrayBase<S, D>,
48    size: usize,
49    mode: BorderMode<A>,
50    origin: isize,
51) -> Array<A, D>
52where
53    S: Data<Elem = A>,
54    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
55    D: Dimension,
56{
57    // We need 2 buffers because
58    // * We're reading neignbors so we can't read and write on the same location.
59    // * The process is applied for each axis on the result of the previous process.
60    // * It's uglier (using &mut) but much faster than allocating for each axis.
61    let mut data = data.to_owned();
62    let mut output = array_like(&data, data.dim(), A::zero());
63
64    for d in 0..data.ndim() {
65        maximum_filter1d_to(&data, size, Axis(d), mode, origin, &mut output);
66        if d < data.ndim() - 1 {
67            std::mem::swap(&mut output, &mut data);
68        }
69    }
70    output
71}
72
73/// Calculate a 1-D maximum filter along the given axis.
74///
75/// See `maximum_filter1d`.
76pub fn maximum_filter1d_to<S, A, D>(
77    data: &ArrayBase<S, D>,
78    size: usize,
79    axis: Axis,
80    mode: BorderMode<A>,
81    origin: isize,
82    output: &mut Array<A, D>,
83) where
84    S: Data<Elem = A>,
85    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
86    D: Dimension,
87{
88    let lower = |a, b| a <= b;
89    let higher = |a, b| a >= b;
90    min_or_max_filter(data, size, axis, mode, origin, higher, lower, output);
91}
92
93/// Calculate a 1-D minimum filter along the given axis.
94///
95/// The lines of the array along the given axis are filtered with a minimum filter of given size.
96///
97/// * `data` - The input N-D data.
98/// * `size` - Length along which to calculate 1D minimum.
99/// * `axis` - The axis of input along which to calculate.
100/// * `mode` - Method that will be used to select the padded values. See the
101///   [`CorrelateMode`](crate::CorrelateMode) enum for more information.
102/// * `origin` - Controls the placement of the filter on the input array’s pixels. A value of 0
103///   centers the filter over the pixel, with positive values shifting the filter to the left, and
104///   negative ones to the right.
105pub fn minimum_filter1d<S, A, D>(
106    data: &ArrayBase<S, D>,
107    size: usize,
108    axis: Axis,
109    mode: BorderMode<A>,
110    origin: isize,
111) -> Array<A, D>
112where
113    S: Data<Elem = A>,
114    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
115    D: Dimension,
116{
117    let mut output = data.to_owned();
118    minimum_filter1d_to(data, size, axis, mode, origin, &mut output);
119    output
120}
121
122/// Calculate a multidimensional minimum filter.
123///
124/// * `data` - The input N-D data.
125/// * `size` - Length along which to calculate 1D minimum.
126/// * `mode` - Method that will be used to select the padded values. See the
127///   [`CorrelateMode`](crate::CorrelateMode) enum for more information.
128/// * `origin` - Controls the placement of the filter on the input array’s pixels. A value of 0
129///   centers the filter over the pixel, with positive values shifting the filter to the left, and
130///   negative ones to the right.
131pub fn minimum_filter<S, A, D>(
132    data: &ArrayBase<S, D>,
133    size: usize,
134    mode: BorderMode<A>,
135    origin: isize,
136) -> Array<A, D>
137where
138    S: Data<Elem = A>,
139    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
140    D: Dimension,
141{
142    // We need 2 buffers because
143    // * We're reading neignbors so we can't read and write on the same location.
144    // * The process is applied for each axis on the result of the previous process.
145    // * It's uglier (using &mut) but much faster than allocating for each axis.
146    let mut data = data.to_owned();
147    let mut output = array_like(&data, data.dim(), A::zero());
148
149    for d in 0..data.ndim() {
150        minimum_filter1d_to(&data, size, Axis(d), mode, origin, &mut output);
151        if d < data.ndim() - 1 {
152            std::mem::swap(&mut output, &mut data);
153        }
154    }
155    output
156}
157
158/// Calculate a 1-D minimum filter along the given axis.
159///
160/// See `minimum_filter1d`.
161pub fn minimum_filter1d_to<S, A, D>(
162    data: &ArrayBase<S, D>,
163    size: usize,
164    axis: Axis,
165    mode: BorderMode<A>,
166    origin: isize,
167    output: &mut Array<A, D>,
168) where
169    S: Data<Elem = A>,
170    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
171    D: Dimension,
172{
173    let lower = |a, b| a <= b;
174    let higher = |a, b| a >= b;
175    min_or_max_filter(data, size, axis, mode, origin, lower, higher, output);
176}
177
178/// MINLIST algorithm from Richard Harter
179fn min_or_max_filter<S, A, D, F1, F2>(
180    data: &ArrayBase<S, D>,
181    filter_size: usize,
182    axis: Axis,
183    mode: BorderMode<A>,
184    origin: isize,
185    f1: F1,
186    f2: F2,
187    output: &mut Array<A, D>,
188) where
189    S: Data<Elem = A>,
190    A: Copy + Num + PartialOrd + ScalarOperand + FromPrimitive,
191    D: Dimension,
192    F1: Fn(A, A) -> bool,
193    F2: Fn(A, A) -> bool,
194{
195    if filter_size == 0 {
196        panic!("Incorrect filter size (0)");
197    }
198    if filter_size == 1 {
199        output.assign(data);
200        return;
201    }
202
203    let size1 = filter_size / 2;
204    let size2 = filter_size - size1 - 1;
205    let mode = mode.to_pad_mode();
206    let n = data.len_of(axis);
207    let pad = vec![origin_check(filter_size, origin, size1, size2)];
208    let mut buffer = Array1::from_elem(n + pad[0][0] + pad[0][1], mode.init());
209
210    #[derive(Copy, Clone, PartialEq)]
211    struct Pair<A> {
212        val: A,
213        death: usize,
214    }
215    let mut ring = VecDeque::<Pair<A>>::with_capacity(filter_size);
216
217    // The original algorihtm has been modfied to fit the `VecDeque` which makes `minpair` and
218    // `last` useless. Moreover, we need to clear the `ring` at the end because there's always
219    // at least one element left. There can be more with greater `filter_size`.
220    Zip::from(data.lanes(axis)).and(output.lanes_mut(axis)).for_each(|input, mut o| {
221        pad_to(&input, &pad, mode, &mut buffer);
222        let buffer = buffer.as_slice_memory_order().unwrap();
223
224        let mut o_idx = 0;
225        ring.push_back(Pair { val: buffer[0], death: filter_size });
226        for (&v, i) in buffer[1..].iter().zip(1..) {
227            if ring[0].death == i {
228                ring.pop_front().unwrap();
229            }
230
231            if f1(v, ring[0].val) {
232                ring[0] = Pair { val: v, death: filter_size + i };
233                while ring.len() > 1 {
234                    ring.pop_back().unwrap();
235                }
236            } else {
237                while f2(ring.back().unwrap().val, v) {
238                    ring.pop_back().unwrap();
239                }
240                ring.push_back(Pair { val: v, death: filter_size + i });
241            }
242            if i >= filter_size - 1 {
243                o[o_idx] = ring[0].val;
244                o_idx += 1;
245            }
246        }
247        ring.clear();
248    });
249}
ndarray_ndimage/filters/min_max.rs

ndarray_ndimage/filters/
min_max.rs
