Cosine-Pruned Medial Axis: A new method for isometric equivariant and noise-free medial axis extraction

Blum [11] first introduced the medial axis as an analogy of a fire propagating with uniform velocity on a grass field. The field is assumed to have the form of a given shape. If one “sets fire” at all boundary points, the medial axis is the set of quench points. There are other equivalent definitions of the medial axis. In this work we use a geometric definition as follows:

The medial axis, together with the associated radius of the maximally inscribed ball, is called the medial axis transform ($\mathbf{MAT}(\Omega)$). The medial axis transform is a complete shape descriptor, meaning that it can be used to reconstruct the shape of the original domain. In some work, $\mathbf{MA}$ and $\mathbf{MAT}$ are also referred to as shape skeletonization. Figure 1 shows an example of a 2D shape and its medial axis as the center of maximal discs. In $\mathbb{R}^{3}$ definition 1 may result in a 2-dimensional medial axis sometimes called the medial surface. We will restrict our examples and analysis to only one-dimensional medial axes.

There are three primary mechanisms to compute the $\mathbf{MA}$: 1) layer by layer morphological erosion also called thinning methods, 2) extraction of the medial axis from the edges of the Voronoi diagram generated by the boundary points, and 3) detection of ridges in distance map of the boundary points. In digital spaces, only an approximation to the “true medial axis” can be extracted.

When using thinning methods  [16, 27, 26, 32], points which belong to $\Omega$ are deleted from the outer boundary first. Later, the deletion proceeds iteratively inside until it results in a single-pixel wide medial axis. The medial axis by thinning can be approximated in terms of erosion and opening morphological operations [55]. Thinning methods are easy to implement in a discrete setting, but they are not robust to isometric transformations.

The most well-known algorithm for thinning skeletonization is perhaps the Zhang Suen [55] algorithm. However, other approaches have been developed using similar principles [52, 26, 32], mainly focused on parallel computation.

Another way to estimate the medial axis works by computing the Voronoi Diagram (VD) of a polygonal approximation of the object’s contour. The contour is expressed as line segments in 2D or a polygonal mesh in 3D. The seed points for the VD are the vertices of the polygonal representation. The medial axis is then computed as the union of all of the edges $e_{ij}$ of the VD, such that the points $i$, and $j$ are not neighbors in the polygonal approximation [33].

Voronoi skeletonization methods preserve the topology of $\Omega$. However, a suitable polygonal approximation of an object is crucial to guarantee the medial axis’ accuracy. Noise in the boundary forms convex areas in the contour, which induce spurious branches on the medial axis. In general, the better the polygonal approximation, the closer the Voronoi skeleton will be to the real medial axis. Nevertheless, this is an expensive process, particularly for large and complex objects [36].

The most common methods to extract the medial axis are those based on the distance transform. Within these methods, the medial axis is computed as the ridges of the distance transform inside the object [3, 40, 2, 22, 14, 35, 41]. This interpretation of the medial axis follows definition 1, because the centers of the maximal balls are located on points $x$ along the ridgeline of the distance transform, and the radius of the balls correspond to the distance value at $x$.

When computed in a discrete framework, distance-based approaches suffer from the same isometric robustness limitations as thinning and Voronoi methods  [36]. Moreover, noise in the contour produced by a low discretization resolution directly affects the medial axis’ computation by introducing artificial ridges in the distance transform and, consequently, spurious branches.

The medial axis’ sensitivity to boundary noise limits its applications to real-life problems [10]. Even negligible boundary noise can cause spurious branches, as shown in Figure 2.

One strategy for removal of spurious branches consists of computing $\mathbf{MA}(\Omega)$ as the medial axis of a smoother version of the shape, $\Omega^{\prime}$ [38, 31]. The main disadvantage of this approach is that, in most cases, the resulting medial axis is not a good approximation. Additionally, the smoothing of $\Omega$ can potentially change its topology. Miklos et al. [30] introduced a slightly different approach they call Scaled Axis Transform (SAT). The SAT involves scaling the distance transform and computing the medial axis of the original un-scaled shape as the medial axis of the scaled one. However, in [35], the authors show that the SAT is not necessarily a subset of the medial axis of the original shape. Instead, they propose a solution that guarantees a better approximation by including additional constraints on the scaled distance transform.

Another method to overcome the noise sensitivity limitation of the medial axis is spurious branch pruning. Pruning methods are a family of regularization processes incorporated in most medial axis extraction algorithms [43]. Effective pruning techniques focus on different criteria to evaluate the significance of medial axis branches. Thus, the decision is whether to remove the branch (and its points) or not. We can say that pruning methods are adequate if the resulting $\mathbf{MA}$ is noticeably simplified while preserving the topology and its ability to reconstruct the original object. Most pruning methods rely on ad hoc heuristic rules, which are invented and often reinvented in a variety of equivalent application-driven formulations [43]. Some authors apply these rules while computing all medial axis points. Others do so by removing branches that are considered useless after the computation [36, 19, 22].

One of the most popular pruning methods is Couprie et al. [14]. They consider the angle $\theta$ formed by a point $x\in\Omega$, and its two closest boundary points denoted by the set $\Pi_{\Omega}(x)$. This solution removes points from $\mathbf{MA}$ for which $\theta$ is lower than a constant threshold. This criterion allows different scales within a shape but generally leads to an unconnected medial axis. Another pruning method found in the state-of-the-art is the work of Hesselink et al. [22]. They introduce the Gamma Integer Medial Axis (GIMA), where a point belongs to the medial axis if the distance between its two closest boundary points is at least equal to $\gamma$.

Many pruning methods rely on the distance transform $\mathcal{D}_{\Omega}(x)$. The distance transform acts as a generator function for the medial axis, such that points $x\in\mathbf{MA}$ if and only if they satisfy some constraint involving their distance to the boundary. However, other authors have proposed alternative generator functions in their pruning strategies [21, 4].

For example, [21] and [4] introduce what they call Poisson skeletons by approximating $\mathcal{D}_{\Omega}(x)$ as the solution of the Poisson equation with constant source function. Poisson skeletons rely on a solid mathematical formulation. Among other concepts, they use the local minimums and maximums of the curvature of $\delta\Omega$. However, when such methodology is applied in a discrete environment, many spurious branches appear due to the need to define the length of a kernel size to estimate these local extreme points.

Let us represent $\Omega$ as a square binary image $\mathcal{I}:\mathbb{N}^{2}\mapsto\left\{0,1\right\}$ with a resolution of $M\times M$ pixels. We start the computation of the CPMA by estimating a set of smoothed versions of the $\mathcal{I}$ via the Discrete Cosine Transform (DCT) and its inverse (IDCT):

where $(u,v)$ are coordinates in the frequency domain, and $(x,y)$ are the spatial coordinates of the Euclidean space where $\Omega$ is defined. The values of $C_{u}$ and $C_{v}$ are determined by:

The DCT is closely related to the discrete Fourier transform of real valued-functions. However, it has better energy compaction properties with only a few of the transform coefficients representing the majority of the energy. Multidimensional variants of the various DCT types follow directly from the one-dimensional definition. They are simply a separable product along each dimension.

Let us now denote by $\mathcal{I}^{(i)}$ with $i=1,2,...,M$ the reconstructions of $\mathcal{I}$ using only the first $i$ frequencies as per equation 2. We seek to obtain a $\mathcal{F}_{\Omega}$ acting as a sort of probability indicating how likely it is for a point $x$ to be in the medial axis of $\Omega$. Points on relevant branches will appear regularly in the smoothed shapes’ medial axis, resulting in high score function values. In contrast, spurious branches will only appear occasionally, resulting in low values.

The score function is defined for all $x$ in the domain of $\mathcal{I}$. Notice that we represent $\mathbf{MAT}$ as another binary image where $\mathbf{MAT}(x)=1$ only when $x$ belongs to the medial axis. Finally, with $\mathcal{F}_{\Omega}$, we have all the elements to present our definition of the CPMA.

We empirically set the value of the threshold to $\tau=0.47$. However, we conducted an additional experiment to show that this value is consistent across different shapes.

Although the CPMA results in a noise-free $\mathbf{MA}$, there is no restriction in its formulation to force the CPMA to create a connected medial axis. We solve this by finding all the disconnected pieces of the CPMA. Later, we connect them using a minimum distance criterion $g(x_{i},x_{j})$, where $x_{i}$ and $x_{j}$ are endpoints of two distinct pieces. However, neither the Euclidean distance nor the geodesic distance are suitable criteria because they lead to connections between nodes that do not follow the medial axis (See figure 4). We instead connect $x_{i}$ and $x_{j}$ with a minimum energy path using an energy function $E_{\Omega}$. We must guarantee that $E_{\Omega}(x)$ has high values when $x$ is close to $\delta\Omega$ and low values when $x$ is close to the medial axis. This way, we enforce the paths to be close to the $\mathbf{MAT}$. We call the result of this strategy the Connected CPMA (C-CPMA). In section 3.3, we provide details of $E_{\Omega}$ computation.

The distance transform-based medial axis depends only on the shape $\Omega$, not on the position or size in the embedding Euclidean space. Therefore the medial axis should be equivariant under isometric transformations.

Moreover, the CPMA depends primarily on $\mathcal{F}_{\Omega}$, which also holds the isometric equivariant property.

However, in a discrete domain, this equivariance is only an approximation because points on both $\Omega$ and $\mathbf{MAT}(\Omega)$ are constrained to be on a fixed, regular grid. In a continuous domain, it is easy to demonstrate that the cosine transform has exact isometric equivariance.

To compute the CPMA enforcing the connectivity, we create a lattice graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$. A point $p$ in the domain of $\mathcal{I}$ is a node of $\mathcal{G}$, if and only if $p\in\Omega$. The node $p$ shares an edge with all its neighbors in the lattice only if the neighbors are also inside $\Omega$. We used an $8$-connectivity neighborhood in 2D and a $26$-connectivity neighborhood in 3D.

In order to determine the minimum energy path between pairs of pixels/voxels, we compute the minimum distance path inside $\mathcal{G}$ using Dijkstra’s algorithm. We assign weights to every edge with values extracted from $E_{\Omega}$. Given $(x,y)\in\mathcal{E}$, we compute the energy of every edge as follows:

This method guarantees connectivity, but it is inefficient because of the minimum energy path’s iterative computation. We sacrifice performance in favor of connectivity. We include the pseudo-code to compute the CPMA and the C-CPMA in algorithms 1 and 2.

The CPMA only relies on one parameter, $\tau$. The value of $\tau$ is the threshold that determines whether a point of $\mathcal{F}_{\Omega}$ is a medial axis point. We empirically set the value of $\tau=0.47$. However, in section 5, we present the result of an additional experiment to show how sensitive the CPMA is to different threshold values.

Another essential consideration when computing the CPMA is the maximum frequency used to reconstruct the original shape through the IDCT. Using less than $M$ frequencies enables a faster computation of the CPMA without losing accuracy. We found that using frequencies greater than $\frac{M}{2}$ does not yield significant improvement for the CPMA.

We chose three extensively used datasets of 2D and 3D segmented objects to evaluate our methodology on medial axis extraction robustness. These datasets are part of the accepted benchmarks in literature, enabling us to compare our results.

To compare the robustness of a medial axis extraction method, we adopt an evaluation strategy similar to [13]. Consequently, we measure the similarity between the medial axis of a shape $\Omega$ and $\Omega^{\prime}$. The shape $\Omega^{\prime}$ derives from a “perturbation” of $\Omega$. We are interested in evaluating how well our methodology responds to induced noise on the contour/surface. We are also interested in assessing how stable the CPMA responds in the presence of rotations of $\Omega$ to test its isometric equivariance.

We employ the Hausdorff distance ($d_{H}$), and Dubuisson-Jain dissimilarity ($d_{D}$) as metrics between shapes. The Dubuisson-Jain similarity is a normalization of the Hausdorff distance [15], which aims to overcome $d_{H}$ sensitivity to outliers. The Dubuisson-Jain similarity between point sets $X$ and $Y$ is defined as:

with

We must first choose a strategy to induce noise to the boundary to evaluate the noise sensitivity. We use a stochastic approach denoted by $\mathcal{E}$, where a random number of points $p\in\delta\Omega$ are deformed by a vector $v$ in the direction orthogonal to the boundary, with a deformation magnitude that is normally distributed, $|v|\sim\mathcal{N}(p,1)$. This noise model is recursively applied $n$ times to every shape in our datasets. We denote as $\mathbf{MAT}(\mathcal{E}(\Omega,k))$ the medial axis of a shape $\Omega$ after applying the noise model $k$ times. In our experiments, we used $k=1,...,20$ noise levels.

For every object in each dataset, we compared the medial axis $\mathbf{MAT}(\Omega)$ with the noisy versions $\mathbf{MAT}(\mathcal{E}(\Omega,k))$ to determine how sensitive a method is to boundary noise.

The medial axis is ideally an isometric equivariant shape representation so that $\textrm{MA}(R(\Omega))=R(\textrm{MA}(\Omega))$. Due to sampling factors, this relationship is only an approximation. However, we can measure the equivariance by comparing $\mathbf{MAT}(R(\Omega))$ with $R(\mathbf{MAT}(\Omega))$. The more similar they are, the more equivariant the method.

Because the translation equivariance is trivial, we evaluate isometric equivariance only with rotations of $\Omega$. We do not use reflections because the properties of rotation matrices we use in this study can be extended to reflection matrices. In our experiments, 2D rotations are counterclockwise in the range $[0,90]$ degrees around the origin. We use $30$ rotations at regular intervals. In 3D, we use a combination of azimuthal ($\theta\in[0,90]$) and elevation ($\phi\in[0,90]$) rotations around the origin. The angles $\theta$ and $\phi$ take values each $18$ degrees.

We chose seven of the most relevant methods in the scientific literature to compare with CPMA extraction results. Each method was selected based on a careful review of the state-of-the-art. These methods illustrate the variety of approaches that authors employ to prune the medial axis. Notice that the first method we used in our study is the un-pruned $\mathbf{MAT}$.

Table 1 summarizes all of the methods in our comparative study. In many cases, the performance of a pruning method depended on its parametrization. We selected parameters for all of the methods following the best performance parametrization reported in the state-of-the-art.

We compared the stability of the CPMA under boundary noise against other approaches in Table 1. We conducted our experiments on Kimia216 and the Animal2000 Dataset for 2D images. Additionally, we used a set of three-dimensional triangular meshes from the Groningen Benchmark for 3D experimentation.

For our noise sensitivity experiments, we applied $20$ times the noise model $\mathcal{E}(\Omega,k)$ to every object of each dataset. We then computed their $\mathbf{MAT}$ using every method in Table 1 with different parameters. Finally, each $\mathbf{MAT}(E(\Omega,k))$ was compared with $\mathbf{MAT}(\Omega)$ using both the Hausdorff distance and Dubuisson-Jain dissimilarity. We report the per method average of each metric over all the elements of each dataset.

First, we tested our medial axis pruning method on the Kimia216 dataset and present the results in Table 2. The table shows that the CPMA and the C-CPMA are competitive against state-of-the-art methods for medial axis extraction. Our results show similar performance to two state-of-the-art methods: the GIMA and SFEMA. The CPMA and C-CPMA also performed better than Poisson Skeletons, SAT, topological thinning, and the MAT itself. For visual comparison, Figure 6 (top row) shows both Hausdorff distance and Dubuisson-Jain dissimilarity against noise level. The figure only displays the best parametrization of every method to improve visualization. It is clear that the CPMA and the C-CPMA are among the three best results when we use the Dubuisson-Jain dissimilarity metric. However, we observe a decrease in performance when we the Hausdorff distance metric. We believe this occurs because of Hausdorff’s distance sensitivity to outliers.

The Animal2000 dataset contains nearly ten times more shapes than Kimia216. This leads to more variation between shapes, and therefore a more challenging setting. Table 3 shows similar results compared to Kimia216, confirming that the noise invariant properties of the CPMA still hold in a more robust dataset. The GIMA and the SFEMA are still the best methods measured with the Dubuisson-Jain dissimilarity, closely followed by both the CPMA and the C-CPMA. Results of using the Dubuisson-Jain dissimilarity as a metric show that the CPMA is close to methods such as BEMA and SFEMA. However, the results are not as good as when using the Hausdorff distance metric. Figure 6 (middle row) depicts the best performance for every method in comparison to ours. Our experiment’s results suggest that the CPMA noise invariant properties generalize across different datasets.

For our 3D experiments, we selected $14$ objects from the Groningen Benchmark, reflecting the most common shapes used in the literature. Each object was voxelized to a binary voxel grid with resolution $150\times 150\times 150$. This resolution offered sufficient details as well as a sufficiently low computational cost. In contrast to the 2D case, we applied $\mathcal{E}(\Omega,k)$ only $10$ times to the 3D object. We did this for two reasons: 1) to reduce computational complexity and 2) noise tends to be more extreme in 3D at the chosen resolution. The results on the Groningen dataset are shown in Table 4 and Figure 6 (bottom row). Notice that both the CPMA and C-CPMA achieved the best results among the other methods when compared with the dissimilarity measure. These results show that our methodology has noise-invariance properties, and it is stable in the presence of small surface deformation. However, the results show unusual patterns when compared with the Hausdorff distance. In fact, for some methods, the metric decreases when the noise level increases. We attribute this behavior to the outlier sensibility of the Hausdorff distance.

We complete the noise stability analysis showing examples of the $\mathbf{MAT}$ computed with our methodology, and compared to the other methods in this study, Figure 7.

We measured the dissimilarity between $\mathbf{MAT}(R(\Omega))$ and $R(\mathbf{MAT}(\Omega))$ for different instances, and different definitions of the medial axis transform. The lower this dissimilarity, the more stable the method is under rotation. The rotation sensitivity analysis on the Kimia216 dataset is summarized in Table 5 and Figure 8 (top row). The results show that the curves of the CPMA and the C-CPMA fall near the average of the rest of the methods achieving state-of-the-art performance. The results also surpassed several methods, such as Poisson skeleton, SAT, and thinning methods. Notice that when using the dissimilarity metric the CPMA, the GIMA, the SFEMA, and the BEMA form a subgroup that performs significantly better compared to the others. Moreover, the performance of these methods oscillates around a value of dissimilarity of around $1$ pixel on average. The intuition for this result is that regardless of the rotation, skeletons computed with these methods vary only at one pixel. Consequently, we can claim that they exhibit rotation equivariance.

We applied the same analysis to the Animal2000 dataset achieving similar results. In this case, the CPMA and the C-CPMA ranked third and fourth, respectively, among all methods when we used the dissimilarity metric. The results for all methods and parameters are presented in Table 6. As before, we also present a summary with the best parametrization for each method in Figure 8 (middle row) to facilitate the interpretation. Notice that due to the larger number of objects in the Animal2000 dataset, the curves for every method appear to be smoother, highlighting stability across different rotation angles and shapes.

Finally, we conducted the rotation sensitivity analysis on the 3D dataset. the results are summarized in Figure 8 (bottom row). The image shows the four 3D medial axis extraction methods we compared in our study for combinations of azimuthal and elevation angles. This figure shows how both the Hausdorff distance and the Dubuisson-Jain dissimilarity become higher when the rotation becomes more extreme, except in the case of C-CPMA. We believe this behavior is due to the connectivity enforcement mitigating the gaps in the medial axis, and reducing the metrics.

We tested our medial axis pruning method on the Kimia216 dataset. Many medial axis pruning methods depend on hyper-parameters to accurately estimate the medial axis [22, 14, 20, 35]. These parameters usually have a physical meaning in the context of the object whose medial axis they seek to estimate. Often, the parameters are distances or angles formed between points inside the object. In other work, some authors create score functions like ours, intending to use its values as a filter parameter to remove points on spurious branches of the $\mathbf{MAT}$. In most cases, however, such parameters are subject to factors like resolution or scale. Thus, we conducted another experiment to test the sensitivity of the CPMA to the pruning parameter $\tau$ at different scale factors of the input object.

Figure 9 shows the average of the Jaccard index as the reconstruction metric vs the values of $\tau$. We compared an object $\Omega$ against its reconstruction $\hat{\Omega}$ over all images in Kimia216 dataset. The figure shows how high values of $\tau$ deteriorate the reconstruction, while lower values do not prune enough spurious branches. From the figure, we can infer that values around $\tau=0.47$ offer a good trade-off between reconstruction and branch pruning. Moreover, around this value, the standard deviation reaches its minimum value suggesting optimal performance regardless of the object. Because the value of $\tau$ is stable for different scale factors, we conclude that scale does not affect the selection of the threshold.