Keywords

1 Introduction

Morphological operations of dilation and erosion are well-established image processing methods, and significant effort has been put into developing efficient algorithms for gray-scale morphology. An important advance, also for our work, is the van Herk/Gil-Werman algorithm [8, 20] which computes 1D dilation and erosion in constant time per pixel – regardless of the length of the structuring element. This opens a possibility for efficient 2D and 3D algorithms, if the structuring element can be decomposed in linear components.

Decompositions and approximations of various 2D structuring elements have been known for decades, while 3D structuring elements have received almost no attention. This also holds true for the 3D sphere, which is a structuring element of large practical importance, as it is invariant to rotation. To our knowledge, approximation of the 3D sphere has not been treated in the context relevant for morphological filtering of digital volumes. Furthermore, when investigating how 3D morphology is handled by software commonly used for image processing, we found no evidence of fully utilizing the approximation of the sphere.

Efficient computation is especially important in 3D image processing, where the volume of a spherical structuring element grows cubically with the radius, and datasets are often large. In particular, applications requiring multiple passes of morphological filtering, like granulometry or porosity analysis, might be seriously affected by the lack of efficient algorithms for 3D dilation and erosion.

In this paper we treat the problem of approximating a 3D sphere of an arbitrary radius with linear components. Our main contribution is an algorithm for approximating a 3D sphere using a zonohedron. The resulting linear components can be used for morphological filtering in any image processing framework. When combined with the van Herk/Gil-Werman algorithm, computational complexity of dilation or erosion is independent of the radius of the sphere. Furthermore, based on our sphere approximation, we implemented an efficient GPU-based dilation/erosion which is available at https://github.com/patmjen/Zonohedra. To exemplify the use of our method we apply gray-scale opening to detect spheres in \(\mu \)CT data of cement and stone wool. And lastly, we perform tests showing the benefits of our sphere approximation, and the efficiency of our implementation.

1.1 Related Work

Research in increasing the computational speed in 2D morphology [18] generally falls in two categories. The first seek to exploit simplifying assumptions, for example by only considering binary or integer data [13, 19] or by assuming a flat structuring element [17]. The other category attempts to decompose or approximate the structuring element with a series of shapes which can be processed quickly. The flat line segment is particularly attractive thanks to the van Herk/Gil-Werman [8, 20] algorithm. They have shown that erosion/dilation with a horizontal, vertical or diagonal structuring element requires only three operations per pixel. This was improved [7] and generalized to operations with discrete lines at arbitrary angles [16].

Many structuring elements cannot be decomposed into line segments and one has to resort to a approximation. As stated by Adams [1], a 2D disk can be approximated via successive dilations/erosions of line segments. This technique was extended to 3D in the context of solid modeling [10]. However, they only deal with the case where all line segments have equal and continuous length, which makes it impractical for digital images defined on a discrete grid. An extensive theoretical study of the approximation of n-dimensional spheres (and other shapes) using dilations of continuous line segments has been performed by Bourgain et al. [3] and Campi et al. [5]. Techniques for exact decomposition have also been explored, e.g. where a structuring element is decomposed into elements of size \(3 \times 3\) [14, 15], or with an explicit focus on spheres [21]. However, these still scale with the size of the structuring element and are therefore still problematic when working with large structures – especially for volumetric data.

2 Method

2.1 Computational Complexity of Mathematical Morphology

Fundamental operations of mathematical morphology are the dilation and erosion. For an image I and a binary (flat) structuring element S, dilation and erosion are given by, respectively,

$$\begin{aligned} (I \oplus S)(\mathbf x ) = \max _\mathbf{y \in S} \left\{ I(\mathbf x + \mathbf y ) \right\} , \end{aligned}$$
(1)

and

$$\begin{aligned} (I \ominus S)(\mathbf x ) = \min _\mathbf{y \in S} \left\{ I(\mathbf x + \mathbf y ) \right\} , \end{aligned}$$
(2)

where \( \mathbf y \in S \) refers to the position of non-zero values in S. A direct implementation of dilation/erosion leads to an algorithm requiring s comparisons per image pixel, s being a number of non-zero values in S. In 3D morphology using a spherical structuring element, the direct algorithm scales with the cube of the sphere radius.

An important reduction in computational complexity can be made in the special case where non-zero voxels are arranged on a line. In this case, the van Herk/Gil-Werman algorithm [8, 20] allows dilation/erosion with only three comparisons per image pixel, independent of the structuring element size.

Since dilation is associative [9] we have

$$\begin{aligned} I \oplus (S_1 \oplus S_2 \oplus ... \oplus S_m) = (...((I \oplus S_1) \oplus S_2) \oplus ....) \oplus S_m, \end{aligned}$$
(3)

with the dual property valid for erosion

$$\begin{aligned} I \ominus (S_1 \oplus S_2 \oplus ... \oplus S_m) = (...((I \ominus S_1) \ominus S_2) \ominus ....) \ominus S_m. \end{aligned}$$
(4)

Thus, if a sphere S can be decomposed such that \( S = S_1 \oplus S_2 \oplus ... \oplus S_m \), dilation/erosion with S may be accomplished via successive dilations/erosions with \( S_1, S_2, ..., S_m \). When the number of non-zero values in S exceeds the sum of non-zero values in the decomposition, the implementation utilizing decomposition of S will be more efficient. Even better, if m is fixed and decomposition is such that dilations/erosions with \( S_1, S_2,...,S_m \) exploit the constant-time algorithm of van Herk/Gil-Werman, we can achieve an implementation with computational complexity independent of the radius of S. This will reduce the computational work per voxel from \(O(r^3) \) to O(1) and allow efficient morphological filtering with large structuring elements.

Our aim is to approximate a 3D spherical structuring element using a series of dilations with line segments. The similar 2D approximation, utilized for 2D image morphology by Adams [1], is illustrated in Fig. 1. Since a series of dilations with line segments corresponds to a Minkowski sum, the resulting shape will, by definition, be a zonohedron [11].

Fig. 1.
figure 1

Successive dilations with line segments. The resulting shape approximates a disk.

Full size image

2.2 Constructing the Zonohedral Approximation

A zonohedron is a convex polyhedron defined as a Minkowski sum (dilation) of a finite number of line segments [11]. For a zonohedron centered around the origin, line segments are usually represented as vectors \( \mathcal {V} = \{ \mathbf v _1, \mathbf v _2,..., \mathbf v _n \}\), \(\mathbf v _i \in \mathbb {R}^3\), called generators. The zonohedron is then given as

$$\begin{aligned} Z = \left\{ \sum _{i=1}^n d_i \mathbf v _i \,\mathrel {\Big |}\, -1/2 \le d_i \le 1/2, \mathbf v _i \in \mathcal {V} \right\} . \end{aligned}$$
(5)

For the purpose of volumetric filtering, generators with a simple voxelization are of highest interest. Therefore, we restrict ourselves to three sets of direction vectors

$$\begin{aligned} \begin{aligned} \mathcal {V}_1&= \left\{ (1,0,0), (0,1,0), (0,0,1) \right\} ,\\ \mathcal {V}_2&= \left\{ (1,1,0), (-1,1,0), (0,1,1), (0,-1,1),(1,0,1),(1,0,-1) \right\} ,\\ \mathcal {V}_3&= \left\{ (1,1,1), (-1,1,1), (1,-1,1), (1,1,-1) \right\} . \end{aligned} \end{aligned}$$
(6)

On a unit cube, vectors in \(\mathcal {V}_1\) connect centers of opposite faces, those in \(\mathcal {V}_2\) connect midpoints of opposite edges, and those in \(\mathcal {V}_3\) connect opposite corners.

To change the shape and the size of a zonohedron, while keeping it roughly spherical, we scale the direction vectors in each set with \( a_j \ge 0 \,, j = 1,2,3\), and consider zonohedra defined by the set of 13 generators

$$\begin{aligned} \mathcal {V} = a_1\mathcal {V}_1 \,\cup \, a_2 \mathcal {V}_2\,\cup \, a_3 \mathcal {V}_3\,, \end{aligned}$$
(7)

where \(a_j \mathcal {V}_j = \{ a_j \mathbf v \mid \mathbf v \in \mathcal {V}_j \}\). Examples of zonohedra generated by these vectors are shown in Fig. 2, which also shows the most general case, a truncated rhombicuboctahedron. Note that for integer values of \(a_1, a_2, a_3\) the generators will all have integer coordinates, and simple voxelization.

Fig. 2.
figure 2

A few zonohedra constructed from generator set \( \mathcal {V} = a_1 \mathcal {V}_1 \,\cup \, a_2 \mathcal {V}_2 \,\cup \, a_3 \mathcal {V}_3 \). (a) Truncated octahedron from \( \mathcal {V}_2 \). (b) Truncated cuboctahedron from \( \mathcal {V}_1 \cup \mathcal {V}_2 \). (c) Truncated cuboctahedron from \( 3\mathcal {V}_1 \cup \mathcal {V}_2 \). (d) Truncated rhombicuboctahedron from \( \mathcal {V}_1 \cup \mathcal {V}_2 \cup \mathcal {V}_3 \). (e) Truncated rhombicuboctahedron from \( 4\mathcal {V}_1 \cup \mathcal {V}_2 \cup 2\mathcal {V}_3 \).

Full size image

2.3 Minimizing Approximation Error

To find a good approximation, we chose to minimize the directed Hausdorff distance \(d_\text {h} \) between the surface of Z, denoted \(\partial Z\), and a sphere with radius r, denoted \(\partial B\). The directed Hausdorff distance is defined by [12]

$$\begin{aligned} d_\text {h}(\partial Z, \partial B) = \sup _\mathbf{x \in \partial Z} \inf _\mathbf{y \in \partial B} d(\mathbf x , \mathbf y ), \end{aligned}$$
(8)

and since \(\partial B\) is a sphere this reduces to

(9)

The directed Hausdorff distance finds the point in \(\partial Z\) which has the greatest distance to \(\partial B\) and then reports that distance. The point in question, in our case, will either be the point in \(\partial Z\) furthest away from the origin (maximizer of \( \Vert \mathbf x \Vert -r\)), or the point closest to the origin (maximizer of \(r-\Vert \mathbf x \Vert \)).

The point in \(\partial Z\) furthest away from the origin has to be a vertex of Z, and since all vertices of our zonohedron have one of two distances from the origin (see Fig. 3(a)), only two vertices need to be considered when maximizing \(\Vert \mathbf x \Vert - r\).

To find the point in \(\partial Z\) closest to the origin, we notice that Z has four different types of faces (see Fig. 3(a)). Three of those (two kinds of octagons and a hexagon) have the point closest to the origin in the center of the face. It can be shown that the last type of face (quadrilateral) never contains a point either furthest or closest to the origin, so these need not be considered. In conclusion, we can reduce (9) to

$$\begin{aligned} d_\text {h}(\partial Z, \partial B) = \max \left\{ \max _{k=1,2} \{ \Vert \mathbf p _k\Vert - r \}, \max _{k = 3,4,5}\{ r-\Vert \mathbf p _k\Vert \} \right\} , \end{aligned}$$
(10)

where \( \mathbf p _k \) are two vertices of Z and centers of three faces of Z as indicated in Fig. 3(b).

Each of five points \(\mathbf p _k\) can be expressed as a linear combination of generators from \(\mathcal {V}\). Due to our grouping of the generators, we can find five matrices \(\mathbf M _k\) such that

$$\begin{aligned} \mathbf p _k = \mathbf M _k \mathbf a \,, \end{aligned}$$

where \(\mathbf a =(a_1,a_2,a_3)^T\). For the points indicated in Fig. 3(b) we identified the five matrices as

(11)
Fig. 3.
figure 3

The symmetry of the truncated rhombicuboctahedron. (Color figure online)

Full size image

To find the values of \( \mathbf a = (a_1, a_2, a_3)^T \) resulting in the best zonohedral approximation, we solve the following optimization problem

$$\begin{aligned} \begin{aligned} \min _\mathbf{a }&\max \left\{ \max _{k=1,2} \{ \Vert \mathbf M _k \mathbf a \Vert - r \}, \max _{k = 3,4,5}\{ r-\Vert \mathbf M _k \mathbf a \Vert \} \right\} \\ \text {s.t. }&0 \le a_j. \end{aligned} \end{aligned}$$
(12)

Rewriting this on epigraph form results in a second order cone problem which is convex [4]. Thus, using a branch and bound method [2], we can find integer solutions which are globally optimal but not necessarily unique.

It should be noted that one can additionally constrain the zonohedron so it is either contained within, or contains, the sphere of radius r, without losing convexity. Furthermore, it is useful to add the constraint that \(a_1 \ge 1\), otherwise the resulting structuring element is not guaranteed to be solid – instead it will have a checkerboard pattern.

3 Results

We first assess the quality of the zonohedral approximation by comparing it with a sphere. Figure 4 shows the result of dilating a single voxel with a discretized sphere and with a corresponding zonohedral approximation. For small radii the discretized sphere and its zonohedral approximation are equal, but for larger radii the difference is more apparent. This is also illustrated in Fig. 5, which shows that the directed Hausdorff distance increases linearly with the radius. When the zonohedron is constrained to be inside or outside the sphere the distance roughly doubles. The jagged lines are due to constraining the zonohedron to have integer side lengths.

Fig. 4.
figure 4

Discretized spheres (top) and their zonohedral approximations (bottom). The radii are (left to right): 3, 6, 9, and 12 voxels.

Full size image
Fig. 5.
figure 5

Directed Hausdorff distances between zonohedral approximation and a sphere at different radii. Distances are shown for the best possible approximation, and two constrained cases.

Full size image
Fig. 6.
figure 6

Morphological opening of a 3D X-ray CT scan of cement containing spherical bubbles. Data has been inverted and thresholded to make the structures of interest more visible. The two top rows show a 3D rendering of the full volume and a zoom on a small region. The bottom rows show a horizontal slice of the full and zoomed volume. Results are shown using both a discretized sphere of radius 16 and its zonohedral approximation as the structuring element.

Full size image
Fig. 7.
figure 7

Morphological opening of a 3D X-ray CT scan of stone wool containing spherical impurities. The data has been thresholded to enhance structures of interest. The two top rows show a 3D rendering of a subset of the full volume and a zoom on a small region. The bottom rows show a horizontal slice of the full and zoomed volume. Results are shown using both a discretized sphere of radius 16 and its zonohedral approximation as the structuring element.

Full size image
Fig. 8.
figure 8

Radius distributions for the detected objects shown in Figs. 6 and 7.

Full size image
Fig. 9.
figure 9

Run times for dilation with a spherical structuring element using a naive implementation and a zonohedral approximation. Benchmarks were run on a volume of unsigned 8-bit integers with size \( 500 \times 500 \times 500\).

Full size image
Fig. 10.
figure 10

Comparison of run times for different implementations of dilation with a spherical structuring element shown in a log-log plot. Benchmarks were run on a volume of size \(100 \times 100 \times 100\).

Full size image
Fig. 11.
figure 11

Comparison of run times for dilation with a spherical structuring element using the zonohedral approximation and the one provided by Avizo. Benchmarks were run on a volume of size \(500 \times 500 \times 500\).

Full size image

We now apply the developed method to two 3D X-ray \(\mu \)CT scans containing large spherical structures. The first is a porous cement sample which contains spherical bubbles of varying size. The second is a network of stone wool fibers, that holds a number of spherical impurities. The volumes consist of \( 2048 \times 2048 \times 2048 \) for the cement sample and \( 958 \times 1011 \times 3642 \) voxels for the stone wool sample. A morphological opening was performed on both volumes using a discretized sphere of radius 16 and a zonohedral approximation, with results shown in Figs. 6 and 7. The zonohedron was constrained to be within the sphere of radius 16, to ensure it could be contained in any object which contains the discretized sphere. With the discretized sphere, the operation took 40 min for the cement sample and 15 min for the stone wool sample. Using the zonohedral approximation, the operation took, respectively, 2 min and 1 min. All experiments were performed on a machine with an Intel® Xeon® E5-2637 v3 CPU, and NVIDIA TITAN X GPU, and 256 GB of RAM. Visually, the results are quite similar; the main difference lies with small objects as illustrated in the bottom right of Fig. 6(c). Furthermore, since the zonohedron is constrained to be inside the discretized sphere it introduces additional objects, as is apparent in the top row of Fig. 7. These objects will mainly be those with a radius similar to that of the structuring element since if any part of the structuring element is outside an object it is removed. To quantify the effect of this, an ellipsoid was fitted to each object (after thresholding) and a radius was then estimated from the length of the major principal axis. From the results in Fig. 8 it is clear that the radius distributions only differ for small objects, as expected. Thus, if the shape of the structuring element is of critical importance, care must be taken before applying the zonohedral approximation.

Finally, we evaluate the computational performance of our method. In order to make the method viable for large volumes we made a GPU implementation of the van Herk/Gil-Werman algorithm [8, 20] based on the description by Domanski et al. [6]. This serves as the basis of the benchmark. First, we compare the run time when using a naive GPU implementation of dilation with a discretized sphere and a zonohedral approximation using the van Herk/Gil-Werman algorithm. Results for different radii are displayed in Fig. 9 and show the expected cubic scaling for the naive method and constant run time for the approximation.

Next, we compare with other known image processing frameworks – both open source and commercial. Since ITK provides a way to do dilations/erosions with line segments, the zonohedral approximation was also implemented in terms of these. Figure 10 shows that all except Avizo and the zonohedral approximation scale with the size of the structuring element. However, for larger volumes Avizo does scale with the radius which is evident from Fig. 11. The reason for ITK’s improvement when using integers is that it changes to a moving histogram algorithm, such as the one proposed by Droogenbroeck et al. [19]. Still, the GPU based zonohedral approximation achieves the best performance out of the tested alternatives.

4 Conclusion

We have presented an approach for approximating a spherical structuring element with a zonohedron and how to compute its generator vectors by solving a convex optimization problem with integer constraints. Morphological operations with the zonohedral approximation can easily be implemented in any image processing framework which supports dilation/erosion with line segments. This allows morphological filtering to be run efficiently with significantly reduced computation times compared to other frameworks. Finally, the run times are independent of the size of the structuring element, which further increases the performance benefit for large structures.