diff --git a/README.md b/README.md index 5806130..e35d377 100644 --- a/README.md +++ b/README.md @@ -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 @@ -29,7 +29,7 @@ Using with Cargo ```toml [dependencies] ndarray = "0.15" -ndarray-ndimage = "0.2" +ndarray-ndimage = "0.4" ``` Contributing diff --git a/src/filters/gaussian.rs b/src/filters/gaussian.rs index e7ff755..ad65392 100644 --- a/src/filters/gaussian.rs +++ b/src/filters/gaussian.rs @@ -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(); diff --git a/src/filters/mod.rs b/src/filters/mod.rs index bc97e52..a853fea 100644 --- a/src/filters/mod.rs +++ b/src/filters/mod.rs @@ -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. diff --git a/src/filters/uniform.rs b/src/filters/uniform.rs new file mode 100644 index 0000000..7879262 --- /dev/null +++ b/src/filters/uniform.rs @@ -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( + data: &ArrayBase, + size: usize, + mode: BorderMode, +) -> Array +where + S: Data, + 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( + data: &ArrayBase, + size: usize, + axis: Axis, + mode: BorderMode, +) -> Array +where + S: Data, + 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(size: usize) -> Vec +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::(0); + } +} diff --git a/src/lib.rs b/src/lib.rs index d61c752..76cf1f8 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -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}; diff --git a/tests/filters.rs b/tests/filters.rs index 2320c47..771cb42 100644 --- a/tests/filters.rs +++ b/tests/filters.rs @@ -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) @@ -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 = (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 = (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 = (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 = (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 = (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 = (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]);