DISTANCE FUNCTION AND EXTENSION IN NORMAL DIRECTION FOR IMPLICITLY DEFINED INTERFACES

. In this paper we present a novel application of extrapolation procedure for three popular numerical algorithms to compute the distance func- tion for an interface that is given only implicitly. The methods include the fast marching method [8], the fast sweeping method [10] and the linearization method [3]. The extrapolation procedure removes the necessity of a special initialization procedure for the grid nodes next to the interface that is used so far with the methods, thus it represents a natural extension of these methods. The extrapolation procedure can be used also for an extension of a function that is deﬁned only locally on the interface in the direction given by the gradient of distance function [2].


1.
Introduction. When using level set methods to describe an evolving curve or surface [9,6], the interface is given only indirectly as a zero level set of some function. Consequently, to obtain the distance function to such interface one needs to solve eikonal equation with Dirichlet boundary conditions defined on implicitly given boundary.
To avoid any explicit reconstruction of the boundary we propose here to use extrapolation procedure near the boundary for three popular numerical algorithms: the fast marching method [8], the fast sweeping method [10], and the linearization of eikonal equation [3]. In such a way, one can skip any initialization of numerical values in the grid nodes next to the interface as it is realized for these methods so far. For the fast marching method with the extrapolation procedure near the interface we can report an improvement when compared with the fast marching method in [1] for the computations of the distance function of implicitly defined interface.
Moreover, these numerical methods using the extrapolation near interface can be applied also for so called extension in normal direction to the interface [11,1,2]. Such procedure is required when some quantity known only on the interface must be extended away from the interface to a surrounded computational domain. The extension can be obtained by solving a stationary linear advection equation with Dirichlet boundary conditions defined on the implicitly given interface.
We note that our aim is not to compare the available numerical algorithms for the computations of distance function or for the solution of stationary linear advection equation with respect to efficiency and accuracy. We recognize that each of them is preferable in some special situations, see e.g. [5].
The paper is organized as follows. In section 2 we present the related mathematical models. In section 3 the numerical scheme of Rouy-Tourin [7] is extended by extrapolation procedure for the grid nodes next to the interface. In section 4 some numerical experiments are presented.
2. Mathematical models. Let Γ be a closed interface contained in D ⊂ R 2 . The domain enclosed by Γ will be denoted by Ω, i.e. Γ ≡ ∂Ω. We suppose that the interface is given only implicitly as the zero level set of some function φ : D → R such that φ(x) < 0 for x ∈ Ω and φ(x) > 0 for x ∈ D \ Ω, so The implicit definition (1) of Γ using the zero level set of φ is typical when some evolving interface is described using level set methods. We avoid in our approach any explicit reconstruction of Γ.
The distance function d of Γ is defined by where | · | denotes the Euclidean norm. The gradient of distance function plays an important role in level set methods. One can show [6] that if for x ∈ D there exists a unique closest point γ ∈ Γ such that d(x) = |x − γ| then ∇d(x) is well defined and |∇d(x)| = 1. The existence of ∇d(x) or the property |∇d(x)| = 1 need not to be true when the closest point on Γ for x is not unique.
Following [9,6], one can find d as a viscosity solution of eikonal equation The problem (3) can be solved independently for two subdomains of D, once for x ∈ Ω and once for x ∈ D \ Ω. For the first subdomain one has Γ ≡ ∂Ω, for the second one the boundary is given by Γ ∪ ∂D, but no boundary conditions are required on ∂D. The equation (3) can be written formally as a stationary advection equation with The formulation (4) can be used for a linearization of eikonal equation [3]. Once the distance function d is available and ∇d can be computed, one can utilize the linear advection equation for an unknown function s, The interpretation of (5) can be seen [2] as an extension (extrapolation) of known values S (that are given only on the interface Γ) to the values s defined in the whole domain D. The extension using (5) is constant in normal direction to the interface, and the value s(x) is equal to S(γ) if γ ∈ Γ is the unique closest point to x. Again, the problem (5) has to be solved for two subproblems.
3. Numerical methods. We present a numerical scheme of Rouy-Tourin proposed in [7] that is popular to use when solving the eikonal equation (3) on rectangular grids. The method seems to give the smallest error among several first order accurate schemes [8,10]. We use a standard notation for rectangular grids. For a simplicity of notation let D be a square (0, L) 2 . Let N > 0 be a chosen number and h = L/N . We aim to find the values d ij that approximate the exact distances d(x i , y j ) with x i = ih and y j = jh, i, j = 0, 1, . . . , N by numerical solution of (3). Analogously, the values s ij ≈ s(x i , y j ) will be found by solving (5) numerically.
The main idea is to choose an appropriate approximation of gradient ∇d ij ≈ ∇d(x i , y j ). It can be obtained by finite difference scheme where the "upwind" differencing shall be used: and analogously When using (7) and (8) for the nodes lying on ∂D, one has to simply skip the unavailable difference in (7) or (8).
Before using the approximation (6) with (7) and (8) for the fast marching [8] or the fast sweeping method [10], the eikonal equation (3) is written in the form |∇d| 2 = 1. Consequently, the numerical scheme to find the approximative distance function can be written in the form When solving the eikonal equation using (4) the scheme takes the form Note that when linearizing (10), one must insure that no division by zero occurs in (10), see later the related discussion in section 3.2. Finally, to solve (5) numerically we use standard first order accurate upwind differencing to obtain 3.1. Extrapolation procedure. For implicitly defined interfaces (as described in section 2) one has to use (9) -(11) for two subdomains. To identify to which subdomain the grid points (x i , y j ) belong, we denote φ ij := φ(x i , y j ). Clearly, (x i , y j ) ∈ Ω when φ ij < 0, and (x i , y j ) ∈ D \ Ω when φ ij > 0. If φ ij = 0 one can use directly the Dirichlet boundary conditions on Γ for d ij and s ij , therefore we exclude these values later from unknowns in algebraic equations.
We say that the grid node (x i , y j ) is next to the interface in x direction if at least one of the following two inequalities is valid: Analogously, we say the node (x i , y j ) is next to the interface in y direction if at least one of following inequalities is true: We comment the computations of distance function for the grid nodes next to the interface. In previous efforts [1,10], the values d ij for the nodes next to the interface are computed by some kind of "brute force" method (or its approximation) using directly the definition (2). This might not always be straightforward to do and needs not to be compatible with the computations of d ij obtained by the numerical scheme (9).
Here we propose to use an extrapolation procedure near the interface, so the values d ij for the grid nodes next to the interface are computed by the numerical scheme (9).
To do so we characterize the position of a point on the interface between two neighboring grid points for which φ changes its sign. Let (x i , y j ) be a grid node next to the interface in x direction. The required point, say x i−α j or x i+α j , will be identified with some α i−1j ∈ (0, 1) or α i+1j ∈ (0, 1) such that Analogously the values α ij±1 ∈ (0, 1) and the point y i j±α shall be treated.
The values α i±1j will be determined from the linear interpolation of φ ij and φ i±1j by requiring that 0 = α i±1j φ i±1j + (1 − α i±1j )φ ij and analogously for α ij±1 . Clearly, Note that the values α i±1j and α ij±1 are specific for each grid node (x i , y j ) that we do not further emphasize in their notation.
The idea of extrapolation near the interface is to use the same procedure to compute the values that are not available in (7) and (8) for the grid nodes next to the interface. Particularly, the value u i±1 j for the node next to the interface in x direction can be extrapolated from the values u ij and u i±α j := u(x i±α j , y j ) with the latter one given by Dirichlet boundary conditions on Γ. In our case, Substituting (15) to (7) and (8) we obtain for the nodes next to the interface in x direction and analogously for the nodes next to the interface in y direction Summarizing, in numerical schemes (9) -(11) one has to use the definitions (16) and (17) for the grid nodes next to the interface, for all other grid nodes the definitions (7) and (8) shall be used.
In next section we introduce numerical methods to solve the algebraic equations (9) -(11).

3.2.
Numerical solution of algebraic systems (9) - (11). Before introducing three algorithms to solve the algebraic system of equations (9) -(11), we discuss the possible forms of these equations in details.
Firstly, we consider (9) used for the fast marching and the fast sweeping method. In a generic case when the both components of ∇d ij computed from (7) and (8) are nonzero, the equation (9) has three unknowns and can be written in the form The choice of signs ± in (18) is determined from the valid cases in conditional definitions (7) or (8). If one component of ∇d ij is zero, the equation (9) has only two unknowns and takes a simpler form Clearly, the correct solutions of quadratic equations (19) are d ij = d i±1 j + h and d ij = d i j±1 + h. Now we discuss the form of (9) for the nodes next to the interface. Firstly, if only one inequality in (12) -(13) is valid, then one component of ∇d ij is computed from (16) or (17) and the other one from (7) or (8). We recognize then two cases. The first case is when the component of ∇d ij computed from (7) or (8) is nonzero, then the equation (9) turns to In the other case one obtains directly Of course, the values d ij in (21) represent the exact distance between x i and x i±α and the exact distance between y j and y j±α . Finally, if the both component of ∇d ij are computed from (16) and (17), the value d ij can be determined from (9) directly by Again, the choice of plus or minus sign in α is given by the corresponding valid cases in conditional definitions (16) and (17). The value d ij in (22) represents the exact distance of (x i , y j ) to the line segment between two points (x i±α , y j ) and (x i , y j±α ). Note that if (9) turns to any particular form (18) -(22), one has that |∇d ij | = 1. We note that the exact distances in (21) and (22) for the grid nodes next to the interface are fixed in the fast marching method as described in [1], and the form (20) of numerical scheme is never used in [1]. We can show later in section on numerical experiments that using the extrapolation procedure (20) -(22) with the fast marching method one can improve slightly the results for the chosen examples comparing with the method given in [1]. Now we are ready to present the solution procedures to solve the algebraic system of equations (9). The basic feature of fast marching and fast sweeping method is to suggest such consecutive steps in solving (9) that one has to solve in each step only one scalar quadratic equation with single unknown d ij . It means that if any particular form (18) -(20) has to be solved then the required neighbor values d i±1j or d ij±1 are known at that moment. It is interesting to remark that the scheme (9) can take different particular form (18) -(22) during different steps of solution procedure.
We begin with the introduction of fast marching method [8] using the extrapolation procedure. We skip the implementation details like a heapsort technique that can be found in literature and we emphasize the algebraic part of the algorithm.
The basic idea is to label all grid nodes by one of the marks [1]: Accepted, Close and Far. At the beginning of fast marching method with the extrapolation procedure the set of Accepted nodes is empty, the set of Close nodes is given by all grid nodes next to the interface, and the rest of nodes are labeled as Far. Note that in the previous implementations of fast marching method [8,1], the set of Accepted nodes consists of all grid nodes next to the interface for which the values of distance function are computed by some kind of brute force method.
The basic idea of fast marching method is to use the conditional definitions (7) and (8) for the node (x i , y j ) only with the neighbors (x i±1 , y j ) and (x i , y j±1 ) that are labeled as Accepted. This can be formally viewed as if the neighbor values d i±1 j and d i j±1 for the nodes (x i , y j ) labeled as Close or Far are set to infinity. Consequently, at the initial part of algorithm, the values d ij for all Close nodes can be computed using only the direct definitions (21) or (22).
The fast marching method consists now in repeating the following step until all nodes are labeled as Accepted: • the minimal value d ij among Close nodes is fixed and the corresponding node is labeled as Accepted; • the neighbors of this node with label Far (if any) are labeled as Close; • the values for all neighbors of this node that are not yet Accepted are recomputed using (18) -(22). In such a way the fast marching method gives the final results in the finite number of steps where the number of steps equals to the number of grid nodes. In each step at most three values associated with Close nodes must be computed by solving one of a quadratic equation (18) -(22) having only one unknown variable d ij . Following [8] the larger solution of two possible solutions to the quadratic equation shall be assigned to d ij , see also the text after (19).
The fast sweeping method [10] is iterative method that is based on nonlinear Gauss-Seidel method with four (in 2D case) different orderings ("sweepings") of updates for unknown values d ij . The initial guess is d 0 ij = C where C is a large constant value that can not be reached by the distance function for the given interface and the domain D. The consecutive iterations d The detailed description for the solution of quadratic equations is given in [10]. It is important to note that the value d k+1 ij is updated in each iteration only if d k+1 ij < d k ij . Standard criteria can be used based on residuum to decide when to stop the iterations. It can be shown that in a theory a finite number of iterations is necessary to reach the convergence [10].
In the last case we consider the solution procedure for the linearization method based on the linearization of numerical scheme (10). Let d (23) The algebraic system (23) is linear and can be solved by the Gauss-Seidel iterations in four alternating directions. In general it may happen that for some inappropriate is undefined if using (10). Therefore the method (23) is suitable only when some good initial guess d 0 ij for the approximative distance function is available.
Finally we comment the solution of linear advection equation (5). Clearly, the proposed numerical scheme (11) is in a character very similar to the left hand side of (23), so one can again use the Gauss-Seidel iterative methods in four alternating directions for the solution of linear algebraic system. No nonlinear (quadratic) equations have to be solved, therefore we do not consider the fast sweeping method for (5).
The fast marching method to solve (5) is considered in [1] where, in fact, the both equations (3) and (5) are solved in parallel using (9) and (11). In each step of the method, the scheme (11) is a linear equation with a single unknown s ij , therefore it is trivial to solve. The order of nodes (x i , y j ) in which the unknowns s ij are determined is given by the order in which the Close nodes are labeled as Accepted in the computations of d ij .
In next section we apply all solution procedures for some standard benchmark examples.
4. Numerical experiments. The purpose of following numerical experiments is to illustrate that the extrapolation procedure for the nodes next to the interface does not decrease the order of accuracy for three solution algorithms. We begin with the computations of approximative distance functions to a smooth interface (a circle), an interface with corners (a square) and a nontrivial interface ( a quatrefoil). For the circle and the square the exact solution is available, for the quatrefoil it is computed by an approximative brute force method. The brute force method is approximative in the sense that the distance d(x i , y j ) to a polygonal interface is computed as the minimal distance to the points that defines the polygon.
In Table 1 we present the discrete L 1 error for four consecutively refined grids computed with three methods with the extrapolation procedure for three interfaces. The experimental order of convergence (EOC) is computed for each method and each interface. The EOC takes approximately the same value for each method, with the value around 1 for the circle and the quatrefoil, and the value around 0.8 for the square. In Figure 1 we plot numerical solutions where the values of computed distance functions inside of interface are multiplied by −1 to obtain so called signed distance functions [9]. In Figure 1 one can see a visual comparison of numerical solution obtained by the fast marching method with the exact solution (for the circle and the square) or the distance obtained by approximative brute force method. We do not compare visually the numerical solutions obtained by the fast sweeping method, because practically no difference is visible. Next we present a comparison when using the fast marching method with an approximative brute force initialization for the grid nodes next to the interface as described in [1] and the fast marching method with the extrapolation procedure (20) -(22). As we commented in section 3.2 the method in [1] uses only the forms (21) and (22) of (9) for the grid nodes next to the interface.
For the computations of distance to the square the identical numerical errors are obtained. For the computations with the circle and the quatrefoil an improvement of results can be reported for the fast marching method with the extrapolation procedure as given in Table 2.  Table 2. The discrete L 1 errors (the magnitude is 10 −3 ) and the EOCs for the distance function to the circle (the columns 2 − 5) and the quatrefoil (the columns 6 − 9). The fast marching method with the extrapolation method is given in the columns with E in headers, the errors obtained by the fast marching method in [1] are given in the columns with [1] in the headers.
Finally, an extension of a function S = x of which the values are taken only on a circle is computed in normal direction and compared with available exact solution [4]. The numerical solutions obtained by the fast marching method and the linearization method are compared in Table 3 and in Figure 2.  Table 3. The discrete L 1 errors and the EOCs for the extension of function S = x defined only on a circle for the fast marching method (M) and the linearization method (L).