Skip to content

Commit

Permalink
Add uniform_filter/1d
Browse files Browse the repository at this point in the history
  • Loading branch information
matsjoyce-refeyn committed Sep 13, 2023
1 parent 1bc8743 commit f86aa12
Show file tree
Hide file tree
Showing 6 changed files with 279 additions and 6 deletions.
6 changes: 3 additions & 3 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@ This crate provides multidimensional image processing for [`ndarray`]'s `ArrayBa
It aims to:
- be a Rust replacement for [`scipy.ndimage`] with some other tools like [`numpy.pad`] and anything else relevant to image processing. We do not want all options and arguments offered by `scipy.ndimage` because some of them are incompatible with Rust. We hope to offer the most used ones.
- be faster or as fast as `scipy.ndimage`. Most of it is cythonized so it's not as easy as it seems. In fact, I'm usually unable to be faster than SciPy but it does happen on some functions.
- avoid using `unsafe`. This is not an unbreakable rule. Its usage will be evaluated and dicussed in the pull requests.
- avoid using `unsafe`. This is not an unbreakable rule. Its usage will be evaluated and discussed in the pull requests.

Currently available routines include:
- Filters: convolve/1d, correlate/1d, gaussian_filter/1d, min/max_filter/1d, median_filter, prewitt, sobel
- Filters: convolve/1d, correlate/1d, gaussian_filter/1d, min/max_filter/1d, uniform_filter/1d, median_filter, prewitt, sobel
- Fourier filters: none. Please use the excellent [`rustfft`] crate
- Interpolation: shift, spline_filter/1d, zoom
- Measurements: label, label_histogram, largest_connected_components, most_frequent_label
Expand All @@ -29,7 +29,7 @@ Using with Cargo
```toml
[dependencies]
ndarray = "0.15"
ndarray-ndimage = "0.2"
ndarray-ndimage = "0.4"
```

Contributing
Expand Down
2 changes: 1 addition & 1 deletion src/filters/gaussian.rs
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ where
let half = weights.len() / 2;

// We need 2 buffers because
// * We're reading neignbors so we can't read and write on the same location.
// * We're reading neighbours so we can't read and write on the same location.
// * The process is applied for each axis on the result of the previous process.
// * It's uglier (using &mut) but much faster than allocating for each axis.
let mut data = data.to_owned();
Expand Down
1 change: 1 addition & 0 deletions src/filters/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ pub mod gaussian;
pub mod median;
pub mod min_max;
pub mod symmetry;
pub mod uniform;

// TODO We might want to offer all NumPy mode (use PadMode instead)
/// Method that will be used to determines how the input array is extended beyond its boundaries.
Expand Down
110 changes: 110 additions & 0 deletions src/filters/uniform.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
use ndarray::{Array, ArrayBase, Axis, Data, Dimension};
use num_traits::{Float, FromPrimitive};

use crate::{array_like, BorderMode};

use super::{con_corr::inner_correlate1d, symmetry::SymmetryStateCheck};

/// Uniform filter for n-dimensional arrays.
///
/// Currently hardcoded with the `PadMode::Reflect` padding mode.
///
/// * `data` - The input N-D data.
/// * `size` - The len
/// * `mode` - Method that will be used to select the padded values. See the
/// [`BorderMode`](crate::BorderMode) enum for more information.
///
/// **Panics** if `size` is zero, or one of the axis' lengths is lower than `size`.
pub fn uniform_filter<S, A, D>(
data: &ArrayBase<S, D>,
size: usize,
mode: BorderMode<A>,
) -> Array<A, D>
where
S: Data<Elem = A>,
A: Float + FromPrimitive + 'static,
for<'a> &'a [A]: SymmetryStateCheck,
D: Dimension,
{
let weights = weights(size);
let half = weights.len() / 2;

// We need 2 buffers because
// * We're reading neighbours so we can't read and write on the same location.
// * The process is applied for each axis on the result of the previous process.
// * It's uglier (using &mut) but much faster than allocating for each axis.
let mut data = data.to_owned();
let mut output = array_like(&data, data.dim(), A::zero());

for d in 0..data.ndim() {
// TODO This can be made to work if the padding modes (`reflect`, `symmetric`, `wrap`) are
// more robust. One just needs to reflect the input data several times if the `weights`
// length is greater than the input array. It works in SciPy because they are looping on a
// size variable instead of running the algo only once like we do.
let n = data.len_of(Axis(d));
if half > n {
panic!("Data size is too small for the inputs (sigma and truncate)");
}

inner_correlate1d(&data, &weights, Axis(d), mode, 0, &mut output);
if d != data.ndim() - 1 {
std::mem::swap(&mut output, &mut data);
}
}
output
}

/// Uniform filter for 1-dimensional arrays.
///
/// * `data` - The input N-D data.
/// * `size` - Length of the uniform filter.
/// * `axis` - The axis of input along which to calculate.
/// * `mode` - Method that will be used to select the padded values. See the
/// [`BorderMode`](crate::BorderMode) enum for more information.
///
/// **Panics** if `size` is zero, or the axis length is lower than `size`.
pub fn uniform_filter1d<S, A, D>(
data: &ArrayBase<S, D>,
size: usize,
axis: Axis,
mode: BorderMode<A>,
) -> Array<A, D>
where
S: Data<Elem = A>,
A: Float + FromPrimitive + 'static,
for<'a> &'a [A]: SymmetryStateCheck,
D: Dimension,
{
let weights = weights(size);
let mut output = array_like(&data, data.dim(), A::zero());
inner_correlate1d(data, &weights, axis, mode, 0, &mut output);
output
}

fn weights<A>(size: usize) -> Vec<A>
where
A: Float + FromPrimitive + 'static,
{
if size == 0 {
panic!("Size is zero");
}
[A::one() / A::from_usize(size).unwrap()].repeat(size)
}

#[cfg(test)]
mod tests {
use super::*;
use approx::assert_relative_eq;

#[test]
fn test_weights() {
assert_relative_eq!(weights(5).as_slice(), &[0.2, 0.2, 0.2, 0.2, 0.2][..], epsilon = 1e-7);
assert_relative_eq!(weights(1).as_slice(), &[1.0][..], epsilon = 1e-7);
}

#[should_panic]
#[test]
fn test_weights_zero() {
weights::<f64>(0);
}
}
1 change: 1 addition & 0 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ pub use filters::{
maximum_filter, maximum_filter1d, maximum_filter1d_to, minimum_filter, minimum_filter1d,
minimum_filter1d_to,
},
uniform::{uniform_filter, uniform_filter1d},
BorderMode,
};
pub use interpolation::{shift, spline_filter, spline_filter1d, zoom};
Expand Down
165 changes: 163 additions & 2 deletions tests/filters.rs
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@ use ndarray::{arr1, arr2, s, Array1, Array2, Axis};

use ndarray_ndimage::{
convolve, convolve1d, correlate, correlate1d, gaussian_filter, maximum_filter,
maximum_filter1d, median_filter, minimum_filter, minimum_filter1d, prewitt, sobel, BorderMode,
Mask,
maximum_filter1d, median_filter, minimum_filter, minimum_filter1d, prewitt, sobel,
uniform_filter, BorderMode, Mask,
};

#[test] // Results verified with SciPy. (v1.9.0)
Expand Down Expand Up @@ -670,6 +670,167 @@ fn test_gaussian_filter_panic() {
let _ = gaussian_filter(&a, 2.0, 0, BorderMode::Reflect, 4);
}

#[test] // Results verified with SciPy. (v1.9.1)
fn test_uniform_filter_1d() {
let a: Array1<f32> = (0..7).map(|v| v as f32).collect();
assert_relative_eq!(
uniform_filter(&a, 1, BorderMode::Reflect),
arr1(&[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]),
epsilon = 1e-7
);
assert_relative_eq!(
uniform_filter(&a.view(), 2, BorderMode::Reflect),
arr1(&[0., 0.5, 1.5, 2.5, 3.5, 4.5, 5.5]),
epsilon = 1e-7
);
assert_relative_eq!(
uniform_filter(&a.view(), 3, BorderMode::Reflect),
arr1(&[0.33333333, 1., 2., 3., 4., 5., 5.66666667]),
epsilon = 1e-7
);
}

#[test] // Results verified with SciPy. (v1.9.1)
fn test_uniform_filter_2d() {
let a: Array1<f32> = (0..70).step_by(2).map(|v| v as f32).collect();
let mut a = a.into_shape((5, 7)).unwrap();
a[(0, 0)] = 17.0;
assert_relative_eq!(
uniform_filter(&a, 4, BorderMode::Reflect),
arr2(&[
[12.25, 12.75, 12.125, 12., 14., 16., 17.5],
[15.75, 16.25, 15.625, 15.5, 17.5, 19.5, 21.],
[24.125, 24.625, 25.0625, 26., 28., 30., 31.5],
[36., 36.5, 38., 40., 42., 44., 45.5],
[46.5, 47., 48.5, 50.5, 52.5, 54.5, 56.]
]),
epsilon = 1e-5
);
let a: Array1<f32> = (0..84).step_by(2).map(|v| v as f32).collect();
let mut a = a.into_shape((6, 7)).unwrap();
a[(0, 0)] = 8.5;
assert_relative_eq!(
uniform_filter(&a, 5, BorderMode::Reflect),
arr2(&[
[14.16, 14.96, 15.88, 17.2, 19.2, 20.8, 21.6],
[19.76, 20.56, 21.48, 22.8, 24.8, 26.4, 27.2],
[30.28, 31.08, 32.34, 34., 36., 37.6, 38.4],
[43.6, 44.4, 46., 48., 50., 51.6, 52.4],
[54.8, 55.6, 57.2, 59.2, 61.2, 62.8, 63.6],
[60.4, 61.2, 62.8, 64.8, 66.8, 68.4, 69.2]
]),
epsilon = 1e-5
);

let a: Array1<f32> = (0..112).step_by(2).map(|v| v as f32).collect();
let mut a = a.into_shape((8, 7)).unwrap();
a[(0, 0)] = 18.2;
assert_relative_eq!(
uniform_filter(&a, 3, BorderMode::Reflect),
arr2(&[
[13.42222222, 10.71111111, 8.66666667, 10.66666667, 12.66666667, 14.66666667, 16.],
[18.71111111, 18.02222222, 18., 20., 22., 24., 25.33333333],
[28.66666667, 30., 32., 34., 36., 38., 39.33333333],
[42.66666667, 44., 46., 48., 50., 52., 53.33333333],
[56.66666667, 58., 60., 62., 64., 66., 67.33333333],
[70.66666667, 72., 74., 76., 78., 80., 81.33333333],
[84.66666667, 86., 88., 90., 92., 94., 95.33333333],
[94., 95.33333333, 97.33333333, 99.33333333, 101.33333333, 103.33333333, 104.66666667]
]),
epsilon = 1e-4
);
}

#[test] // Results verified with SciPy. (v1.9.1)
fn test_uniform_filter_3d() {
let a: Array1<f32> = (0..720).map(|v| v as f32 / 50.0).collect();
let mut a = a.into_shape((10, 9, 8)).unwrap();
a[(0, 0, 0)] = 0.2;
a[(3, 3, 3)] = 1.0;

let g = uniform_filter(&a, 6, BorderMode::Reflect);
assert_relative_eq!(
g.slice(s![0, .., ..]),
arr2(&[
[1.62740741, 1.63074074, 1.64074074, 1.6537037, 1.67, 1.69, 1.70666667, 1.71666667],
[
1.65407407, 1.65740741, 1.66740741, 1.68037037, 1.69666667, 1.71666667, 1.73333333,
1.74333333
],
[
1.73407407, 1.73740741, 1.74740741, 1.76037037, 1.77666667, 1.79666667, 1.81333333,
1.82333333
],
[1.8637037, 1.86703704, 1.87703704, 1.89185185, 1.91, 1.93, 1.94666667, 1.95666667],
[2.02, 2.02333333, 2.03333333, 2.05, 2.07, 2.09, 2.10666667, 2.11666667],
[2.18, 2.18333333, 2.19333333, 2.21, 2.23, 2.25, 2.26666667, 2.27666667],
[2.34, 2.34333333, 2.35333333, 2.37, 2.39, 2.41, 2.42666667, 2.43666667],
[2.47333333, 2.47666667, 2.48666667, 2.50333333, 2.52333333, 2.54333333, 2.56, 2.57],
[2.55333333, 2.55666667, 2.56666667, 2.58333333, 2.60333333, 2.62333333, 2.64, 2.65]
]),
epsilon = 1e-4
);
assert_relative_eq!(
g.slice(s![9, .., ..]),
arr2(&[
[11.46, 11.46333333, 11.47333333, 11.49, 11.51, 11.53, 11.54666667, 11.55666667],
[
11.48666667,
11.49,
11.5,
11.51666667,
11.53666667,
11.55666667,
11.57333333,
11.58333333
],
[
11.56666667,
11.57,
11.58,
11.59666667,
11.61666667,
11.63666667,
11.65333333,
11.66333333
],
[11.7, 11.70333333, 11.71333333, 11.73, 11.75, 11.77, 11.78666667, 11.79666667],
[11.86, 11.86333333, 11.87333333, 11.89, 11.91, 11.93, 11.94666667, 11.95666667],
[12.02, 12.02333333, 12.03333333, 12.05, 12.07, 12.09, 12.10666667, 12.11666667],
[12.18, 12.18333333, 12.19333333, 12.21, 12.23, 12.25, 12.26666667, 12.27666667],
[
12.31333333,
12.31666667,
12.32666667,
12.34333333,
12.36333333,
12.38333333,
12.4,
12.41
],
[
12.39333333,
12.39666667,
12.40666667,
12.42333333,
12.44333333,
12.46333333,
12.48,
12.49
]
]),
epsilon = 1e-4
);
}

#[should_panic]
#[test] // Results verified with SciPy. (v1.9.1)
fn test_uniform_filter_panic() {
let a: Array1<f32> = (0..7).map(|v| v as f32).collect();

let _ = uniform_filter(&a, 0, BorderMode::Reflect);
}

#[test] // Results verified with SciPy. (v1.9.0)
fn test_prewitt() {
let a = arr1(&[2.0, 8.1, 0.5, 4.0, 1.1, 9.0, 9.0, 0.8]);
Expand Down

0 comments on commit f86aa12

Please sign in to comment.