Integration over curves and surfaces defined by the closest point mapping

We propose a new formulation for integrating over smooth curves and surfaces that are described by their closest point mappings. Our method is designed for curves and surfaces that are not defined by any explicit parameterization and is intended to be used in combination with level set techniques. However, contrary to the common practice with level set methods, the volume integrals derived from our formulation coincide exactly with the surface or line integrals that one wish to compute. We study various aspects of this formulation and provide a geometric interpretation of this formulation in terms of the singular values of the Jacobian matrix of the closest point mapping. Additionally, we extend the formulation - initially derived to integrate over manifolds of codimension one - to include integration along curves in three dimensions. Some numerical examples using very simple discretizations are presented to demonstrate the efficacy of the formulation.


Introduction
This paper provides simple formulations for integrating over manifolds of codimensions one, or two in R 3 , when the manifolds are described by functions that map points in R n (n = 2, 3) to their closest points on curves or surfaces using the Euclidean distance. The idea for the formulation originated in [7] where the authors proposed a formulation for computing integrals of the formˆ∂ in the level set framework, namely when the domain Ω is represented implicitly by the signed distance function to its boundary ∂Ω. Typically in a level set method [11,12,14], to evaluate an integral of the form of (1) where ∂Ω is the zero level set of a continuous function ϕ, it is necessary to extend the function v defined on the boundary ∂Ω to a neighborhood in R n . The extension of v, denotedṽ, is typically a constant extension of v. The integral is then approximated by an integral involving a regularized Dirac-δ function concentrated on ∂Ω, namelŷ ∂Ω v(x(s))ds ≈ˆR nṽ (x)δ (ϕ(x))|∇ϕ(x)|dx.
Various numerical approximations of this delta function have been proposed, see e.g. [3,2,15,16,17]. In [7], with the choice of ϕ = d ∂Ω being a signed distance function to ∂Ω, the integral (1) is expressed as an average of integrals over nearby level sets of d ∂Ω , where these nearby level sets continuously sweep a thin tubular neighborhood around the boundary ∂Ω of radius . Consequently, (1) is equivalent to the volume integral shown on the right hand side below: where δ is an averaging kernel, x * is the closest point on ∂Ω to x and J(x; d ∂Ω ) accounts for the change in curvature between the nearby level sets and the zero level set. Now suppose that ∂Ω is a smooth hypersurface in R 3 and assume that x is sufficiently close to Ω so that the closest point mapping x * = P ∂Ω (x) = argmin y∈∂Ω |x − y| is continuously differentiable. Then the restriction of P ∂Ω to ∂Ω η is a diffeormorphism between ∂Ω η and ∂Ω, where ∂Ω η := {x : d ∂Ω (x) = η}. As a result, it is possible to write integrals over ∂Ω using points on ∂Ω η as: where J(x, η) comes from the change of variable defined by P ∂Ω restricted on ∂Ω η . Averaging the above integrals respectively with a kernel, δ , compactly supported in Formula (2) then follows from the coarea formula [5] applied to the integral on the right hand side. This paper is motivated by the recent success in the closest point methods for solving partial differential equations on surfaces [8,9,10,13] as well as the need to process data sets that contain unstructured points sampled from some underlying surfaces. In the following section, we show that in three dimensions the Jacobian J in (2) is the product of the first two singular values, σ 1 and σ 2 , of the Jacobian matrix of the closest point mapping ∂P ∂Ω ∂x ; namely, To motivate the new approach using singular values, we consider Cartesian coordinate systems with the origin placed on points sufficiently close to the surface, and the z direction normal to the surface. Thus the partial derivatives of the closest point mapping in the z direction will yield zero and the partial derivative in the other two directions naturally correspond to differentiation in the tangential directions. Thus we see that one of the singular values should be 0 while the other two are related to the surface area element. We also derive a similar formula for integration along curves in three dimensions (codimension 2). The advantages of this new formula include the ease for constructing higher order approximations of J via e.g. simple differencing, even in neighborhoods of surface boundaries where curvatures become unbounded. For clarity of the exposition in the rest of the paper, we will now denote the distance function simply by d.

Integration using the closest point mapping
In this section, we relate the Jacobian J in (2) to the singular values of the Jacobian matrix of the closest point mapping from R 2 or R 3 to Γ, where Γ denotes the curves or surfaces on which integrals are defined. We assume that in three dimensions, if Γ is not closed, it has smooth boundaries.

Codimension 1
We consider a C 2 compact curve or surface Γ that can either be closed or not. If Γ is closed, then it is the boundary of a domain Ω so that Γ can be denoted ∂Ω. If Γ is not closed, we assume that it has smooth boundaries. We define d : R n → R ∪ {0} to be the distance function to Γ and P Γ to be the closest point mapping P Γ : R n → Γ (for n = 2, 3) defined as We let d 0 be the distance function to Γ if it is open and d s be the signed distance function to Γ = ∂Ω if it is closed. The signed distance function is defined as Then we define d as follows: Lemma 1. Let d be the distance function to Γ defined in (4). For |η| sufficiently close to 0, the Gaussian curvature at a point on the η level set Γ η := {ξ : d(ξ) = η} can be expressed as Proof. Starting with the definition of the Gaussian curvature G for a surface (see [6]), we can obtain an expression for the Gaussian curvature of its η-level set in terms of d as We show that this expression is the same as (5) by rearranging the terms above and using the fact that close to Γ the distance function satisfies |∇d| = 1. First we rearrange the terms in G: and rewrite each of the first six terms in terms of |∇d| 2 , e.g.
Thus we have Using (7) and rearranging the rest of the terms in (6) we obtain G = G η .
Proposition 2. Consider a C 2 compact surface Γ ⊂ R n (n = 2, 3) of codimension 1 and let d be defined as in (4). Define the closest point projection map P Γ as in (3) for x ∈ R n . For |η| sufficiently close to zero, let Γ η be the η level set of d Define the Jacobian J η as where κ η is the signed curvature of Γ η in 2D, and H η and G η are its Mean curvature and Gaussian curvature respectively in 3D.
Then if ∂P Γ ∂x is the Jacobian matrix of P Γ we have where σ 1 , σ 2 are the first two singular values of the Jacobian matrix ∂P Γ ∂x . Proof. The distance function d satisfies the property d(x) = 0 for x ∈ Γ. Also, since Γ is C 2 , its distance function d belongs to C 2 (R n , R); see e.g. [1,4]. It follows that the order of the mixed partial derivatives does not matter. In addition, the normals to a smooth interface do not focus right away so that the distance function is smooth in a tubular neighborhood T around Γ, and is linear with slope one along the normals. Therefore we have |∇d| = 1 for all x ∈ T.
The third important fact is that the Laplacian of d at a point x gives (up to a constant related to the dimension) the mean curvature of the isosurface of d passing through where H(x) is the Mean curvature of the level set {y : d(y) = d(x)}. Differentiating (10) with respect to each variable gives the following equations in three dimensions: In particular the two dimensional case can be derived by assuming that the distance function is constant in z.

Two dimensions.
In that case the Jacobian matrix ∂P Γ ∂x of the closest point projection map is Since Schwartz' Theorem holds, we have d xy = d yx making ∂P Γ ∂x a real symmetric matrix. It is therefore diagonalizable with eigenvalues 0 and 1 − d∆d. Indeed, we have Since ||v|| = 1, v is an eigenvector corresponding to the eigenvalue λ = 1−d∆d. Thus, for x such that d(x) = η we have that the eigenvalue λ of ∂P Γ ∂x satisfies λ = 1 − η∆d = 1 + ηκ η by (11). Since 1 + ηκ η ≥ 0, it follows that λ coincides with the singular value of ∂P Γ ∂x and hence Three dimensions. Since for |η| sufficiently close to 0 the distance function is C 2 , the Jacobian matrix is a real symmetric matrix which is diagonalizable with one zero eigenvalue and two other eigenvalues λ 1 and λ 2 . Indeed using (12), (13), (14) and (10) we can show that (5). Since the other two eigenvalues of ∂P Γ ∂x are the solutions of the quadratic equation where σ 1 and σ 2 are singular values of ∂P Γ ∂x . This leads to the following proposition: where δ is an averaging kernel and Σ(x)is defined as If Γ is open there is a little more to show since Equation (2) was only derived for closed manifolds. Before we state the result, it is necessary to understand how Γ η defined in (8) (an η−level set of d) looks like for an open curve in two dimensions and for a surface with boundaries in three dimensions. In two dimensions, Γ η consists of a flat tubular part on either side of the curve and two semi circles at the two ends of the curve. See Figure 1. In three dimensions Γ is in general made up of three distinct parts: the interior part, the edges of the boundary and the corners. If we assume that Γ has N edges then Figure 1: An example of an open curve Γ (black curve) and its η-level set Γ η (red curve). Γ η consists of a tubular part and two semi circles at the two ends.
where I η is the inside portion of Γ η , T η i is the cylindrical part of Γ η representing the set of points located at a distance η from the i -th edge E i , and finally S η i is the spherical part of Γ η representing the set of points located at a distance η from the i -th corner C i . See Figure 2. Figure 2: An example of a surface with boundaries viewed from different angles (2(a), 2(b) and 2(c)) and its corresponding η-level set Γ η viewed from the same angles (2(d), 2(e) and 2(f)). Figure 2(f) shows the surface and Γ η .
In both cases we need to integrate over Γ η and then subtract the two semi circles at the two end points of the curve (in two dimensions) or subtract the portions of sphere at the corners of the surface and the portions of cylinders at the edges of the surface (in three dimensions). However, it turns out that the subtraction is unnecessary since Σ(x) = 0 on each of the subtracted pieces as shown below.
• Two dimensions. On the semi-circle around the end point of a curve, the closest point mapping is constant since all points on the semi-circle Γ η map to the end point. As a result, the singular values of the Jacobian matrix of the closest point mapping are all zeros and thus Σ(x) = 0 on the semi-circles around the end points of a curve.
• Three dimensions. As in two dimensions, on the portions of sphere around a corner point of a surface, the closest point mapping is constant and thus Σ(x) = 0. On the portion of cylinders, the closest point mapping is constant along the radial dimension (one of the principal directions or singular vector) resulting of the singular value along that direction to be zero. Since Σ(x) is the product of the singular values, it follows that Σ(x) = 0 on the portion of cylinders as well.
Consequently, Equation (15) holds for any C 2 curve or surface with C 2 boundaries of codimension 1.

Codimension 2
We consider a C 2 curve in R 3 denoted by Γ and let γ(s) be a parameterization by arclength of Γ. We denote by d : R 3 → R + ∪ {0} the distance function to Γ and let P Γ : R 3 → Γ be the closest point mapping to Γ. We consider a parameterization of the tubular part of the level surface for η ∈ [0, ] defined as where T = dγ ds , N and B constitute the Frenet frame for γ as illustrated in Figure 3. As in the previous section, if Γ is closed then d is the signed distance function to Γ . If we project a point x on the tubular part of the level surface Γ η defined in (8), we have P Γ (x(s, θ, η)) = γ(s). If L is the length of the curve it follows that 2π 0ˆL 0 g(P Γ (x(s, θ, η)))|x s × x θ |dsdθ =ˆ2 π 0ˆL 0 g(γ(s))η(1 − ηκ(s) cos θ)dsdθ, = ηˆL 0 g(γ(s))ˆ2 π 0 (1 − ηκ cos θ)dθds, = 2πηˆg(γ(s))ds.
Note that the tubular part of the level surface Γ η does not contain the two hemispheres of Γ η which are located at the two end points of the curve Γ. Thus, where C 1 and C 2 are the two hemispheres of the level surface Γ η located at the two end points of the curve Γ. Consequently, for sufficiently small and by the coarea formula we obtain Γ g(γ(s))ds = 1 2πˆ 0 where K is a C 1 averaging kernel supported in [0, ] and χ (C 1 ∪C 2 ) c (x) is the characteristic function of the set (C 1 ∪ C 2 ) c . In our numerical simulations we consider the kernel Since the formulation above does not use the two hemispheres located at both end points of the curve, in order to integrate over the tubular part of Γ η only, it is necessary to subtract the integration over each of the hemispheres C 1 and C 2 . The result can be summarized in the following proposition: Proposition 4. Consider a single C 2 curve Γ in R 3 parameterized by γ(s) where s is the arclength parameter, and let d be the distance function to Γ. We define K to be a C 1 averaging kernel compactly supported in [0, ] and P Γ : R 3 → Γ to be the closest point mapping to Γ.
If g is a continuous function defined on Γ then for sufficiently small > 0 we havê Γ g(γ(s))ds = 1 2πˆR3 where x η is a point on a sphere of radius η.
Note that for the computation of the length of a curve, the correction terms given by integrating over both C 1 and C 2 iŝ This simple correction is, however, not suitable for more general cases that contain multiple curve segments and several integrands. We shall derive a more elegant and seamless way to perform such correction in the following section. Now if we consider a C 2 curve in three dimensions and let P Γ be its closest point mapping, we have the following proposition: Theorem 5. Let σ(x) be the nonzero singular value of ∂P Γ ∂x and let g be a continuous function defined on Γ. If γ(s) is the arclength parameterization of Γ and if is sufficiently small we havê Γ g(γ(s))ds = 1 2πˆR3 where d is the distance function to Γ.
Case 1: x is on the spherical part of Γ η corresponding to the η-distance to either of the two end points of the curve Γ. WLOG we assume that x is at a distance η from the first end point C 1 parameterized by γ(0). The result is the same if x is on the other sphere, i.e. at a distance η from the other end point C 2 . In that case, P Γ (x) = γ(0) for all x on the spherical part so that the Jacobian matrix ∂P Γ (x) ∂x = 0 . Therefore, for x on the spherical part of Γ η , all singular values of the Jacobian matrix are zero.
Case 2: x is on the tubular part of Γ η . In that case, if we use the Frenet frame centered at the point x = x(s, θ, η) ∈ Γ η , we can write x in the new coordinate system ( T, N, B) as where u = 0 is the coordinate of x along T, v is the coordinate along N and w is the coordinate along B. Since the projection P Γ (x) = γ(s) does not depend on v nor w (since the plane ( N, B) is normal to the curve Γ) is follows that On the other hand, we have where ∂s ∂u is the variation of the arclength parameter s with respect to u when the point x is moving on Γ η along the tangential direction T. Since u is the arclength parameter along the tangential direction T, it follows that we have a unit speed parameterization along T giving the identity ∂x ∂u · T = 1.
In addition, where κ is the curvature of Γ at γ(s) and τ is the torsion of the curve Γ at the point γ(s). Since the level surface Γ η is a tube of radius η, its intersection with the normal plane ( N, B) is a circle of radius η. Hence if we use polar coordinates on the normal plane, we obtain v = η cos θ and w = η sin θ . It follows that ∂x ∂s · T = 1 − κη cos θ. Therefore, in the Frenet frame, the Jacobian matrix of the closest point projection map can be written as where 1 1−κη cos θ is the nonzero eigenvalue of the Jacobian of the closest point mapping. Since 1 1−κη cos θ ≥ 0, it is also its singular value σ(x). Therefore we have Now using (16) and (17) we obtain It follows that for K a C 1 averaging kernel compactly supported in [0, ], for sufficiently small and by the coarea formula, we havê

Numerical simulations
In this section we investigate the convergence of our numerical integration using simple Riemann sum over uniform Cartesian grids. In our computations we use the cosine kernel K cos (η) = χ [0, ] (η) 1 2 1 + cos πη for integration on surfaces of codimension 1, and the kernel K 1,1 defined in (18) for codimension 2. Unless stated otherwise, the singular values are computed from the matrix whose elements are computed by the standard central difference approximations of the Jacobian matrix ∂P Γ ∂x . With these compactly supported kernels, formulas (15) and (20) can be considered integration of functions defined on suitable hypercubes, periodically extended. In such settings, simple Riemann sums on Cartesian grids are equivalent to sums suing Trapezoidal rule, and the order of accuracy will be related in general to the smoothness of the integrands; exception can be found when the normals of the surfaces are rationally dependent on the step sizes used in the Cartesian grids.

Integration of codimension one surfaces
We tested our numerical integration on two different portions of circle, a torus, a quarter sphere and a three quarter sphere. We computed their respective lengths or surface areas by integrating the constant 1 over the curve or surface. Each of these tests were designed to exhibit the convergence rate of our formulations on cases with varying difficulty. In particular, the convergence rate of our formulation depends on the smoothness of the closest point mapping inside the tubular neighborhood of the surfaces.
The results for the portions of circle are given in Tables 1 and 2. In the first convergence studies (Table 1) the line where the closest point mapping has a jump discontinuity is parallel to the grid lines. In this case we see a second order convergence rate using central differencing to compute the Jacobian matrix ∂P Γ ∂x . In the second test case however, the portion of circle is chosen so that the line where the closest point mapping has a jump discontinuity is not parallel to the grid lines. In that case the normal to the curve is rationally dependent on the step size of the Cartesian grid and the convergence rate reduces to first order even though we used central differencing to compute ∂P Γ ∂x . We note that in these two tests, we chose (the half width of the tubular neighborhood around the curve) small enough so that the line where the closest point mapping is discontinuous is outside of it.
In three dimensions we first tested our method on a torus (closed smooth surface). The results for the torus are reported in Table 3. In this case the closest point mapping is very smooth and we see third order convergence when using a third order difference scheme to approximate ∂P Γ ∂x . For surfaces with boundaries we tested the method on a quarter sphere and a three quarter sphere. The three quarter sphere case is illustrated in Figure 4. Figure 4: The three quarter sphere (4(a)) and its corresponding η-level set Γ η (4(b)).
The reason for choosing these two cases is because the closest point mapping has a different degree of smoothness for each of these surfaces. For the quarter sphere the closest point mapping is smooth enough, but for the three quarter sphere, the tubular neighborhood around the surface contains the line where the closest point mapping has a jump discontinuity. In that latter case, it is therefore necessary to use an adequate one sided discretization to compute ∂P Γ ∂x accurately. The one-sided discretization that we used is reported in Section 3.3. The test for the quarter sphere still uses central differencing to compute ∂P Γ ∂x . The results for the portions of sphere are reported in Tables 4 and 5. Table 1: Relative errors in the numerical approximation of the length of a planar curve, which is a portion of circle of radius R = 0.75 centered at 0. The width for the tubular neighborhood of the curve is = 0.2. In this computation, the closest point mapping has a jump discontinuity along a straight-line which is arranged to be parallel to the grid lines.     Table 5: Relative errors in the numerical approximation of the surface area of a 3/4 sphere with radius R = 0.75 centered at 0 (this is the portion of a sphere that misses half of a hemisphere). In this computation, we summed up grid points that are within = 0.2 distance from the surface. Due to this setup, the closest point mapping has a discontinuity that stems out from the boundary of the surface. We used the discretization described in Section 3.3 to compute each entry of the Jacobian matrix

Integrating along curves in three dimensions
In codimension 2, we tested our numerical integration on a coil wrapped around the helix defined parametrically as x(t) = (r cos(t), r sin(t), bt) , with r = 0.75 and b = 0.25. The coil is then wrapped around the helix at a distance of 0.2 from the helix. As our test case, we computed the length of the coil by integrating 1 along the curve. The results are reported in Table 6.

One-sided discretization of the Jacobian matrix
Here for completeness, we describe the one-sided discretization used in computing results reported in Table 5. For simplicity, we will describe the one-sided discretization for a uniform Cartesian grid, namely for P Γ (x i,j ) = (U i,j , V i,j ) with x i,j = (ih, jh), i, j ∈ Z and h > 0 being the step size. The Jacobian matrix will be approximated by simple finite differences defined below: The discretization of U and V have to be defined together because the two functions are not independent of each other. With and the smoothness indicator x ) i,j , otherwise, and (V x ) i,j is defined according to the choice of stencil based on S ± (U i,j ) The discretization of U y and V y is defined similarly with the choice of the stencil based on S ± (V i,j ).

Summary
In this paper, we presented a new approach for computing integrals along curves and surfaces that are defined either implicitly by the distance function to these manifolds or by the closest point mappings. We are motivated by the abundance of discrete point sets sampled from surfaces using devices such as LIDAR, the need to compute functionals defined over the underlying surfaces, as well as many applications involving the level set method or the use of closest point methods. Contrary to most other existing approximations using either smeared out Dirac delta functions or locally obtained parameterized patches, we derive a volume integral in the embedding Euclidean space which is equivalent to the desired surface or line integrals. This allows for easy construction of higher order numerical approximations of these integrals. The key components of this new approach include the use of singular values of the Jacobian matrix of the closest point mapping, which can be computed easily to high order even by simple finite differences.