|Home | About | Journals | Submit | Contact Us | Français|
We present a novel surface smoothing framework using the Laplace-Beltrami eigenfunctions. The Green’s function of an isotropic diffusion equation on a manifold is constructed as a linear combination of the Laplace-Beltraimi operator. The Green’s function is then used in constructing heat kernel smoothing. Unlike many previous approaches, diffusion is analytically represented as a series expansion avoiding numerical instability and inaccuracy issues. This proposed framework is illustrated with mandible surfaces, and is compared to a widely used iterative kernel smoothing technique in computational anatomy. The MATLAB source code is freely available at http://brainimaging.waisman.wisc.edu/~chung/lb.
In medical image processing and analysis, anatomical manifolds are commonly represented as triangular meshes. Procedures such as image segmentation, mesh construction and geometric computations on the meshes introduce geometric noise to the mesh coordinates. Therefore, it is imperative to reduce the mesh noise while preserving the geometric details of the object for various applications.
Many approaches have been proposed for smoothing surface data and coordinates. Probably the most widely used method is to numerically solve diffusion equations on anatomical surfaces [1,8]. However, these approaches use discretization schemes which tend to suffer numerical instability and inaccuracy. Iterated kernel smoothing is also a widely used method in computational anatomy [3,6] and computer vision . Kernel weights are spatially adapted to follow the shape of the heat kernel in a discrete fashion along a manifold. In the tangent space, the heat kernel is approximated linearly using the Gaussian kernel. However, this process compounds the linearization error at each iteration.
In this paper, we propose to construct the heat kernel analytically using the eigenfunctions of the Laplace-Beltrami operator, avoiding the need for the linear approximation used in [3,6]. Although solving for the eigenfunctions requires the finite element method, the proposed method is analytic in a sense that heat kernel smoothing is formulated as a series expansion explicitly. Thus, the proposed method avoids the numerical instability associated with solving the diffusion equations numerically [1,8]. Although there are many papers on solving diffusion equations on arbitrary triangular meshes [1,8,14], this is the first paper that explicitly and correctly constructed heat kernel for an arbitrary surface and solved heat diffusion using the eigenfunctions of the Laplace-Beltrami operator.
Consider a real-valued C2 measure Y defined on a manifold 3. We assume the following additive model on Y:
where θ(p) is the unknown mean signal and ε(p) is a zero-mean Gaussian random field. We may assume further Y L2(), the space of square integrable functions on with the inner product f, g = ∫ f(p)g(p)dμ(p), where μ is the Lebesgue measure such that μ() is the volume of . Solving
for the Laplace-Beltrami operator Δ on , we order eigenvalues 0 = λ0 < λ1 ≤ λ2 ≤ ··· and corresponding eigenfunctions ψ0, ψ1, ψ2, ···. Then, the eigenfunctions ψj form an orthonormal basis in L2() .
Using the eigenfunctions, heat kernel Kσ(p, q) is analytically written as
where βj = Y, ψj are Fourier coefficients. This is taken as the estimate for θ. Since the expansion (4) is a unique solution to isotropic heat diffusion [3,12], we can avoid the need to solve the diffusion equation using less stable numerical schemes such as the finite difference method [1,8,14].
In this new analytic framework, we need to compute the terms in (4). We first solve for the eigensystem (2) and obtain λj and ψj. To estimate βj, we let σ = 0 in (4). Then Kσ becomes the Dirac-delta function and we obtain
This is the usual Fourier expansion of Y. The finite series expansion of (5) will be then used in estimating the Fourier coefficients in the least squares fashion.
Since the closed form expression for the eigenfunctions of the Laplace-Beltrami operator on an arbitrary curved surface is unknown, the eigenfunctions are numerically calculated by discretizing the Laplace-Beltrami operator. To solve the eigensystem (2), we need to discretize it on a triangular mesh using the Cotan formulation [4,11]. In particular, Qiu et al.  presented a similar Cotan discretization and constructed splines on a manifold using the eigenfunctions; however, there is no direct mathematical relation between splines and heat kernel smoothing.
Using the Cotan formulation, (2) is simplified as the generalized eigenvalue problem:
where C is the stiffness matrix, A is the mass matrix, and ψ is the unknown eigenfunction evaluated at mesh vertices. Because C and A are large sparse matrices, we have solved the problem using the Implicitly Restarted Arnoldi Method [7,9] without consuming large amount of memory and time for zero entries. The MATLAB code is given at http://brainimaging.waisman.wisc.edu/~chung/lb. Fig. 1 shows the first few eigenfunctions for a mandible surface.
Let be the subspace spanned by up to k-th degree basis. Then an arbitrary measurement Y is estimated in the subspace by minimizing the sum of squared residual:
Consider the triangular mesh for with Nv nodes, and let β = (β0, ···, βk )′ and Y = (Y(p1), ···, Y(pNv))′, for k ≤ Nv. Then, we can represent (7) as the normal equation,
where Ψ = (Ψ0, ···, Ψk) and Ψj = (ψj(p1), ···, ψj(pNv))′, and β are estimated in the least squares fashion, = (Ψ′Ψ)−1Ψ′Y. Since the size of matrix Ψ′Ψ can become fairly large when there is a need to obtain large number of basis, it may be difficult to directly invert the matrix. So we have adopted a more general iterative strategy to overcome possible computational bottleneck for large k.
The Fourier coefficients are estimated based on an iterative procedure that utilizes the orthonormality of the eigenfunctions . Decompose the subspace into smaller subspaces as the direct sum = ··· , where each subspace is the projection of along the j-th eigenfunction. Instead of directly solving the normal equation (8), we project the normal equations into a smaller subspace and find the corresponding βj in an iterative fashion from increasing the degree from 0 to k.
At degree k = 0, we write Y = Ψ0β0+r0, where r0 is the residual of estimating Y in subspace . Then, we estimate β0 by minimizing the residual in the least squares fashion:
At degree j, we have
where the previous residual rj−1 = Y − Ψ00 ··· −Ψj−1j−1. The next residual rj is then estiamted as . The optimal stopping rule is determined if the decrease of the root mean squared errors (RMSE), , is statistically significant using the F-test (Fig. 2). We compute the F statistic at each degree, and find the degree of expansion where corresponding p-value first becomes bigger than the pre-specified significance 0.01.
The CT images were obtained using several different models of GE multi-slice helical CT scanners. The CT scans were acquired directly in the axial plane with a 1.25 mm slice thickness, matrix size of 512 × 512 and 15–30 cm FOV. Image resolution varied and was in the range of 0.29 to 0.59 mm as determined by the ratio of field of view (FOV) divided by the matrix. CT scans were converted to DICOM format and subsequently Analyze 8.1 software package (AnalyzeDirect, Inc., Overland Park, KS) was used in segmenting binary mandibular structure based on histogram thresholding. By checking the Euler characteristic, holes in mandible images were automatically filled up using morphological operations. This process was necessary to make the mandible binary volume to be topologically equivalent to a solid sphere and and produces the surface mesh that is topologically equivalent to a sphere.
We applied the proposed method in smoothing a mandibular surface. The optimal eigenfunction expansion was determined using the F-test at α = 0.01. Fig. 2 shows the plot of the RMSE for varying degrees between 5 to 200. As the degree k increases, the RMSE for each coordinate rapidly decreases and starts to flatten out at a certain degree.
The numerical implementation was done with MATLAB 7.9 in 2 × 2.66 GHz Quad-Core Intel Xeon processor MAC PRO with 32 GB memory. For a mesh with 22050 vertices, the entire process took approximately 105 seconds: 85 seconds for setting up the generalized eigenvalue problem (6), 10 seconds for actually solving (6), 0.1 seconds for the IRF, and 9 seconds for finding the optimal degree.
The proposed heat kernel smoothing was compared against widely used iterated kernel smoothing [3,6] of which the MATLAB code is given in http://www.stat.wisc.edu/mchung/softwares/hk/hk.html. In iterated kernel smoothing, the weights of the kernel are spatially adapted to follow the shape of heat kernel in discrete fashion along a surface mesh. Smoothing with large bandwidth is broken into iterated smoothing with smaller bandwidths:
For small σ, using the parametrix expansion , we can approximate heat kernel locally using the Gaussian kernel for small bandwidth:
where d(p, q) is the geodesic distance between p and q. For sufficiently small bandwidth, all the kernel weights are concentrated near the center, so the first neighbors of a given vertex in a mesh is used in the approximation.
By taking the proposed heat kernel smoothing as the ground truth, we were able to determine the performance of iterated kernel smoothing. Due to the lack of the ground truth, there are no studies in the literature validating the performance of iterated kernel smoothing except . For heat kernel smoothing, we used the bandwidth σ = 1 and eigenfunctions up to k = 139 degree. For iterated kernel smoothing, we varied the number of iterations 1 ≤ m ≤ 200 with the correspondingly smaller bandwidth 1/m to have the effective bandwidth of 1. The performance of the iterated kernel smoothing depended on the number of iterations, as shown in the plot of RMSE of x-coordinate over the number of iterations (Fig. 3 right). The RMSE was up to 0.5430 and it did not decrease even when we increase the number of iterations. This comparison directly demonstrates for the first time, the limitation of iterated heat kernel smoothing which does not converge to heat diffusion.
In another comparison (Fig. 4), we numerically constructed a heat kernel with a small bandwidth 0.2 as a sample data. Then we performed the additional iterated kernel smoothing with σ = 0.2 and m = 49 on the sample data to obtain the effective smoothing bandwidth of 10. The result was compared with the proposed heat kernel smoothing with the bandwidth 9.8 on the sample data making the effective bandwidth of 10. From Fig. 4, we immediately see that the shapes of two kernels are different. This visually demonstrates iterated kernel smoothing substantially diverges from heat kernel smoothing.
We present a novel surface data smoothing framework where a smoothed surface measure is represented as a weighted linear combination of Laplace-Beltrami eigenfunctions. The expansion analytically solves isotropic heat diffusion. Taking the expansion as the ground truth, the proposed method is compared against widely used iterated kernel smoothing. The proposed method outperforms iterated kernel smoothing in accuracy. The divergence of iterated kernel smoothing is visually represented as kernel shapes on the mandible confirming the superiority of this proposed framework.
This work was supported by NIH-grant R01 DC6282 from the National Institute of Deafness and other Communicative Disorders, a core grant P-30 HD03352 to the Waisman Center from the National Institute of Child Health and Human Development, WCU-grant from the government of Korea to the Department of Brain and Cognitive Sciences, Seoul National University. We thank Dongjun Chung, Lindell R. Gentry, Mike S. Schimek, Katelyn J. Kassulke and Reid B. Durtschi for assistance with image acquisition and segmentation.