 # A Cylindrical Basis Function for Solving Partial Differential Equations on Manifolds

Numerical solutions of partial differential equations (PDEs) on manifolds continues to generate a lot of interest among scientists in the natural and applied sciences. On the other hand, recent developments of 3D scanning and computer vision technologies have produced a large number of 3D surface models represented as point clouds. Herein, we develop a simple and efficient method for solving PDEs on closed surfaces represented as point clouds. By projecting the radial vector of standard radial basis function(RBF) kernels onto the local tangent plane, we are able to produce a representation of functions that permits the replacement of surface differential operators with their Cartesian equivalent. We demonstrate, numerically, the efficiency of the method in discretizing the Laplace Beltrami operator.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 0.1 Introduction

Many applications in the natural and applied sciences require the solution of partial differential equations on manifolds. Such applications arise in areas such as computer graphics[18, 10, 1], image processing[20, 2, 9, 21], mathematical physics, biological systems[17, 8], and fluid dynamics[11, 12, 19].A lot more interest, especially in computer graphics, has been generated around solution of PDEs on closed 2D manifolds, as these arise as boundaries of 3D objects. The development of high resolution 3D scanning devices, which capture these surfaces as point clouds, have made numerical methods that can be applied directly to point clouds very attractive.

A class of numerical methods that have been developed to solve PDEs on closed surfaces, involve the expression of differential operators as projections of their Cartesian equivalents onto local tangent planes via a projection operator (). The resulting operators are then discretized using, for example, finite element methods. Recently, the projection method has been extended to PDE’s defined on manifolds represented as point clouds . For such methods, functions defined on the manifold are represented using radial basis functions (RBF). The surface differential operators are obtained by applying a projection operator to the RBF discretization of the Cartesian equivalent.

While the projection methods modify Cartesian differential operators, another class of methods embed the surface PDE into so that solutions to the embedded problem when restricted to the surface provide the solution on the surface. Because these methods result in embeded PDEs posed in , spatial complications arise when the PDEs have to be solved on a restricted surface domain. The closest point method, developed in , attempts to resolve this problem by extending the functions into

in a way that makes them constant in the normal direction. This allows the simple replacement of surface differential operator by their Cartesian equivalent. The method however requires a high order interpolation at each time step in order to obtain solutions on the surface.

Recently, the orthogonal gradient method presented in, was introduced to extend this idea to point clouds. Here, additional nodes are introduced during the construction of a distance function, to force the function defined on the surface to be constant in the normal direction. The nodes are chosen according to an offset parameter which controls their distance from the surface. The accuracy of the method and condition number of the resulting differential matrix is however sensitive to the choice of . The use of additional nodes, to enforce derivative constraints, increases the computational complexity of forming the interpolation and differential matrices.

In this work we propose a modified RBF kernel, the Cylindrical Basis Function (CBF), which is intrinsically constant in the normal direction at each point of a surface. The modified kernel, when used in the representation of smooth functions defined on a closed surface, allows surface differential operators to be replaced by their Cartesian equivalent without the need to impose additional constraints on the function or have an implicit representation of the surface. We also avoid inherent challenges in performing higher order interpolation, typical of embedding techniques, by discretizing operators directly on the manifold. The proposed method is simple and requires the solution of a much smaller linear system, compared to the orthogonal gradient method.

## 0.2 The Radial Basis Function

Given function data at the node locations , the RBF interpolant to the data is given as

 s(x)=N∑i=1λiϕ(∥x−xi∥), (1)

where is the RBF kernel centered at the node , and are coefficients chosen to satisfy the interpolation conditions

 s(xi)=f(xi)i=1⋯N, (2)

which is equivalent to solving the linear system

 A→λ=→f. (3)

is the matrix with entries , usually referred to as the interpolation matrix. The norm is taken to be the Euclidean norm. Some common kernels are presented in Table 1.

## 0.3 RBF Discretization of Cartesian Differential Operators

Given a point cloud and data from a smooth function defined on these points, an RBF interpolation of satisfies,

 f(xi)=N∑j=1λjϕ(||xi−xj||)i=1,2,⋯,N. (4)

Let be a differential operator acting on and the value of at the point then,

 (5)

which defines a linear system and can be represented in matrix form as,

 B→λ=→g, (6)

where the differential matrix has entries .

The interpolation matrix A in Eq.(3) is non-singular (see ) and permits the substitution of into Eq.(6) leading to . The differentiation matrix gives the RBF discretization of L with respect to the point cloud P.

## 0.4 Construction of Surface Differential Operators

Let be a smooth function defined on an arbitrary surface . The gradient of expanded in the normal, first and second tangent orthogonal coordinates at some point can be expressed as

 ∇f=∂nf→n+∂t1f→t1+∂t2f→t2. (7)

The surface gradient operator is the projection of the regular gradient onto the local tangent plane at . Thus,

 ∇Γf=∂t1f→t1+∂t2f→t2. (8)

Two approaches are typically used in the literature [4, 13] to obtain surface operators from regular operators. The first involves projecting the Cartesian gradient operator onto the local tangent plane of the surface at some surface point. The surface gradient becomes,

 ∇Γf =∇f−∂nf→n =(I−→n→nT)∇f.

Such methods are classified as projection methods. The second approach involves extending the function f, into

such that it is constant in the normal direction at each point on the surface. As pointed out by  under such conditions the surface gradient and Cartesian gradient agree on the surface. Surface differential operators can then be replaced with the simpler Cartesian operators. The Orthogonal Gradient Method enforces this requirement by extending an RBF approximation of the function outside of the surface , originally having N points, using additional points. This increases the complexity of the resulting linear system from to . Also the accuracy of the method is influenced by the choice of an offset parameter which controls the proximity of the points to the surface.

In order to avoid the increased complexity of introducing 2N additional points to enforce the null gradient condition we propose a modified kernel which projects all radial vectors in the traditional RBF kernel onto the local tangent plane. We refer to this as the Cylindrical Kernel. Our modified kernel is intrinsically constant in the normal direction and greatly simplifies the construction of surface differential operators.

## 0.5 Cylindrical Basis Function (CBF)

###### Definition 1.

Given a point cloud sampled from a smooth manifold, and data from a smooth function defined at these points, a CBF interpolant of f is defined as

 c(x)=N∑j=1λjϕ(rj(x)) (9)

where,

 rj(x)=∥(x−xj)−→nj[(x−xj)⋅→nj]∥. (10)

and are the interpolation coefficients and kernel respectively, is the unit normal vector at the point and Eq.(10) is the cylindrical distance with respect to the Euclidean norm.

###### Lemma 1.

Let be a point cloud on a 2D-manifold embedded in then the CBF kernel satisfies for all

 ∇ϕ(ri(x))⋅→nyi=0∀i=1,2,⋯,N and x∈R3. (11)

Consider the Laplacian of f on ,

 △f=(∂nf→n+∂t1f→t1+∂t2f→t2)⋅(∂nf→n+∂t1f→t1+∂t2f→t2) (12)

expanding out the operator it is clear that the surface Laplacian, is equivalent to the regular Laplacian if . Therefore the following corollary is a consequence of Lemma 1

###### Corollary 1.

On a smooth 2D manifold the CBF interpolation, of a smooth function, satisfies that,

 ∇c(x) =∇Γc(x) (13) △c(x) =△Γc(x) (14)

This implies that surface differential operators can now be discretized by simply discretizing their corresponding Cartesian operators. We illustrate the great simplicity of this approach by computing the Laplacian of the cylindrical kernel. We present the main result as follows; Given , the cylindrical kernel, the surface Laplacian is given as,

 △ϕ(rij)=ϕ′′(||rij||)(1−(^rij⋅→nj)2)+ϕ′(||rij||)1+(^rij⋅→nj)2||rij|| (15)

where is the projected radial vector computed from the point to the center of the basis function . is the normalized projected vector.

## 0.6 Implementation

We outline the algorithm for discretizing the Laplace-Beltrami Operator on a 2D-manifold.

1. Obtain the surface normals. Surface normals can be computed for every point in the point cloud

using Principal Component Analysis(PCA). At any point

the goal is to find the coordinate system consisting of the normal and two tangential directions to the local coordinate plane within a neighborhood of . Let be a matrix whose rows represent the k nearest neighbors of . Denote by the mean of all the neighbors i.e

 c=1kk∑j=1xj

and let be row of the matrix . It is well know that the eigen vectors of the covariance matrix form an orthogonal basis which are the principal components . The eigen vector corresponding to the smallest eigen value is the desired normal direction to the local plane at .

2. Assemble the collocation matrix,

 A=⎡⎢ ⎢ ⎢ ⎢⎣ϕ(r11)ϕ(r12)…ϕ(r1n)ϕ(r21)ϕ(r22)…ϕ(r2n)4ϕ(rn1)ϕ(rn2)…ϕ(rnn)⎤⎥ ⎥ ⎥ ⎥⎦

where

 ϕ(rij):=ϕ(||(xi−xj)−→nj[(xi−xj)⋅→nj]||)

corresponds to the basis function centered at the node and evaluated at node .

3. Assemble the differential matrix,

 B=⎛⎜ ⎜ ⎜ ⎜⎝△ϕ(r11)△ϕ(r12)…△ϕ(r1n)△ϕ(r21)△ϕ(r22)…△ϕ(r2n)4△ϕ(rn1)△ϕ(rn2)…△ϕ(rnn)⎞⎟ ⎟ ⎟ ⎟⎠

where is as defined in Eq.(15)

4. The discrete LB operator is then given as

## 0.7 Numerical Experiments

### 0.7.1 Eigen Values of LB Operator

The eigen values and eigen functions of the LB operator provide intrinsic global information that can be used in characterizing the structure of surfaces[15, 14]. We therefore tested the performance of our method in obtaining the spectra of the LB operator on the unit sphere by solving the eigen value problem

 ΔSu=−λu, (16)

and compared our first 100 eigen values to the exact values. Our discrete operator was computed using a uniform sampling of the unit sphere with nodes ranging from 258-16386. We used a local implementation of the Gaussian radial basis function by choosing the shape parameter c, such that the function vanished outside the support radius. We performed the same experiment using Finite Element method and compared the convergence of both methods using the Euclidean distance as a measure of errors (Figure 2). As we refine the nodes we see from Figure 1 that our computed eigen values line almost perfectly with the true values. With 4098 nodes our computed eigen values are at a distance of 1.4 from the true values while the FEM values are 11.8 units away. Our results indicate that 4 times the number of nodes used for CBF would be required for FEM in order to achieve similar accuracy. We also observed a decline in the accuracy of our solution at 16386 nodes from 4098 nodes. This we believe is as a result of the significantly higher condition number of the collocation matrix. We had to consistently increase the size of the support domain as the number of nodes increased to achieve optimum results. We show the computational time with increasing numbers of nodes. Computations were performed using Matlab R2013a with an intel core i3-4130 8.40 GHZ processor with 8GB of ram. Figure 1: first 100 Eigen values of the Laplace Operator on a unit sphere. Number of nodes from left to right 258, 1026, 4098, 16386 Figure 2: Superior convergence of CBF and Computational Efficiency. CPU time was measured in secs.

### 0.7.2 Heat Diffusion on Sphere and Molecule

To verify that our discrete operator accurately captures the heat diffusion of the LB operator we solved the heat equation

 ut=ϵΔΓu (17)

on a unit sphere as well as a more complicated molecular surface. We used a Forward Euler time discritization scheme with . In Figure 3 we show that the diffusion, using our CBF discrete operator, matches well with the exact solution obtained using the spherical harmonics across a time duration of 1 sec. Figure 3: Heat Evolution on unit sphere with initial condition, the spherical harmonic Y01 and diffusion coefficient ϵ=1 over 1sec duration. (Above) Numerical solution using 4098 nodes. (Below) Exact evolution of sphereical harmonic Y01.

We used total of 7718 point on the molecule surface and allowed the heat to diffuse for 5 secs. As shown in Figure 4, The heat diffuses as expected from the heat source to neighboring regions. The initial heat distribution was chosen as a Gaussian bell, centered at a point . Figure 4: Heat equation on a molecule with the gaussian bell as initial distribution centered at one of its nodes

### 0.7.3 Laplace Smoothing

Laplace smoothing is a common technique used to smooth surface meshes . Here the Laplace operator acts as a weighted average of the coordinates of each point and its one ring neighbors to draw them towards the barycenter of their common region. To test the applicability of our CBF Laplacian, we used it in smoothing a noisy meshed sphere having 258 nodes. The nodes of the mesh were used as initial condition for the heat equation and updated for 60 iterations. The smoothing equation used was in the form where acts as the smoothing scale factor with being the diffusion coefficient and the time step size. is our CBF approximation to the Laplace Beltrami Operator. In Figure 5 we show the smoothing of the triangular mesh on sphere over 60 steps. To measure the effectiveness of our smoothing technique we measured the angle between the point normals at each stage of smoothing to the true normals of the perfectly smooth sphere. Both normals were computed using PCA with 6 neighbours. We choose not to use the position vectors of the point clouds on the perfect sphere as true normals in order to permit a more uniform comparison with the smoothed sphere. We see in Figure 6(a) and Figure 6(b) that the normals converge to the perfect sphere normals as we iterate to 60 steps. We didn’t notice any remarkable improvement in the smoothing after 60 steps. We find the method scales very well with size of mesh nodes as demonstrated in Figure 6(c).

## 0.8 Conclusion

In this paper we have introduced a new technique for discretizing surface differential operators on closed manifolds. By projecting the radial vector used in standard RBF kernels onto the local tangent plane of the manifold, we produce a modified kernel which is constant in the normal direction to the surface. Our modified kernel, when used to represent functions defined on manifolds, permits the simple replacement of surface differential operators by Cartesian operators. We have demonstrated through numerical experiments the superior performance of CBF in discretizing the Laplace Beltrami operator on the sphere by comparing the first 100 eigen values with the exact values. We also solved the heat equation on the sphere and a molecule and also showed that the discrete operator can be used to effectively smooth noisy surface meshes. Here we worked primarily with the Gaussian kernel. We will consider the effect of other basis functions such as the Matern kernel in future efforts. We will also investigate the approximation of other differential operators such as divergence operator.

## References

•  B.Levy, editor.

Laplace-Beltrami eigenfunctions towards an algorithm that ’understands’ geometry

, Matsushima, June 2006. International conference on shape modeling and applications, IEEE.
•  U. Diewald, T. Preusser, and M Rumpf. Anisotropic diffusion in vector field visualization on euclidean domains and surfaces. Visualization and Computer Graphics, IEEE Transactions on, 6(2):139–149, 2000.
•  David A. Field. Laplace smoothing and delaunay triangulations. Communications in Applied Numerical Methods, 4:709–712, 1988.
•  Natasha Flyer and Grady B. Wright. A radial basis function method for the shallow water equations on a sphere. Proc R. Soc. A, 465:1949–1967, 2009.
•  Yiming Fu, Jin Zhanga, and Liangjiao Wana. Application of the energy balance method to a nonlinear oscillator arising in the microelectromechanical system (mems). Current Applied Physics, 11(3):482–485, 2011.
•  G.Dziuk and C.Elliot. Surface finite elements for parabolic equations. Journal of Computational Mathematics, 25(4):385–407, 2007.
•  John B. Greer. An improvement of a recent eulerian method for solving pdes on general geometries. Journal of Scientific Computing, 29(3):321–352, 2006.
•  Michael Hofer and Helmut Pottmann. Energy-minimizing splines in manifolds. ACM Trans. Graph., 23(3):284–293, 2004.
•  J.Dorsey and P.Hanrahan. Digital materials and virtual weathering. Scientific American, 282:64–71, 2000.
•  Rongjie Lai, Yonggang Shi, K. Scheibel, S. Fears, R. Woods, A.W. Toga, and T.F Chan. Metric induced optimal embedding for instrinsic 3d shape analysis. In

Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on

, pages 2871–2878, 2010.
•  J. C. MISRA, A. SINHA, and G. C. SHIT. A numerical model for the magnetohydrodynamic flow of blood in a porous channel. Journal of Mechanics in Medicine and Biology, 11(03):547–562, 2011.
•  T.G Myers and J.P.F.Charpin. A mathematical model for atmospheric ice accretion and water flow on a cold surface.
•  Cecile Piret. The orthogonal gradients method: A radial basis function method for solving partial differential equations on arbitrary surfaces. Journal of Computational Physics, 231(14):4662–4675, 2012.
•  Martin Reuter, Silvia Biasotti, Daniela Giorgi, Giuseppe Patanè, and Michela Spagnuolo. Discrete laplace-beltrami operator for shape analysis and segmentation. Computers & Graphics, 33(3):381 – 390, 2009.
•  Martin Reuter, Franz-Erich Wolter, and Niklas Peinecke. Laplace-beltrami spectra as ’shape-dna’ of surfaces and solids. Computer Aided Design, 38(4):342–366, 2006.
•  Steven J. Ruuth and Barry Merriman. A simple embedding method for solving partial differential equations on surfaces. Computational Physics, 227(3):1943–1961, 2008.
•  Klaus Schittkowski. Parameter identification and model verification in systems of partial differential equations applied to transdermal drug delivery. Mathematics and Computers in Simulation, 79(3):521–538, 2008.
•  Jian Sun, Maks Ovsjanikov, and Leonidas Guibas. A consise and provably informative multi-scale signature based on heat diffusion. In Proceedings of the Symposium on Geometry Processing, pages 1383–1392, Aire-la-Ville, Switzerland, Switzerland, 2009. Eurographics Association.
•  J.P.F.Charpin T.G Myers. The flow and solidification of a thin fluid film on an arbitrary three-dimensional surface.
•  Greg Turk. Generating textures on arbitrary surfaces using reaction-diffusion. SIGGRAPH Comput. Graph., 25(4):289–298, 1991.
•  Andrew Witkin and Michael Kass. Reaction-diffusion textures. SIGGRAPH Comput. Graph., 25(4):299–308, 1991.
•  Hongkai Zhao. Geometric understanding of point clouds using laplace-beltrami operator. In Proceedings of the 2012 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 214–221, Washington, DC, USA, 2012. IEEE Computer Society.