diff --git a/src/lib.rs b/src/lib.rs index f5f8209..7b0d4f1 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -3,7 +3,7 @@ //! The `ndarray-image` crate provides multidimensional image processing for `ArrayBase`, //! the *n*-dimensional array data structure provided by [`ndarray`]. -use ndarray::{Array, Array3, ArrayBase, Data, Dimension, Ix3, ShapeBuilder}; +use ndarray::{arr3, Array, Array3, ArrayBase, ArrayView3, Data, Dimension, Ix3, ShapeBuilder}; mod filters; mod interpolation; @@ -30,12 +30,63 @@ pub use pad::{pad, pad_to, PadMode}; pub type Mask = Array3; /// 3D common kernels. Also called Structuring Element. -#[derive(Clone, Copy, Debug, PartialEq)] -pub enum Kernel3d { +#[derive(Clone, Copy, PartialEq)] +pub enum Kernel3d<'a> { /// Diamond/star kernel (center and sides). + /// + /// Equivalent to `generate_binary_structure(3, 1)`. Star, - /// 3x3 cube. + /// Ball kernel (center and sides). + /// + /// Equivalent to `generate_binary_structure(3, 2)`. + Ball, + /// 3x3x3 cube. + /// + /// Equivalent to `generate_binary_structure(3, 3)`. Full, + /// Generic kernel of any 3D size. + /// + /// The generic kernels are incredibly slower on all morphological operations. + Generic(ArrayView3<'a, bool>), +} + +impl<'a> std::fmt::Debug for Kernel3d<'a> { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match *self { + Kernel3d::Star => write!(f, "Star {:?}", self.dim()), + Kernel3d::Ball => write!(f, "Ball {:?}", self.dim()), + Kernel3d::Full => write!(f, "Full {:?}", self.dim()), + Kernel3d::Generic(k) => write!(f, "Generic {:?}", k.dim()), + } + } +} + +impl<'a> Kernel3d<'a> { + /// Return the kernel dimension. + pub fn dim(&self) -> (usize, usize, usize) { + match *self { + Kernel3d::Star | Kernel3d::Ball | Kernel3d::Full => (3, 3, 3), + Kernel3d::Generic(k) => k.dim(), + } + } + + /// Return the actual 3D array. + pub fn array(&self) -> Array3 { + match *self { + Kernel3d::Star => arr3(&[ + [[false, false, false], [false, true, false], [false, false, false]], + [[false, true, false], [true, true, true], [false, true, false]], + [[false, false, false], [false, true, false], [false, false, false]], + ]), + Kernel3d::Ball => arr3(&[ + [[false, true, false], [true, true, true], [false, true, false]], + [[true, true, true], [true, true, true], [true, true, true]], + [[false, true, false], [true, true, true], [false, true, false]], + ]), + Kernel3d::Full => Array3::from_elem((3, 3, 3), true), + Kernel3d::Generic(k) => k.to_owned(), + } + } } /// Utilitary function that returns a new *n*-dimensional array of dimension `shape` with the same diff --git a/src/morphology.rs b/src/morphology.rs index 33db7b2..423c403 100644 --- a/src/morphology.rs +++ b/src/morphology.rs @@ -2,42 +2,54 @@ use ndarray::{s, ArrayBase, ArrayView3, ArrayViewMut3, Data, Ix3, Zip}; use crate::{array_like, dim_minus, Kernel3d, Mask}; -impl Kernel3d { +impl<'a> Kernel3d<'a> { /// Erode the mask when at least one of the values doesn't respect the kernel. /// /// An erosion is defined either as `all(!(!w & k))` or `!any(!w & k)`. Thus, an empty kernel /// will always produce a full mask. #[rustfmt::skip] fn erode(&self, from: ArrayView3, into: ArrayViewMut3) { + let windows = from.windows(self.dim()); match *self { - Kernel3d::Full => { - // TODO This is incredibly ugly but this is much faster (3x) than the Zip - //Zip::from(from.windows((3, 3, 3))) - // .map_assign_into(into, |w| Zip::from(w).all(|&m| m)); - Zip::from(from.windows((3, 3, 3))).map_assign_into(into, |w| { - w[(0, 0, 0)] && w[(0, 0, 1)] && w[(0, 0, 2)] - && w[(0, 1, 0)] && w[(0, 1, 1)] && w[(0, 1, 2)] - && w[(0, 2, 0)] && w[(0, 2, 1)] && w[(0, 2, 2)] - && w[(1, 0, 0)] && w[(1, 0, 1)] && w[(1, 0, 2)] - && w[(1, 1, 0)] && w[(1, 1, 1)] && w[(1, 1, 2)] - && w[(1, 2, 0)] && w[(1, 2, 1)] && w[(1, 2, 2)] - && w[(2, 0, 0)] && w[(2, 0, 1)] && w[(2, 0, 2)] - && w[(2, 1, 0)] && w[(2, 1, 1)] && w[(2, 1, 2)] - && w[(2, 2, 0)] && w[(2, 2, 1)] && w[(2, 2, 2)] - }) - } - Kernel3d::Star => Zip::from(from.windows((3, 3, 3))).map_assign_into(into, |w| { + Kernel3d::Full => Zip::from(windows).map_assign_into(into, |w| { + // This is incredibly ugly but this is much faster (3x) than the Zip + // |w| Zip::from(w).all(|&m| m) + w[(0, 0, 0)] && w[(0, 0, 1)] && w[(0, 0, 2)] + && w[(0, 1, 0)] && w[(0, 1, 1)] && w[(0, 1, 2)] + && w[(0, 2, 0)] && w[(0, 2, 1)] && w[(0, 2, 2)] + && w[(1, 0, 0)] && w[(1, 0, 1)] && w[(1, 0, 2)] + && w[(1, 1, 0)] && w[(1, 1, 1)] && w[(1, 1, 2)] + && w[(1, 2, 0)] && w[(1, 2, 1)] && w[(1, 2, 2)] + && w[(2, 0, 0)] && w[(2, 0, 1)] && w[(2, 0, 2)] + && w[(2, 1, 0)] && w[(2, 1, 1)] && w[(2, 1, 2)] + && w[(2, 2, 0)] && w[(2, 2, 1)] && w[(2, 2, 2)] + }), + Kernel3d::Ball => Zip::from(windows).map_assign_into(into, |w| { + w[(0, 0, 1)] + && w[(0, 1, 0)] && w[(0, 1, 1)] && w[(0, 1, 2)] + && w[(0, 2, 1)] + && w[(1, 0, 0)] && w[(1, 0, 1)] && w[(1, 0, 2)] + && w[(1, 1, 0)] && w[(1, 1, 1)] && w[(1, 1, 2)] + && w[(1, 2, 0)] && w[(1, 2, 1)] && w[(1, 2, 2)] + && w[(2, 0, 1)] + && w[(2, 1, 0)] && w[(2, 1, 1)] && w[(2, 1, 2)] + && w[(2, 2, 1)] + }), + Kernel3d::Star => Zip::from(windows).map_assign_into(into, |w| { // This ugly condition is equivalent to // *mask = !w.iter().zip(&kernel).any(|(w, k)| !w & k) - // but it's around 5x faster because there's no branch misprediction - w[(0, 1, 1)] - && w[(1, 0, 1)] - && w[(1, 1, 0)] - && w[(1, 1, 1)] - && w[(1, 1, 2)] - && w[(1, 2, 1)] - && w[(2, 1, 1)] + // but it's extremely faster because there's no branch misprediction + w[(0, 1, 1)] + && w[(1, 0, 1)] + && w[(1, 1, 0)] && w[(1, 1, 1)] && w[(1, 1, 2)] + && w[(1, 2, 1)] + && w[(2, 1, 1)] }), + Kernel3d::Generic(kernel) => Zip::from(windows).map_assign_into(into, |w| { + // TODO Maybe use Zip::any when available + // !Zip::from(w).and(kernel).any(|(&w, &k)| !w & k) + Zip::from(w).and(kernel).all(|&w, &k| !(!w & k)) + }) } } @@ -47,34 +59,47 @@ impl Kernel3d { /// mask. #[rustfmt::skip] fn dilate(&self, from: ArrayView3, into: ArrayViewMut3) { + let windows = from.windows(self.dim()); match *self { - Kernel3d::Full => { - // TODO This is incredibly ugly but this is much faster (6x) than the Zip - //Zip::from(from.windows((3, 3, 3))).map_assign_into(into, |w| w.iter().any(|&w| w)) - Zip::from(from.windows((3, 3, 3))).map_assign_into(into, |w| { - w[(0, 0, 0)] || w[(0, 0, 1)] || w[(0, 0, 2)] - || w[(0, 1, 0)] || w[(0, 1, 1)] || w[(0, 1, 2)] - || w[(0, 2, 0)] || w[(0, 2, 1)] || w[(0, 2, 2)] - || w[(1, 0, 0)] || w[(1, 0, 1)] || w[(1, 0, 2)] - || w[(1, 1, 0)] || w[(1, 1, 1)] || w[(1, 1, 2)] - || w[(1, 2, 0)] || w[(1, 2, 1)] || w[(1, 2, 2)] - || w[(2, 0, 0)] || w[(2, 0, 1)] || w[(2, 0, 2)] - || w[(2, 1, 0)] || w[(2, 1, 1)] || w[(2, 1, 2)] - || w[(2, 2, 0)] || w[(2, 2, 1)] || w[(2, 2, 2)] - }) - } - Kernel3d::Star => Zip::from(from.windows((3, 3, 3))).map_assign_into(into, |w| { + Kernel3d::Full => Zip::from(windows).map_assign_into(into, |w| { + // This is incredibly ugly but this is much faster (6x) than the Zip + // |w| w.iter().any(|&w| w)) + w[(0, 0, 0)] || w[(0, 0, 1)] || w[(0, 0, 2)] + || w[(0, 1, 0)] || w[(0, 1, 1)] || w[(0, 1, 2)] + || w[(0, 2, 0)] || w[(0, 2, 1)] || w[(0, 2, 2)] + || w[(1, 0, 0)] || w[(1, 0, 1)] || w[(1, 0, 2)] + || w[(1, 1, 0)] || w[(1, 1, 1)] || w[(1, 1, 2)] + || w[(1, 2, 0)] || w[(1, 2, 1)] || w[(1, 2, 2)] + || w[(2, 0, 0)] || w[(2, 0, 1)] || w[(2, 0, 2)] + || w[(2, 1, 0)] || w[(2, 1, 1)] || w[(2, 1, 2)] + || w[(2, 2, 0)] || w[(2, 2, 1)] || w[(2, 2, 2)] + }), + Kernel3d::Ball => Zip::from(windows).map_assign_into(into, |w| { + w[(0, 0, 1)] + || w[(0, 1, 0)] || w[(0, 1, 1)] || w[(0, 1, 2)] + || w[(0, 2, 1)] + || w[(1, 0, 0)] || w[(1, 0, 1)] || w[(1, 0, 2)] + || w[(1, 1, 0)] || w[(1, 1, 1)] || w[(1, 1, 2)] + || w[(1, 2, 0)] || w[(1, 2, 1)] || w[(1, 2, 2)] + || w[(2, 0, 1)] + || w[(2, 1, 0)] || w[(2, 1, 1)] || w[(2, 1, 2)] + || w[(2, 2, 1)] + }), + Kernel3d::Star => Zip::from(windows).map_assign_into(into, |w| { // This ugly condition is equivalent to - // *mask = w.iter().zip(&kernel).any(|(w, k)| w & k) + // |(w, k)| w & k // but it's around 5x faster because there's no branch misprediction - w[(0, 1, 1)] - || w[(1, 0, 1)] - || w[(1, 1, 0)] - || w[(1, 1, 1)] - || w[(1, 1, 2)] - || w[(1, 2, 1)] - || w[(2, 1, 1)] + w[(0, 1, 1)] + || w[(1, 0, 1)] + || w[(1, 1, 0)] || w[(1, 1, 1)] || w[(1, 1, 2)] + || w[(1, 2, 1)] + || w[(2, 1, 1)] }), + Kernel3d::Generic(kernel) => Zip::from(windows).map_assign_into(into, |w| { + // TODO Use Zip::any when available + // Zip::from(w).and(kernel).any(|idx, &w, &k| w & k) + w.iter().zip(&kernel).any(|(w, k)| w & k) + }) } } } diff --git a/tests/morphology.rs b/tests/morphology.rs index d63f2af..74a5f02 100644 --- a/tests/morphology.rs +++ b/tests/morphology.rs @@ -49,10 +49,27 @@ fn test_binary_erosion_hole() { gt.slice_mut(s![5, 5, 4..7]).fill(false); assert_eq!(gt, binary_erosion(&mask, Kernel3d::Star, 1)); + assert_eq!(gt, binary_erosion(&mask, Kernel3d::Generic(Kernel3d::Star.array().view()), 1)); } #[test] // Results verified with the `binary_erosion` function from SciPy. (v1.7.0) -fn test_binary_erosion_cube_kernel() { +fn test_binary_erosion_ball_kernel() { + let mut mask = Mask::from_elem((11, 11, 11), true); + mask[(5, 5, 5)] = false; + + let mut gt = Mask::from_elem((11, 11, 11), false); + let (width, height, depth) = dim_minus(&mask, 1); + gt.slice_mut(s![1..width, 1..height, 1..depth]).fill(true); + // Remove the ball shape in the image center. + gt.slice_mut(s![4..7, 4..7, 4..7]).fill(false); + gt.slice_mut(s![4..7; 2, 4..7; 2, 4..7; 2]).fill(true); + + assert_eq!(gt, binary_erosion(&mask, Kernel3d::Ball, 1)); + assert_eq!(gt, binary_erosion(&mask, Kernel3d::Generic(Kernel3d::Ball.array().view()), 1)); +} + +#[test] // Results verified with the `binary_erosion` function from SciPy. (v1.7.0) +fn test_binary_erosion_full_kernel() { let mut mask = Mask::from_elem((11, 11, 11), true); mask[(5, 5, 5)] = false; @@ -63,6 +80,7 @@ fn test_binary_erosion_cube_kernel() { gt.slice_mut(s![4..7, 4..7, 4..7]).fill(false); assert_eq!(gt, binary_erosion(&mask, Kernel3d::Full, 1)); + assert_eq!(gt, binary_erosion(&mask, Kernel3d::Generic(Kernel3d::Full.array().view()), 1)); } #[test] // Results verified with the `binary_dilation` function from SciPy. (v1.7.0) @@ -98,18 +116,29 @@ fn test_binary_dilation_plain() { gt.slice_mut(s![2..w - 1, h - 1, 2..d - 1]).fill(true); assert_eq!(gt, binary_dilation(&mask.view(), Kernel3d::Star, 1)); + assert_eq!(gt, binary_dilation(&mask, Kernel3d::Generic(Kernel3d::Star.array().view()), 1)); + + let mut mask = Mask::from_elem((w, h, d), false); + mask.slice_mut(s![4, 4, 4..]).fill(true); + let mut gt = Mask::from_elem((w, h, d), false); + gt.slice_mut(s![3..6, 3..6, 3..]).fill(true); + gt.slice_mut(s![3..6; 2, 3..6; 2, 3]).fill(false); + assert_eq!(gt, binary_dilation(&mask.view(), Kernel3d::Ball, 1)); + assert_eq!(gt, binary_dilation(&mask, Kernel3d::Generic(Kernel3d::Ball.array().view()), 1)); let mut mask = Mask::from_elem((w, h, d), false); mask[(4, 4, 4)] = true; let mut gt = Mask::from_elem((w, h, d), false); gt.slice_mut(s![2.., 2.., 2..]).fill(true); assert_eq!(gt, binary_dilation(&mask.view(), Kernel3d::Full, 2)); + assert_eq!(gt, binary_dilation(&mask, Kernel3d::Generic(Kernel3d::Full.array().view()), 2)); let mut mask = Mask::from_elem((w, h, d), false); mask[(4, 5, 5)] = true; let mut gt = Mask::from_elem((w, h, d), false); gt.slice_mut(s![1.., 2.., 2..]).fill(true); assert_eq!(gt, binary_dilation(&mask.view(), Kernel3d::Full, 3)); + assert_eq!(gt, binary_dilation(&mask, Kernel3d::Generic(Kernel3d::Full.array().view()), 3)); } #[test] // Results verified with the `binary_dilation` function from SciPy. (v1.7.0)