A local Mesh Method for Solving PDEs on Point Clouds

In this work, we introduce a numerical method to approximate diﬀerential operators and integrals on point clouds sampled from a two dimensional manifold embedded in R n . Global mesh structure is usually hard to construct in this case. While our method only relies on the local mesh structure at each data point, which is constructed through local triangulation in the tangent space obtained by local principal component analysis (PCA). Once the local mesh is available, we propose numerical schemes to approximate diﬀerential operators and deﬁne mass matrix and stiﬀness matrix on point clouds, which can be used to solve partial diﬀerential equation (PDE) and variational problem on point clouds. As numerical examples, we use the local mesh method and variational formulation to solve the Laplace-Beltrami eigenproblem and solve the Eikonal equation to compute distance map and geodesics on point clouds.


Introduction
In many problems in science and engineering, data is commonly represented as a collection of points, referred as a point cloud, embedding in a certain space. In practice, a common feature for many point cloud data is that those data points sit on a low dimensional manifold M in an ambient space R n possible with a large dimension. Analyzing, processing and characterizing point cloud become important tasks in many applications, such as computer visions, data mining and machine learning [1,2,3].
A typical example is 3D modeling in computer graphics and computer vision, where geometric objects, usually 2-dimensional surfaces, are represented as a set of 3D points obtained by a laser scanner. In practice, a triangulated surface or an implicit representation, such as level set function, is constructed to approximate the point cloud data. Based on these representations, partial differential equation (PDE) based methods and variational methods provide powerful tools to extract intrinsic geometric information either locally or globally for the underlying manifolds. Here, we propose a local mesh method that can be applied to solving PDE and variational problem on point clouds without requiring the construction of a global mesh or representation of the point cloud, which may be difficult and computationally expensive in three and higher dimensions.
To solve variational problems or differential equations on point clouds, we first need to approximate integrals and differential operators. Previously, Belkin et al. [4] proposed to approximate the Laplace-Beltrami operator on point clouds through the heat diffusion in the ambient Euclidean space, and then the Laplace-Beltrami eigenproblem (Helmholtz equation) on point clouds is discussed. Construction of a local mesh in the tangent space is needed to approximate integration of the heat kernel. Their method suffers from the low order approximation and is only limited to the approximation of the Laplace-Beltrami operator on point clouds. After that, Luo et al. [5] proposed to approximate integrals on point clouds through computing summation with voronoi weight or principle eigenvector weight. More recently, Liang et al. [6,7] proposed a more systematic way to approximate differential operators on point clouds. In their method, discrete approximation of differential operators are constructed by local least square approximations of the manifold and hence the metric using K nearest neighbors (KNN). Since the differential operators are intrinsically approximated from more accurate local manifold reconstruction, their method can achieve high order accuracy and enjoy more flexibility since no mesh is needed at all. It can be applied to manifolds with arbitrary dimensions and co-dimensions with or without boundary. Moreover, the computational complexity depends on the true dimension of the manifold rather than the dimension of the embedded space. However, their method can not deal with integrals and variational problems easily due to the lack of meshes.
In this paper, we propose an approach that requires only local mesh in the tangent space and can approximate both differential operators and integrals on point clouds. Hence, this method applies to cases where the tangent space can be locally triangulated. Our method is based on several simple observations as follows. 1 . Local mesh construction is much easier than global mesh construction for general point clouds; 2 . Definition of differential operators at any point only depends on the local structure of the point cloud; 3 . Integration can be represented as an inner product weighted by the mass matrix, which can be obtained from local information of point clouds by choosing a locally supported basis. To clearly represent our approach in this paper, we only restrict our discussion on point clouds sampled from two dimensional manifolds in an ambient space R n , while all techniques we discussed here can be easily adapted to problems on point clouds sampled from k−dimensional manifolds in R n . We first use local PCA to estimate the tangent space at each point, then the local mesh can be obtained from the Delaunay triangulation in the tangent space.
After the local mesh is obtained, differential operators can be approximated. With the constructed local connectivity, we propose a simple numerical approximation of the mass matrix and stiffness matrix similar to those using finite elements on triangulated surfaces. Then integrals, such as area and volume, and variational problem can be approximated numerically. Based on the variation formulation and our local mesh method, we solve the Laplace-Beltrami eigen-problem. Moreover, we also numerically solve the Eikonal equation to compute distance map and geodesics on point clouds based on either the fast sweeping or the fast marching method.
The rest of the paper is organized as follows. In section 2, we introduce the idea of local mesh construction for point clouds. Then differential operators and integrals can be approximate for given point clouds based on the constructed local mesh. Selected applications to the local mesh method are discussed section 3. We propose an analogue of finite element method to solve Laplace-Beltrami eigenproblems on point clouds. Moreover, we discuss numerical approximation of distance maps and geodesics on point clouds by solving the Eikonal equation based on the fast marching or the fast sweeping methods. After that, section 4 reports all numerical experiments and comparisons to demonstrate the proposed methods. Conclusions are made in section 5.

Local Mesh method
Given a point cloud P = {p i ∈ R n |i = 1, . . . , N } sampled from a smooth surface/manifold M, it is a challenge to extract its global information without a global mesh or parametrization. For instance, a elementary geometric quantity such as the area of M is not easy to obtain. One way is to first approximate the manifold M by constructing a global triangulation or implicit representation. Then, many differential geometry, PDE and variational tools can be applied to M. However, surface reconstruction for point clouds may be difficult and time consuming, especially for a large amount of data in three or higher dimensions with complicated geometry and topology plus possible noise and nonuniform sampling. On the other hand, local geometry and connectivity at each point can be easily obtained. This motivates our method for approximating differential operator, PDEs and variational formulation on point clouds based on local mesh. Moreover, one can obtain useful global information for point clouds by solving PDE problems on point cloud [6].

Local connectivity construction for point clouds
Let's denote the indices set of K-nearest neighborhood (KNN) of each point p i ∈ P by I(i). Local PCA on KNN is well-known and widely used for local linear approximation for point clouds [8].
One can determine the local tangent space and normal space using local PCA. For convenience, we We can estimate local linear structure near p i using the covariance matrix P i of N (i), defined by: Here, c i is the local barycenter c i = 1 K k∈N (i) p k . The eigenvectors (e 1 i , e 2 i , e 3 i ) of P i form an orthogonal frame associated with If the point cloud is sampled from a two dimensional manifold, and the local sampling is dense enough to resolve local feature size, then i . e 1 i ≥ e 2 i provide two orthogonal directions in the tangent plane, and e 3 i is the normal direction of P at the point p i .
To construct the local connectivity and a mesh at the point p i , we project its K-nearest neighborhood N (i) on the tangent plane T i of p i spanned by e 1 i , e 2 i . Namely, we have the following construction:p

Numerical approximation of differential operators on point clouds
Since differential operators, such as gradient, divergence and Laplacian, are defined locally, we can obtain point-wise numerical approximation of these differential operators using the local mesh constructed in section 2.1. Here, we directly borrow the idea of approximating differential operators on triangle meshes [9,10,11,12,13] on the first ring of each point of the point cloud P.
To make the paper self-contained, we write down details of the numerical approximation of gradient, divergence and Laplace operators on the local mesh at a point p i . First, the discretization is defined on on a single triangle by their coordinate invariant definition in differential geometry [14,15]. Then an average in the first ring is taken weighted by the area of each triangle.
First of all, we show the discretization on a single triangle T with three vertices {p i ∈ R n | i = 0, 1, 2}. In the discrete case, we have a function f = {f (p 0 ), f (p 1 ), f (p 2 )} and a vector field } defined on each vertex respectively. With the barycentric coordinates Then we have ∂ u 1 = p 1 − p 0 , ∂ u 2 = p 2 − p 0 , and the metric matrix of T would be: where · is the dot product in R n . We then have the following discretization: We next discretize the divergence operator on the triangle T . Suppose Differentiate both sides of above equality, we have: Since √ G is constant on each triangle, we can obtain the discretization of the divergence operator on triangle T as follows: Now, we consider the discretization of the gradient and divergence operators on the local mesh at a point. We take a weighted average on the first ring of p i in terms of triangle areas.
Namely, for any function f and vector field − → V defined on the given point cloud P, the gradient and divergence operators at each point are approximated as follows: Moreover, we can also approximate the Laplace-Beltrami operator as follows: where , α 1 ij and α 2 ij are the two angles opposite to the edge p i p j , V(i) is the first ring neighborhood of the vertex p i .

Remark 1
Here, we only consider approximate the differential operators using the first ring on local mesh obtained from section 2.1. More accurate approximation can be obtained using the second ring structure.

Integral approximation
Integration on point clouds, such approximation of surface area, enclosed volume, and variational formula, is an important task in many applications. One way is to reconstruct the surface first, which can be difficult and expensive itself. More direct methods have been proposed before. Li et al. [16] propose to approximate areas of surfaces in R 3 by estimating the intersections of lines in R 3 based on Cauchy-Crofton formula. More recently, Lui et al. [5] propose to estimate area of point clouds in R n using voronoi cells of local mesh reconstruction from point clouds. Here we propose a method for integral approximation using local mesh.
First we construct the mass matrix on a point cloud P through which we can approximate integration on P as follows. Given a point cloud P = {p i ∈ R n | i = 1, · · · N } with local connectivity Then, any function Consider the linear interpolation of e i on its first ring structure, we define the inner product of e i , e j on P: We define the N × N mass matrix M = (m ij ) for the point cloud P by: Here M is a sparse and symmetric matrix. Let's denote the constant one function on P by 1 = (1, · · · , 1) T . For any two functions f = (f (p 1 ), · · · , f (p N )) T and g = (g(p 1 ), · · · , g(p N )) T , we can approximate the following integrations on the point cloud P by: For example, we can approximate the area of the underlying surface of P by: If the underlying surface of the given point cloud P is a closed surface M in R 3 , we can also estimate the 3D volume of the manifold D bounded by M using Stokes' theorem. Assume each We define three global coordinates functions of P by: The Stokes' theorem tells us: where ν = (ν 1 , ν 2 , ν 3 ) is the outward normal vector of M. Therefore, we approximate the volume of D as follows: 3 Solving PDEs on point clouds As long as differential operators and integrals can be approximated on the given point clouds based on the local mesh method, many problems related to differential equations on point cloud can be further studied. As applications, we propose an analogue of the finite element method on point cloud to solve a typical elliptic eigenvalue problem, Laplace-Beltrami eigenvalue problem. After that, we demonstrate how to use the proposed local mesh method to solve the Eikonal equation, a typical nonlinear partial differential equation, on the given point cloud by adapting the fast marching and fast sweeping methods [17,18,19,20]. Using the resulting distance maps, we also discuss the construction of geodesics on point clouds.

Laplace-Beltrami eigenproblems
As an intrinsic differential operator defined on a surface, the Laplace-Beltrami (LB) operator can be used as a bridge to connect local and global geometry of the surface. Its eigenvalues and eigenfunctions can be used to study and characterize the surface. For example, generic LB eigenfunctions of surfaces are morse functions [21], which enable LB eigenfunctions as descriptors for the topologies of the surfaces. Moreover, LB eigenfunctions can be viewed as either global or local descriptors to analyze surface geometric structures [22,23,24,25,26]. Furthermore, the LB operator is also closely related to harmonic maps between two surfaces. Recently a few numerical methods [4,6,7] have been proposed to compute LB eigenproblems on point clouds. These methods are all based on pointwise discretization of the LB operator. Here, we design a variational formulation using the local mesh method, which can be used to solve other type of elliptic PDEs on point clouds as well.
Given a point cloud P = {p i ∈ R n | i = 1, · · · N } sampled from a Riemannian surface M embedded in R n , a discrete analogue of LB eigen-system on P can also be studied. Based on our local mesh method, we can use the numerical approximation of the LB operator on P given by (11) and compute the spectrum of the corresponding matrix. We can also compute LB eigensystem on point clouds based on the variational formulation in (21) Consider {e i } N i=1 as the set of linear elements defined on P given in (12), and write E = Then the discrete version of LB eigensystem on P can be defined by the following weak formulation: In fact, if we can write φ = N i v i e i , the essential part of the above problem is a numerical approximation of P ∇ P e i · ∇ P e j , which forms the stiffness matrix S = (S ij ) N ×N for the given point cloud P. Ideally, we would like S to have properties, such as symmetric and nonnegative definite, same as those for the stiffness matrix for a triangulated surface [29]. However, these are global properties which may not be possible to achieve due to the use of local mesh only. To be more 1 If M is a closed surface, then ∂M = Ø specific, the first ring structure of p i is not necessary compatible with the first ring structure of p j although p j belongs to the first ring of p i . To have a numerical approximation of the stiffness matrix, we first define where ∇ T e i and ∇ T e j are computed by (5) and α k ij (p i ), k = 1, 2 are the angles opposite to the edge connecting points p i and p j in the first ring of p i . This is the approximation of P ∇ P e i · ∇ P e j based on the first ring structure at p i . Note that A ij may not be equal to A ji due to the possible incompatibility of the first ring structures of p i and p j . One simple symmetrized definition of the stiffness matrix is the following: The above definition of the diagonal element is to enforce the consistence condition, i.e., constant function is an eigenfunction with zero eigenvalue. In particular, if all triangles in the first ring structure are acute, off-diagonal elements are non-positive and diagonal elements S ii = | k =i S ik | are positive. Hence all eigenvalues are real and non-negative. When the distribution of points is reasonably uniform, this definition of stiffness matrix works quite well. However, when the points are distributed non-uniformly, the first ring structure of p i is more likely incompatible with the first ring structure of p j . Numerical tests suggest that the following definition for off-diagonals elements produces better results (see Figure. 5): The intuition for this definition of off-diagonals is that (a) no direct coupling between p i and p j unless they are both in each other's first ring, (b) use better triangles in the first ring structures at p i and p j to approximate P ∇ P e i · ∇ P e j when they are both in each other's first rings. The first case in (25) means we use triangles with acute angles α k ij or α k ji further away from 0 degree. The second case in (25) means we use triangles with obtuse angles α k ij or α k ji further away from 180 degree. The third case in (25) means we use triangles with acute angles α k ij or α k ji rather than those with obtuse angles.
Then the solutions to (22) can be computed from the following generalized matrix eigen-problem: Note that the stiffness matrix S and the mass matrix M are symmetric and sparse matrices of size N × N . Thus, all eigenvalues and eigenfunctions computed from the above method are real.
Computationally, there are a variety of softwares, such as Matlab, to solve the above problem (26).
We would like to point out that this approach is not necessary restricted to solve LB eigensystem.
Similar approach can be used to solve other types of PDE, such as the poisson equation, the heat equation etc., on point clouds.

Distance maps and geodesics on point clouds
Once local connectivity is available on a point cloud, one can use appropriate local solver on triangular mesh to solve many types of PDEs on the point cloud. Here we specifically demonstrate how to solve a special type of nonlinear hyperbolic PDE, the Eikonal equation (27), to compute distance maps and geodesics on point clouds using local mesh method. Here is the Eikonal equation for the distance map to a given set Γ on a manifold M: In particularly, a distance map d p to a given point p can be computed by setting Γ = p. Once the distance map d p is obtained, one can construct the geodesic from any point q ∈ M to p by solving the ODE: where X(s) traces out the geodesic path between two points p, q ∈ M and d p (·) is the distance function with boundary condition d p (p) = 0.
To solve the Eikonal equation (27) on a point cloud P, the first step is to design an appropriate discretization at each point of P. Since we have a local mesh, i.e., the first ring structure, we can directly adopt the monotone upwind discretization for triangular mesh [18,20], which is equivalent to an approximation of the dynamic programming principle on the first ring. As the result of the discretization, at each point p i ∈ P we have a nonlinear equation that couples the value at p i with the values of its first ring neighbors. Hence we have to solve a system of nonlinear equations numerically. Fortunately, we have efficient ways to solve it either using the fast marching method [17] or the fast sweeping method [19] which are also generalized to triangular mesh in [18] and [20] respectively, where only first ring structures of the given triangular mesh are needed. Thus, the proposed local mesh method can be directly adapted to either fast marching method or fast sweeping method to solve the Eikonal equation on point clouds. This is a more direct and more accurate approach to obtain distance maps than the previous method proposed in [45], where the Eikonal equation is solved in the tubular neighborhood of the point cloud P in R n to approximate the true Eikonal equation defined on P.
Once we have computed the distance map d p (x) on the point cloud P, we can find the shortest path from p to q on P by solving (28). To have more accurate approximation of geodesics and avoid zigzag paths, we can not require the nodes of discretized geodesic path strictly passing through points of P. The question is how to add new points in the way that the geodesic can truly reflect found in [6] and [7]). We construct the Delaunay triangulation of the projectionsp c ,p 1 ,p 2 , · · · ,p K and find the first ring R = {T 1 c , · · · , T l c } of p c , which is the same as we did in Section 2.1. Suppose v has local coordinate (v 1 , v 2 , v 3 ) in p c ; e 1 , e 2 , e 3 . We find the intersection of line segment starting at p c with the direction (−v 1 , −v 2 , 0) and the first ring, notice that this computation is done within the tangent space of p c . Denote the intersection as p 0 = (x 0 , y 0 , 0), we then project it back to the approximated surface to obtain the next point on the geodesic path p n = (x 0 , y 0 , z(x 0 , y 0 )). This process is illustrated in Figure 2.
We start from the target point and set it as the current point on the geodesic path. Using the algorithm introduced above, we trace back the distance field to get the next point on the geodesic path (notice that it may not belong to the point cloud P) and set it to be the current one. We repeat the process until the current point is close enough to the source point. Numerical examples will be shown in Section 4.

Numerical Results
In this section, numerical experiments are presented to illustrate applications of the proposed local mesh method to point clouds. First, we report numerical results of the area and volume approximation for point clouds sampled from different surfaces. Second, we demonstrate the proposed method in section 3.1 to solve the Laplace-Beltrami eigenproblem for point cloud data. In addi-

Area and Volume approximation
In general accuracy of the area and volume approximation depends on the size and distribution of the point cloud sampled from the given closed surface. Larger size and better distribution sampling lead to better numerical approximation. To demonstrate this point, we first test our method on different point clouds sampled from the unit sphere S 2 in R 3 . We choose two groups of point clouds on S 2 . One group of point clouds are uniformly sampled from S 2 with the number of points about geometrically increasing from 1k to 16k, and another group of point clouds are non-uniformly sampled with the same number of points as the first group. The first row of Figure 3 plots two uniformly sampled and two non-uniformly sampled point clouds used in this experiment.
As a comparison, we also compute the area approximation using the voronoi weighting scheme introduced in [5]. The numerical results and the corresponding relative errors for this two groups of data are reported in Table. 4.1. According to our results, we can observe that our method provides highly accurate approximation of the area and volume for the unit sphere in the case of uniformly sampled data set, while the approximation results using nonuniformly sampled data are not as good as the ones obtained from uniformly sampled data set. In both cases, our results are more accurate than the results using the voronoi weighting scheme in [5]. In addition, we can also observe that the relative errors halve as the number of point doubles for uniformly sampled points. However, this behavior can not be observed for nonuniformly sampled points. Similar phenomenons can be observed for point clouds sampled from a flat torus T = {(cos u, sin u, Next apply our method to a few real data sets illustrated in the second row of Figure 3. Since there is no true area and volume for these surfaces, we use the results computed from their triangular mesh representation as the ground truth. The numerical approximation and the relative errors are reported in the last four rows of

Laplace-Beltrami eigenproblems
To test the accuracy of the proposed method for solving LB eigenproblems of point clouds, we would like to first apply our method to point clouds data sampled from the unit sphere S 2 . For the unit sphere, its exact value of the n-th LB eigenvalue is given by λ n = n(n + 1), with multiplicity 2n + 1. Thus, we can compare the numerical results with these true values. In other words, we compute the relative error: |λ n,i − λ n | λ n whereλ n,i 's are the approximated LB eigenvalues to λ n computed from the proposed method introduced in section 3.1 , and i runs over λ n 's multiplicity. By the definition, E max,n represents the worst possible error in computing λ n .
In our experiments, we apply the two definitions of the stiffness matrix in (24) and (25) to the same two groups of spherical point clouds used in section 4.1. As comparisons, we also compute the LB eigenproblem using other point cloud based methods, such as the MLS (using quadratic approximation) method [6], the PCD method [4], and a mesh based finite element method [11] for the same two groups of data. We show E max,n for λ = 20 and 72 of the uniformly and non-uniformly sampled data in Figure 4 and in Figure 5, respectively. For the uniformly sampled data set, the numerical errors of our local mesh method using (24) and (25) are the same up to 10 −10 , thus we only list one group of results using the local mesh method. In this case, it can be seen that the error of our approach is of the same order as the mesh based method, while they are slightly better that the results obtained by MLS approach and much better than PCD Laplacian method.
The results for nonuniformly sampled point clouds are not as accurate as results obtained from the uniformly sampled data, while they are also comparable with the MLS approach and more accurate than PCD Laplacian method. In addition, it is also clear to see that the results using the stiffness matrix defined in (25) are more accurate than those using (24) for the non-uniform sampled data.
Our next experiment reports numerical results of geodesic tracing on point clouds using the method discussed in section 3.2. We choose point clouds sampled from the unit sphere S 2 ⊂ R 3 .
They can also be viewed as point clouds in R n by simply adding the last n − 3 coordinates as zero.
After that, by applying a randomly chosen orthonormal n × n matrix as a rotation and a randomly chosen vector as a translation, we can obtain point clouds sampled from the isometrically embedded S 2 in R n . We apply our developed method to compute distance maps and geodesics on these point clouds to demonstrate the capability in handling point clouds in high dimension (n=10). Table. 2 shows our results of distance between the north pole and south pole. It is clear to see that our method provides highly accurate estimation for geodesics on point clouds sampled from the unit sphere S 2 embedded in R 10 . In addition, we also report numerical results of geodesic tracing on point clouds sampled from more complicated surfaces in the first row of Figure

Conclusion
In this work, we propose a local mesh method that can be used to approximate differential operators, integrals, variational formulation and to solve PDEs on point clouds. The local connectivity and mesh at each point is constructed through local approximation of the tangent space, which is much easier and cheaper than the construction of a global mesh, especially in high dimensions. Examples of the area and volume approximation, solving eigenvalue problem of Laplace-Beltrami operator using variational formulation, and computing distance maps and geodesics on point clouds are presented.