ndarray_ndimage/filters/
min_max.rs

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

ndarray_ndimage/filters/
min_max.rs
