An ADMM-based scheme for distance function approximation

A novel variational problem for approximating the distance function (to a domain boundary) is proposed. It is shown that this problem can be efficiently solved by ADMM. A review of several other variational and PDE-based methods for distance function estimation is presented. Advantages of the proposed distance function estimation method are demonstrated by numerical experiments. Applications of the method to the problems of surface curvature estimation and computing the skeleton of a binary image are shown.


Introduction
Fast and accurate estimation of the distance to a surface (the distance to a curve in 2D) is important for a number of applications including redistancing or reinitialization for level-set methods [17], wall distance models in turbulence modeling [23,29], heterogeneous material modeling in computational mechanics [6], medial axis transform and meshing [31], FEM extensions [2,16], robot navigation and mapping [9], and various computer graphics studies [10,34,35]. In particular, variational and PDEbased methods for distance function estimation are currently a subject of intensive research [3,4,11,22,24].
In this paper, we propose a variational method for accurate distance function estimation. The core of our approach consists of formulating a new variational problem, whose solution is the distance function and solving this variational problem numerically by ADMM [19]. We review several other variational and PDE-based methods for distance function estimation, such as the recent geodesics-in-heat method [11,12]. We show how to discretize the variational problem by the finite element method, and by using finite difference. In particular, we present the Matlab code for implementing our method (as well as the other variational methods) for computing the distance transform of a 2D binary image. We compare the results obtained by our approach with those of the other methods, and also show how our method can be used for estimating surface curvatures and approximating the skeleton of a 2D binary image.

Screened Poisson distance function approximations
A simple PDE-based approach to estimate d(x) consists of considering a Dirichlet boundary value problem for a screened Poisson equation where t is a small positive parameter. As shown in [30,Theorem 2.3], Thus, √ t ln[1/w(x)] approximates d(x) and the parameter t controls the approximation accuracy. This approximation of the distance function was used in [28] for detecting skeletal structures in grayscale images.
A simple heuristic behind (4) can be found, for example, in [20] and uses the Hopf-Cole transformation [13]. It is easy to check that substituting w( which for small t can be considered as a regularized version of eikonal (1).
The geodesics-in-heat method [11,12], a highly popular distance function approximation, is a close relative of the above approach extended to non-flat geometry. The method is based on the observation that w(x) and v(x) have the same level sets and construct a distance function approximation u(x) as the solution of the following minimization problem where m = ∇w/|∇w| and w(x) is the solution to (3). Practical implementation of the geodesics-in-heat method includes numerical solutions to the screened Poisson (3) and the Poisson equation corresponding to the variational problem (5).
It is also interesting that (3) combined with either (4) or (5) can serve as a warm start to iterative methods for computing the distance function.
Although both the methods, (4) and (5), are based on screened Poisson (3), to distinguish them we refer to the former as the screened Poisson method and to the latter as the heat method. In Section 4, we deal with approximate distance transforms for 2D binary images and provide the reader with simple MATLAB-based implementations of both methods.

Proposed distance function approximation
We now turn to the description of our approach. A key observation is given by the following proposition. Proof Indeed, in order to maximize φ dx, the function φ(x) has to grow as fast as possible and its gradient magnitude |∇φ| achieves at each point x its maximal allowed value. Thus, the solution to (6) satisfies |∇φ| = 1 in and φ = 0 on ∂ and, therefore, is the distance function d(x).
The constrained energy minimization problem (6) can be reformulated as the following unconstrained minimization problem where We solve (7) numerically by ADMM. The corresponding augmented Lagrangian has the form where σ is the vector of Lagrange multipliers and r > 0 is a regularization parameter. Now ADMM for (8) yields the following iterative process: The iterative approach presented here shares a lot of similarities with the numerical method proposed in [15], where ADMM is used to numerically compute the solution to a p-Poisson problem − p u = 1. The limit as p → ∞ of the solution to this p-Poisson problem, with vanishing Dirichlet boundary condition, is d(x), the distance to the boundary function [5]. The main difference with the approach proposed here lies in the second step of minimizing (8) by ADMM, where a simple projection can now be used instead of numerically computing the root of a polynomial.

Numerical experiments: distance to a polygonal mesh using FEM and application to curvature computation
We consider first the problem of computing the distance function to ∂ , a surface (curve in 2D), represented by a triangle mesh (a polygonal chain in 2D).
As described in Section 2, we use ADMM to minimize (7). This gives an iterative process, where the first step involves solving a Poisson problem of the form: − φ = f . This Poisson problem is discretized and solved by the finite elements method. The computational domain bounded by ∂ is represented by a triangle mesh in 2D or a tetrahedral mesh in 3D. Linear basis functions are used at each node of the triangulation. The solution to the Poisson problem is obtained from numerically solving a linear system A = b, where the sparse matrix A, corresponding to a discretization of the Laplacian, is the same for each iteration of ADMM and can thus be prefactored (for example with the Cholesky decomposition). Figure 1 illustrates the result obtained by our approach on a 2D polygonal domain with complex geometry. The left image visualizes the exact distance obtained by computing the minimum distance to any boundary segment. The middle image presents the distance computed by solving (7). The right image demonstrates how the relative residual error φ k+1 − φ k 2 / φ k 2 decreases with each iteration when solving (7) by ADMM. One can observe that (7) solved numerically by ADMM demonstrates a good convergence and is capable to deliver an accurate approximation of the distance function. As seen in Fig. 2, a similar performance of (7) with ADMM is achieved for 3D models. In this figure, we visually compare on a planar slice of the function obtained by numerically solving (7) and the exact distance, which is obtained by computing at a given node of the tetrahedral mesh the minimal distance to the set of triangles corresponding to the input surface ∂ .

Comparison to other approaches
In Figs. 3 and 4, we compare the distance function estimation results achieved by solving (7) with those obtained by using the heat method (3) and (5). The results are shown on a planar slice of the domains.
As seen in Fig. 3, the distance function approximation obtained from minimizing (7) delivers the best result. This is also confirmed by Fig. 4, where the point-wise relative error w.r.t. the exact distance is visualized on a planar slice of the domain. Here u(x) is the distance function approximation obtained by either (7) with ADMM, or the heat method (3) and (5).
One can notice that (7) delivers the best approximation. Some small regions with relatively high errors are likely due to the numerical solution of the ADMM

Application to surface curvature estimation
We can further use the computed distance to approximate the curvature of the boundary surface ∂ . Surface curvature estimation is important for many computer graphics and geometric modeling applications [8] and remains a subject of intensive research [21,27] (see also references therein). It turns out that the distance function from a surface contains full information about the surface curvatures. Namely, the eigenvalues of the Hessian H of the distance function d at a boundary point y are {κ 1 , κ 2 , 0} and the corresponding eigenvectors are {t 1 , t 2 , n}, where t 1 , t 2 are the principal curvature directions at y and n the outward normal to the surface at y. See, for example, [18, Section 14.6].

Fig. 4
Point-wise relative error w.r.t. the exact distance for the distance approximation obtained by the heat method [11] (two left images) and our variational problem (7) solved numerically by ADMM (two right images) visualized on a planar slice of the domain Given an approximation of the distance, computed by numerically solving (7), and sampled on a tetrahedral mesh, we first approximate the gradient of the distance function ∇d at each node of the tetrahedral mesh as the normalized weighted average of the gradients in each incident tetrahedron (the volume of the tetrahedron is used as the weight). By repeating this procedure for each component of ∇d, we obtain an expression of the Hessian of the distance function at each node of the tetrahedral mesh. We then compute the non zero eigenvalues of the Hessian at each boundary node (nodes of the tetrahedral mesh that are on the surface boundary). Figure 5 illustrates results obtained by this approach for computing the mean curvature, H = (κ 1 + κ 2 )/2 of some surfaces.

Approximate distance transforms for 2D binary images
Constructing exact and approximate distances from boundaries of objects is also important for a number of image processing applications [14], where it is usually called the distance transform. In particular, shape skeletonization, which itself has numerous applications in image processing and computer graphics [26], can be used as a quality indicator for distance function approximations.
In this section, we focus on computing approximate distance functions for the boundaries of objects represented by 2D binary images. Since the most popular way to encode 2D images consists of representing them using grids of numbers (pixels), we deal with finite difference approximations. We present MATLAB-based implementations of the screened Poisson method (3) and (4), heat method (5) of Crane, Weischedel, and Wardetzky [11,12], and our method (7). We demonstrate that typically our method produces a more accurate approximation of the distance function. Finally, we present a simple way to extract skeletons of binary shapes.

Fig. 5
Mean curvature for the bunny and the fertility models, computed from the eigenvalues of the Hessian of the distance function MATLAB allows for a very simple implementation of the screened Poisson based and heat distance function approximations discussed in Section 2.2. Since in this section we use them in our numerical experiments and comparison, we present their implementations below. First, our MATLAB implementation of (3) and (4) for approximating the distance to the shape of a binary image looks as follows: Next, our MATLAB implementation of the heat method [11,12] solving (3) and minimizing (5) reads as follows: Finally, a simple implementation of the ADMM-based minimization of (7) in MATLAB, relying on its built-in functions, is shown below In the numerical experiments presented below, we compare these three methods (the screened Poisson (3) and (4), the heat method (3) and (5), and the proposed approach (7)) for computing the distance to the boundary of binary shapes. We set t = 0.6 for both the screened Poisson and heat methods, as choosing a smaller value may lead to instabilities. We set r = 10 and use 50 iterations for our method.
In our experiments, we use two simple geometric shapes (mushroom and keyhole) and four more complex shapes from the Kimia-99 dataset [25]. Figure 6 illustrates the results obtained by these methods on the test shapes. The bottom three rows show the point-wise absolute error of the screened Poisson, the heat method, and our approach, respectively, compared with the exact distance. Quantitative results are provided in Table 1, where the RMS and maximum error for each of the test shapes is provided.
The heat method (5) demonstrates a better stability than the screened Poisson one (3) and (4). In particular, the latter collapses for some images tested above if we Fig. 6 Distance estimation on test images. First row (from top to bottom): input binary images. Second row: the exact distance from the boundary. Third row: absolute error for the screened Poisson based distance function approximations (4). Fourth row: absolute error for the heat distance function approximations (5). Fifth row: absolute error for the proposed distance function approximations (7) set t = 0.5, while the heat method demonstrates slightly better performance. On the other hand, our method (7) can also benefit from a more accurate selection of the regularization parameter r. For example, setting r = 7 makes our method the absolute RMS and maximum error winner for all the images considered above.

2D binary image skeletonization
The proposed approach for distance estimation can be used to compute skeletons of binary images. The skeleton, or medial axis, is an important shape descriptor, which was introduced more than fifty years ago [7] and which remain a topic of active research, see, for example, [1,26,32,33]. The skeleton of a 2D shape can be defined as the locus of centers of inner bitangent circles. Alternatively, given the distance function to the boundary of a shape, the skeleton can be defined as the set of inner distance function singularities. In this study, we use a smooth distance function for constructing an approximate skeleton of a binary image.
Since the magnitude of the gradient of the distance function is equal to one everywhere except at the distance function singularities, it is natural to expect that the points where the gradient of a smooth distance function is small form a thick/fuzzy version of the skeleton. Once such a thick/fuzzy skeleton is extracted, one can get a one-pixel-wide skeleton by applying MATLAB's bwmorph function to the thick/fuzzy skeleton.
In Fig. 7, we demonstrate how this gradient-of-smooth-distance-function approach works and compare it with the standard thinning-based binary image skeletonization procedure consisting of applying bwmorph to the whole binary image. Given a binary image, we start from using (7) for computing a smooth distance function u(x). Then, the fuzzy-skeleton function and its normalized version are evaluated. Here, the factor √ u(x) is used to suppress image skeletal structures appearing due to small perturbations of the image boundary. As demonstrated by the Fig. 7 Using the distance approximation (7) for skeleton extraction. Top row: our distance function approximation (7) is used to compute a gradient-based function S n (x) defined by (9). Middle row: the skeletons obtained by applying simple thresholding and thinning to S n (x). Bottom row: the skeletons obtained by thining the original binary images with bwmorph top row images of Fig. 7, S n (x) delivers a good detection of skeletal structures of the binary images. Applying simple thresholding and thinning by bwmorph yields one-pixel-wide skeletons (the middle row of Fig. 7), which are of higher quality than those generated by thinning the original binary images with bwmorph (the bottom row).

Conclusion
In this paper, we have proposed a new variational problem for the distance-fromsurface function (7) and showed how it can be efficiently solved using ADMM. We have presented the results of our numerical experiments when the problem is discretized by FEM or finite differences. We have demonstrated advantages of (7) over the heat method [11] and shown that our approach can be used for surface curvature estimation and skeletonization of 2D binary images.