1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

use ndarray::{s, ArrayBase, ArrayView3, ArrayViewMut3, Data, Ix3, Zip};

use crate::{array_like, dim_minus, Kernel3d, Mask};

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<bool>, into: ArrayViewMut3<bool>) {
        let windows = from.windows(self.dim());
        match self {
            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 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::GenericOwned(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))
            }),
            Kernel3d::GenericView(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))
            })
        }
    }

    /// Dilate the mask when at least one of the values respect the kernel:
    ///
    /// A dilation is defined as `any(w & k)`. Thus, an empty kernel will always produce an empty
    /// mask.
    #[rustfmt::skip]
    fn dilate(&self, from: ArrayView3<bool>, into: ArrayViewMut3<bool>) {
        let windows = from.windows(self.dim());
        match self {
            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
                // |(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)]
            }),
            Kernel3d::GenericOwned(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)
            }),
            Kernel3d::GenericView(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)
            })
        }
    }
}

/// Binary erosion of a 3D binary image.
///
/// * `mask` - Binary image to be eroded.
/// * `kernel` - Structuring element used for the erosion.
/// * `iterations` - The erosion is repeated iterations times.
pub fn binary_erosion<S>(mask: &ArrayBase<S, Ix3>, kernel: &Kernel3d, iterations: usize) -> Mask
where
    S: Data<Elem = bool>,
{
    // By definition, all borders are set to 0
    let (width, height, depth) = dim_minus(mask, 1);
    let mut eroded_mask = Mask::from_elem(mask.dim(), false);
    let zone = s![1..width, 1..height, 1..depth];
    eroded_mask.slice_mut(zone).assign(&mask.slice(zone));

    kernel.erode(mask.view(), eroded_mask.slice_mut(zone));

    if iterations > 1 {
        let mut previous = eroded_mask.clone();
        for it in 1..iterations {
            let (width, height, depth) = dim_minus(mask, it - 1);
            let from = previous.slice(s![it - 1..width, it - 1..height, it - 1..depth]);
            let (width, height, depth) = dim_minus(mask, it);
            let zone = s![it..width, it..height, it..depth];

            kernel.erode(from, eroded_mask.slice_mut(zone));

            if it != iterations {
                previous = eroded_mask.clone();
            }
        }
    }

    eroded_mask
}

/// Binary dilation of a 3D binary image.
///
/// * `mask` - Binary image to be dilated.
/// * `kernel` - Structuring element used for the dilation.
/// * `iterations` - The dilation is repeated iterations times.
pub fn binary_dilation<S>(mask: &ArrayBase<S, Ix3>, kernel: &Kernel3d, iterations: usize) -> Mask
where
    S: Data<Elem = bool>,
{
    let (width, height, depth) = mask.dim();
    let crop = s![1..=width, 1..=height, 1..=depth];
    let mut new_mask = array_like(mask, (width + 2, height + 2, depth + 2), false);
    new_mask.slice_mut(crop).assign(mask);

    let mut previous = new_mask.clone();
    kernel.dilate(previous.view(), new_mask.slice_mut(crop));

    for _ in 1..iterations {
        previous = new_mask.clone();
        kernel.dilate(previous.view(), new_mask.slice_mut(crop));
    }

    new_mask.slice(crop).to_owned()
}

/// Binary opening of a 3D binary image.
///
/// The opening of an input image by a structuring element is the dilation of the erosion of the
/// image by the structuring element.
///
/// * `mask` - Binary image to be opened.
/// * `kernel` - Structuring element used for the opening.
/// * `iterations` - The erosion step of the opening, then the dilation step are each repeated
///   iterations times.
pub fn binary_opening<S>(mask: &ArrayBase<S, Ix3>, kernel: &Kernel3d, iterations: usize) -> Mask
where
    S: Data<Elem = bool>,
{
    let eroded = binary_erosion(mask, kernel, iterations);
    binary_dilation(&eroded, kernel, iterations)
}

/// Binary closing of a 3D binary image.
///
/// The closing of an input image by a structuring element is the erosion of the dilation of the
/// image by the structuring element.
///
/// * `mask` - Binary image to be closed.
/// * `kernel` - Structuring element used for the closing.
/// * `iterations` - The dilation step of the closing, then the erosion step are each repeated
///   iterations times.
pub fn binary_closing<S>(mask: &ArrayBase<S, Ix3>, kernel: &Kernel3d, iterations: usize) -> Mask
where
    S: Data<Elem = bool>,
{
    let dilated = binary_dilation(mask, kernel, iterations);
    binary_erosion(&dilated, kernel, iterations)
}