Unified Heat Kernel Regression for Diffusion, Kernel Smoothing and Wavelets on Manifolds and Its Application to Mandible Growth Modeling in CT Images

Unified Heat Kernel Regression for Diffusion, Kernel Smoothing and Wavelets on Manifolds and Its Application to Mandible Growth Modeling in CT Images

Moo K. Chung Anqi Qiu Seongho Seo Houri K. Vorperian University of Wisconsin, Madison National University of Singapore Seoul National University, Korea

We present a novel kernel regression framework for smoothing scalar surface data using the Laplace-Beltrami eigenfunctions. Starting with the heat kernel constructed from the eigenfunctions, we formulate a new bivariate kernel regression framework as a weighted eigenfunction expansion with the heat kernel as the weights. The new kernel method is mathematically equivalent to isotropic heat diffusion, kernel smoothing and recently popular diffusion wavelets. The numerical implementation is validated on a unit sphere using spherical harmonics. As an illustration, the method is applied to characterize the localized growth pattern of mandible surfaces obtained in CT images between ages 0 and 20 by regressing the length of displacement vectors with respect to a surface template.

Heat kernel regression, Laplace-Beltrami eigenfunctions, Diffusion wavelets, Surface-based morphometry, Mandible growth
journal: Medical Image Analysis

1 Introduction

In medical imaging, anatomical surfaces extracted from MRI and CT are often represented as triangular meshes. The image segmentation and surface extraction processes themselves are likely to introduce noise to the mesh coordinates. It is therefore imperative to reduce the mesh noise while preserving the geometric details of the anatomical structures for various applications.

Diffusion equations have been widely used in image processing as a form of noise reduction since 1990 (Perona and Malik, 1990). Numerous techniques have been developed for surface fairing and mesh regularization (Sochen et al., 1998; Malladi and Ravve, 2002; Tang et al., 1999; Taubin, 2000) and surface data smoothing (Andrade et al., 2001; Chung et al., 2001; Chung and Taylor, 2004; Cachia et al., 2003a, b; Chung et al., 2005; Joshi et al., 2009). Isotropic heat diffusion on surfaces has been introduced in brain imaging for subsequent statistical analysis involving the random field theory (RFT) that assumes an isotropic covariance function as a noise model (Andrade et al., 2001; Chung and Taylor, 2004; Cachia et al., 2003a, b). Since then, isotropic diffusion has been the standard smoothing technique.

Iterated kernel smoothing has been another widely used method in approximately solving diffusion equations on surfaces (Chung et al., 2005; Han et al., 2006). It is often used in smoothing anatomical surface data including cortical curvatures (Luders et al., 2006b; Gaser et al., 2006), cortical thickness maps (Luders et al., 2006a; Bernal-Rusiel et al., 2008), hippocampus surfaces (Shen et al., 2006; Zhu et al., 2007) and magnetoencephalography (MEG) (Han et al., 2007) and functional-MRI (Hagler Jr. et al., 2006; Jo et al., 2007) data on the brain surface. Due to its simplicity, it is the most widely used form of surface data smoothing in brain imaging. In iterated kernel smoothing, kernel weights are spatially adapted to follow the shape of the heat kernel in a discrete fashion along a manifold. In the tangent space of the manifold, the heat kernel with a small bandwidth can be approximated linearly using the Gaussian kernel. The heat kernel with a large bandwidth is then constructed by iteratively applying the Gaussian kernel with the small bandwidth. However, this process compounds the linearization error at each iteration.

We propose a new kernel regression framework that constructs the heat kernel analytically using the eigenfunctions of the Laplace-Beltrami (LB) operator, avoiding the need for the linear approximation used by Chung et al. (2005) and Han et al. (2006). Although a few studies have introduced the heat kernel in computer vision and machine learning, they mainly used the heat kernels to compute shape descriptors or to define a multi-scale metric (Belkin et al., 2006; Sun et al., 2009; Bronstein and Kokkinos, 2010; de Goes et al., 2008). These studies did not use the heat kernels in regressing functional data on manifolds. This is the first study to use the heat kernel in the form of regression for the subsequent statistical analysis. There have been significant developments in kernel methods in the machine learning community (Schölkopf and Smola, 2002; Nilsson et al., 2007; Shawe-Taylor and Cristianini, 2004; Steinke and Hein, 2008; Yger and Rakotomamonjy, 2011). However, to the best of our knowledge, the heat kernel has never been used in such frameworks. Most kernel methods in machine learning deal with the linear combination of kernels as a solution to penalized regressions. On the other hand, our kernel regression framework does not have a penalized cost function.

Wavelets have recently been popularized for surface and graph data. For instance, spherical wavelets were used on brain surface data already mapped onto a sphere (Nain et al., 2007; Bernal-Rusiel et al., 2008). Since the wavelet basis has local supports in both space and scale, the wavelet coefficients provide shape features that describe local shape variation at a variety of scales and spatial locations. However, spherical wavelets require a smooth mapping from the surface to a unit sphere, thus introducing a serious metric distortion that compounds subsequent statistical parametric maps (SPM). Furthermore, such basis functions are only orthonormal for data defined on the sphere and result in a less parsimonious representation for data defined on other surfaces compared to the intrinsic LB-eigenfunction expansion (Seo and Chung, 2011). To remedy the limitations of spherical wavelets, the diffusion wavelet transform on graph data structures has been proposed (Antoine et al., 2010; Coifman and Maggioni, 2006; Hammond et al., 2011; Kim et al., 2012).

The primary methodological contribution of this study is the establishment of a unified regression framework that combines the diffusion-, kernel- and wavelet-based methods in a coherent mathematical fashion for scalar data defined on manifolds. We unify the apparent differences between the methods while providing detailed theoretical justifications. This paper extends the work by Kim et al. (2011), which introduced heat kernel smoothing to smooth out surface noise in the hippocampus and amygdala. Although the idea of diffusion wavelet transform for surface meshes was explored by Kim et al. (2012), the relationship between the wavelet transform and the proposed kernel regression was not investigated. For the first time, the mathematical equivalence between the two constructs is explained.

The proposed kernel regression framework was subsequently applied in characterizing the growth pattern of the mandible surfaces obtained in CT and identifying the regions of the mandible that show the most significant localized growth. The length of the displacement vector field was regressed over the mandible surface to increase the signal-to-noise ratio (SNR) and hence statistical sensitivity. To our knowledge, this is the first growth modeling of the mandible surface in a continuous fashion without using anatomic landmarks.

2 Methods

2.1 Isotropic Diffusion on Manifolds

Consider a functional measurement observed at each point on a compact manifold . We assume the following linear model on :


where is the unknown mean signal to be estimated and is a zero-mean Gaussian random field. We may assume further , the space of square integrable functions on with the inner product:


where is the Lebesgue measure such that is the total area of .

Imaging data such as electroencephalography (EEG), magnetoencephalography (MEG) (Han et al., 2007) and functional-MRI (Hagler Jr. et al., 2006; Jo et al., 2007), and anatomical data such as cortical curvatures (Luders et al., 2006b; Gaser et al., 2006), cortical thickness (Luders et al., 2006a; Bernal-Rusiel et al., 2008) and surface coordinates (Chung et al., 2005) can be considered as functional measurements. Functional measurements are expected to be noisy and require filtering to boost signal.

Surface measurements have often been filtered using the isotropic diffusion equation of the form (Andrade et al., 2001; Chung et al., 2001; Cachia et al., 2003a; Rosenberg, 1997)


where is the Laplace-Beltrami operator defined on manifold . The diffusion time controls the amount of smoothing. It can be shown that the unique solution of (3) is given by a kernel convolution as follows. A Green’s function or a fundamental solution of the Cauchy problem (3) is given by the solution of the following equation:


where is the Dirac delta function. The heat kernel is a Green’s function of (4) (Evans, 1998), i.e.,

Since the differential operators are linear in (4), we can further convolve the terms with the initial data such that


Hence is a solution of (3).

Figure 1: The heat kernel shape with bandwidths 0.025 (left), 1.25 (middle) and 5 (right) on a mandible surface. The level sets of the heat kernel form geodesic circles.

2.2 Diffusion Smoothing

Isotropic diffusion (3) has been solved by various numerical techniques (Chung, 2001; Andrade et al., 2001; Cachia et al., 2003a, b; Chung and Taylor, 2004). For diffusion smoothing, the diffusion equation needs to be discretized using the cotan formulation (Chung, 2001; Chung and Taylor, 2004; Qiu et al., 2006). Since there are many different cotan formulations, we follow the formulation introduced in Chung (2001). The diffusion equation (3) is discretized as


where is the vector of measurements over all mesh vertices at time . is the stiffness matrix and is the global coefficient matrix, which is the assemblage of individual element coefficients. The sparse matrices and are explicitly given as follows.

Let and denote two triangles sharing the vertex and its neighboring vertex in a mesh. Let the two angles opposite to the edge containing and be and respectively for and . The off-diagonal entries of the stiffness matrix are

if and are adjacent and otherwise. denotes the area of a triangle. The diagonal entries are summed as The off-diagonal entries of the global coefficient matrix are

if and are adjacent and otherwise. The diagonal entries are similarly given as the sum

The ordinary differential equation (5) is then further discretized at each point using the forward finite difference scheme:


where is the estimated Laplacian obtained from the -th row of . For the forward Euler scheme to converge, we need to have a sufficiently small step size (Chung, 2001).

2.3 Iterated Kernel Smoothing

The diffusion equation (3) can be approximately solved by iteratively performing Gaussian kernel smoothing (Chung et al., 2005). The weights of the kernel are spatially adapted to follow the shape of the heat kernel along a surface mesh. Heat kernel smoothing with a large bandwidth can be broken into iterated smoothing with smaller bandwidths (Chung et al., 2005):


Then using the parametrix expansion (Rosenberg, 1997; Wang, 1997), we approximate the heat kernel with the small bandwidth locally using the Gaussian kernel:


where is the geodesic distance between and . For sufficiently small bandwidth , all kernel weights are concentrated near the center, so the first neighbors of a given mesh vertex are sufficient for approximation. Unfortunately, this approximation compounds error at each iteration. For numerical implementation, we used the normalized truncated kernel given by


where are neighboring vertices of . Denote the truncated kernel convolution as


The iterated heat kernel smoothing is then defined as

2.4 Heat Kernel Regression

We present a new regression framework for solving the isotropic diffusion equation (3). Let be the Laplace-Beltrami operator on . Solving the eigenvalue equation


we order the eigenvalues

and corresponding eigenfunctions (Rosenberg, 1997; Chung et al., 2005; Lévy, 2006; Shi et al., 2009). The first eigenvalue and eigenfunction are trivially given as and . It is possible to have multiple eigenfunctions corresponding to the same eigenvalue.

The eigenfunctions form an orthonormal basis in . There is extensive literature on the use of eigenvalues and eigenfunctions of the Laplace-Beltrami operator in medical imaging and computer vision (Lévy, 2006; Qiu et al., 2006; Reuter et al., 2009; Reuter, 2010; Zhang et al., 2007, 2010). The eigenvalues have been used in caudate shape discriminators (Niethammer et al., 2007). Qiu et al. used eigenfunctions in constructing splines on cortical surfaces (Qiu et al., 2006). Reuter used the topological features of eigenfunctions (Reuter, 2010). Shi et al. used the Reeb graph of the second eigenfunction in shape characterization and landmark detection in cortical and subcortical structures (Shi et al., 2008, 2009). Lai et al. used the critical points of the second eigenfunction as anatomical landmarks for colon surfaces (Lai et al., 2010). Since the direct application of eigenvalues and eigenfunctions as features of interest is beyond the scope of this paper, we will not pursue the issue in detail here.

Using the eigenfunctions, the heat kernel is defined as


where is the bandwidth of the kernel. Figure 1 shows examples of the heat kernel with different bandwidths. Then the heat kernel regression or heat kernel smoothing of functional measurement is defined as


where are Fourier coefficients (Chung et al., 2005) (Figure 2). The kernel smoothing is taken as an estimate for the unknown mean signal in (1). The degree of truncation of the series expansion can be automatically determined using the forward model selection procedure.

Figure 2: Schematic of heat kernel smoothing. Given functional data on a surface, we first compute the eigenfunctions and the Fourier coefficients . Then, we combine all the terms and reconstruct the functional signal back.

Unlike previous approaches to heat diffusion (Andrade et al., 2001; Chung and Taylor, 2004; Joshi et al., 2009; Tasdizen et al., 2006), the proposed method avoids the direct numerical discretization of the diffusion equation. Instead we discretize the basis functions of the manifold by solving for the eigensystem (11) to obtain and .

2.5 Diffusion Wavelet Transform

We can establish the equivalence between the proposed kernel regression and recently popular diffusion wavelets. This mathematical equivalence eliminates the need for constructing wavelets using complicated computational machinery as has often been done in previous studies (Antoine et al., 2010; Hammond et al., 2011; Kim et al., 2012), and offers a simpler but more unified alternative.

Consider a wavelet basis obtained from a mother wavelet with scale and translation parameters and respectively in a Euclidean space:

Generalizing the idea of scaling a mother wavelet in Euclidean space to a curved surface is trivial. The difficulty arises when one tries to translate a mother wavelet on a curved surface since it is unclear how to define translation along the surface. If one tries to modify the existing spherical wavelet framework to an arbitrary surface (Nain et al., 2007; Bernal-Rusiel et al., 2008), one immediately encounters the problem of establishing regular grids on an arbitrary surface. Recent work based on diffusion wavelets bypass this problem by taking the bivariate kernel as a mother wavelet (Antoine et al., 2010; Hammond et al., 2011; Mahadevan and Maggioni, 2006; Kim et al., 2012).

For some scale function that satisfies the admissibility conditions, diffusion wavelet at position and scale is given by

where and are eigenvalues and eigenfunctions of the Laplace-Beltrami operator. The wavelet transform is then given by


By letting , we have the heat kernel as the wavelet, i.e.,

The bandwidth of the heat kernel is the scale parameter, while the translation is achieved by shifting one argument in the bivariate heat kernel. Subsequently, wavelet transform (14) can then be rewritten as


with . The expression (15) is the finite truncation of the heat kernel regression in (12). Hence, diffusion wavelet analysis can be simply performed within the proposed heat kernel regression framework without any additional wavelet machinery. We therefore do not distinguish between the heat kernel regression and the diffusion wavelet transform.

Although the heat kernel regression is constructed using the global basis functions , surprisingly the kernel regression at each point coincides with the wavelet transform at that point. Hence, it also inherits the localization property of wavelets at that point. This is clearly demonstrated in the simulation study shown in Figure 3, where a step function of value 1 in the circular band (angle from the north pole) and of value 0 outside of the band was constructed. Note that, on a sphere, the Laplace-Beltrami operator is the spherical Laplacian and its eigenfunctions are spherical harmonics of degree and order . Then the step function was reconstructed using the spherical harmonic series expansion

where the spherical harmonic coefficients were obtained by the least squares estimation (LSE). On the unit sphere, we used the heat kernel regression of the form

with small bandwidth and degree . The spherical harmonic expansion shows severe ringing artifacts compared to the kernel regression, which inherits the localization power of wavelets. Thus the Gibbs phenomenon was not significantly visible.

Figure 3: Gibbs phenomenon (ringing artifacts) is visible in the spherical harmonic series expansion with degree 78 via the least squares estimation (LSE) of the step function defined on a sphere. In contrast, the heat kernel regression with the same degree and bandwidth 0.0001 shows fewer visible artifacts.

2.6 Parameter Estimation in Heat Kernel Regression

Since the closed form expression for the eigenfunctions of the Laplace-Beltrami operator on an arbitrary surface is unknown, the eigenfunctions are numerically computed by discretizing the Laplace-Beltrami operator. To solve the eigensystem (11), we need to discretize it on mandible triangular meshes using the cotan formulation (Chung, 2001; Chung and Taylor, 2004; Shi et al., 2009; Qiu et al., 2006; Lévy, 2006; Reuter et al., 2006, 2009; Rustamov, 2007; Zhang et al., 2007; Vallet and Lévy, 2008; Wardetzky, 2008).

Among the many different cotan formulations used in computer vision and medical image analysis, we used the formulation given in Chung (2001) and Qiu et al. (2006). It requires discretizing (11) as the following generalized eigenvalue problem:


where the global coefficient matrix is the assemblage of individual element coefficients and is the stiffness matrix. We solved (16) using the Implicitly Restarted Arnoldi Method (Hernandez et al., 2006; Lehoucq et al., 1998) without consuming large amounts of memory and time for sparse entries. Figure 4 shows the first few eigenfunctions for a mandible surface.

Figure 4: Eigenfunctions of various degrees for a sample mandible surface. The eigenfunctions are projected on the surface smoothed by the proposed heat kernel smoothing with bandwidth and degree . The smoothed surface is obtained by heat kernel smoothing applied to the coordinates of the surface mesh with the same parameter while preserving the topology of mesh. The first eigenfunction is simply . The color scale is thresholded for better visualization.

Once we obtain the eigenfunctions numerically, we estimate the kernel regression parameters using the least squares estimation (LSE) technique. Note , the Fourier coefficients with respect to basis . are then estimated simultaneously by minimizing the sum of squared residual:


The least squares method is often used in estimating the coefficients in spherical harmonic expansion (Shen et al., 2004; Styner et al., 2006; Chung et al., 2008). Suppose we have mesh vertices . Let

be the surface measurements over all vertices. Denote the -th eigenfunction evaluated at vertices as

The minimum in (17) is achieved when


where is the matrix of size . The LSE estimation of coefficients is then given by


Since it is expected that the number of mesh vertices is substantially larger than the number of eigenfunctions to be used, is well conditioned and invertible. The numerical implementation is available at http://www.stat.wisc.edu/~mchung/mandible with the full data set used in the study.

2.7 Random Field Theory (RFT)

Once we have smoothed functional data on a surface, we apply the statistical parametric mapping (SPM) framework for analyzing and visualizing statistical tests performed on the template surface that is often used in structural neuroimaging studies (Andrade et al., 2001; Lerch and Evans, 2005; Wang et al., 2010; Worlsey et al., 1995; Yushkevich et al., 2008). Since test statistics are constructed over all mesh vertices on the mandible, multiple comparisons need to be accounted for using the RFT (Taylor and Worsley, 2007; Worlsey et al., 1995; Worsley et al., 2004). The RFT assumes the measurements to be a smooth Gaussian random field. Heat kernel smoothing will make data smoother and more Gaussian and enhance the SNR (Chung et al., 2005). The proposed kernel regression framework can then be naturally integrated into the RFT-based statistical inference approach (Taylor and Worsley, 2007; Worsley et al., 2004; Worlsey et al., 1995).

We assume is an unknown group level signal and is a zero-mean unit-variance Gaussian random field in (1). The model assumptions are not as restrictive as it seems since we can always normalize the data in this fashion. We further assume the random field is the convolution of heat kernel on Gaussian white noise with bandwidth , i.e., . Previously, the smoothness of noise, i.e., kernel bandwidth , was estimated using the RESEL (resolution element) technique, which requires estimating the quantity along mesh surfaces, which can introduce a bias (Worsley et al., 1999; Hayasaka et al., 2004; Kilner and Friston, 2010). Thus, surface data is often smoothed with bandwidth that is sufficiently larger than so that any high frequency noise smaller than is masked out. This provides a motivation for developing the heat kernel regression framework.

In (1), we are interested in determining the significance of , i.e.,


Note that any point that gives is considered as signal. The hypothesis (20) is an infinite dimensional multiple comparisons problem for continuously indexed hypotheses over the manifold . The underlying group level signal is estimated using the proposed heat kernel regression. Subsequently, a test statistic is often given by a T- or F-field (Worsley et al., 2004; Worlsey et al., 1995).

The multiple comparisons corrected -value is then computed through by the RFT (Adler, 1981; Cao and Worsley, 2001; Taylor and Worsley, 2007; Worsley, 2003). For the -field with and degrees of freedom defined on 2D manifolds , it is known that


for a sufficiently large threshold , where is the -th Minkowski functional of and is the -th Euler characteristic (EC) density of (Worsley et al., 1998). The Minkowski functionals are given by

The EC-density for -field is then given by

where is the cumulative distribution function of -stat with and degrees of freedom. The second order term dominates the expression (21) and it explicitly has the bandwidth of the kernel regression, thus incorporating the proposed kernel framework into the RFT.

3 Experiments

3.1 CT Image Preprocessing

We applied the proposed smoothing method to CT images of mandibles obtained from several different models of GE multi-slice helical CT scanners. The CT scans were acquired directly in the axial plane with mm slice thickness, matrix size of and cm field of view (FOV). Image resolution varied as voxel size ranged from 0.25 mm to 0.49 mm as determined by the ratio of FOV divided by the matrix. CT scans were converted to DICOM format and Analyze software package (AnalyzeDirect, Inc., Overland Park, KS) was then used in segmenting binary mandible structure based on histogram thresholding.

Image acquisition and processing artifacts and partial voluming produce topological defects such as holes and handles in any medical image. In CT images of mandibles, unwanted cavities, holes and handles in the binary segmentation mainly result from differences in CT intensity between relatively low density mandible and teeth and more dense cortical bone and the interior trabecular bone (Andresen et al., 2000; Loubele et al., 2006). In mandibles, these topological noises can appear in thin or cancellous bone, such as in the condylar head and posterior palate (Stratemann et al., 2010). An example is shown in Figure 5, where the tooth cavity forms a bridge over the mandible. If we apply the isosurface extraction on the topologically defective segmentation results, the resulting surface will have many tiny handles (Wood et al., 2004; Yotter et al., 2009). These handles complicate subsequent surface mesh operations such as smoothing and parameterization. It is thus necessary to correct the topology by filling the holes and removing handles. If we correct such topological defects, it is expected that the resulting isosurface is topologically equivalent to a sphere.

Figure 5: Upper panel: Topological correction on mandible binary segmentation and surface. Disjointed tiny speckles of noisy components are removed by labeling the largest connected component, and holes and handles are removed by the 2D morphological closing operations applied sequentially to each image dimension (a b c d). Lower panel left: Surface reconstruction showing holes and handles in the teeth regions. The isosurface has Euler characteristic . Lower panel right: After the correction with .

Various topology correction techniques have been proposed in medical image processing. Rather than attempting to repair the topological defects of the already extracted surfaces (Wood et al., 2004; Yotter et al., 2009), we performed the topological simplification on the volume representation directly using morphological operations (Guskov and Wood, 2001; Van Den Boomgaard and Van Balen, 1992; Yotter et al., 2009). The direct correction on surface meshes can possibly cause surfaces to intersect each other (Wood et al., 2004). By checking the Euler characteristic, the holes were automatically filled up using morphological operations to make the mandible binary volume topologically equivalent to a solid sphere. All areas enclosed by the higher density bone included in the mandible definition are morphed into being in the definition of the mandible object. The hole-filled images were then converted to surface meshes via the marching cubes algorithm.

Figure 6: Cavity patching by topological closing operations. Left: Surface model of the binary volume that simulates a tooth cavity. Middle: The 3D image volume based closing operation does not properly patch the cavity region. Right: The 2D image slice based closing operation patches the cavity region properly.

In our semi-automated algorithm, we first removed the speckles of noise components by identifying the largest connected component in the binary volume. The resulting binary mandible volume was a single connected component with many small holes and handles. Then we applied the morphological closing operation in each 2D slice of CT images one by one in all three axes. Recombining the topology-corrected 2D slices resulted in topologically correct surface meshes (Figure 5). We used 2D topological closing operations mainly because of better performance and relatively simpler implementation than 3D topological closing operations. In 2D topological operations, we need to consider only 8 neighboring pixels compared to 26 neighboring voxels in a 3D image volume. There are many large concave regions left out by teeth and fillings. These regions may not be closed with a single 3D closing operation but can be easily patched up with a sequence of 2D closing operations, which put more constraints on the underlying topology. Instead of performing a single 3D closing operation that may not work, we sequentially performed 2D closing operations in each image slice in the x- , y- and z-directions. In Figure 5 upper panel, (a) is the binary volume before any closing operation, (b) is after the 2D closing operation in the x-direction, (c) is after the 2D closing operation in the y-direction and (d) is after the 2D closing operation in the z-direction. Each time we apply the 2D closing operation, the holes get smaller. Figure 6 shows a simulated cavity example that was not patched by the 3D closing operation (Van Den Boomgaard and Van Balen, 1992) but was easily patched by the sequential application of 2D closing operations. Note that any 3D object, whose every 2D cross-section is topologically equivalent to a solid disk, is topologically equivalent to a solid sphere. The problem of 3D topology correction can be thus reduced to a much simpler problem of 2D topology correction of multiple slices. Unfortunately, we cannot perform the closing operations to infinitely many possible 2D cross-sections in 3D image volumes. Therefore, we applied the 2D operations in the three axial directions. So there is a small chance the operation may not work in practice. Therefore, at the end of the processing, we performed a visual inspection of the processed volume. Further, we double checked the Euler characteristic of the resulting surface meshes. Note that for each triangle, there are three edges. For a closed surface topologically equivalent to a sphere, two adjacent triangles share the same edge. The total number of edges is thus , where is the number of faces. If is the total number of vertices, the Euler characteristic of a sphere is given by . Thus, we checked if the resulting mesh satisfies the condition . 77 binary mandible volumes used in the study produced the topologically correct surfaces without exception. Figure 5 lower panel shows an example of before and after the topology correction.

3.2 Validation of Heat Kernel Smoothing on Mandible Surfaces

The accuracy of the heat kernel construction using LB-eigenfunctions on a unit sphere using the ground truth can be found in Kim et al. (2011) so the results are not reproduced here. In this paper, we compared the performance of the proposed kernel regression against iterated kernel smoothing and diffusion smoothing using the real mandible data. The comparison results are similar for all 77 mandible surfaces used in the study, so only the results for one representative mandible surface is shown.

Figure 7: (a) Plot of the RMSE of the heat kernel regression against iterated kernel smoothing for coordinates (middle), (top) and (bottom) over the number of iterations up to 200. For the heat kernel regression, and are used. Iterated kernel smoothing does not converge to heat diffusion. RMSE between the kernel regression and the diffusion smoothing is smaller than 0.0046 so they are not displayed in the plot. (b) The squared difference between the kernel regression and the iterated kernel smoothing. The difference is mainly localized in high curvature areas, where the Gaussian kernel used in the iterated kernel smoothing fails to approximate the heat kernel. (c) The squared difference between the kernel regression and the diffusion smoothing. There are almost no visible differences.

The x, y and z coordinates of a mandible surface are treated as functional measurements on the original surface and smoothed with both methods. For the comparison of performance between the smoothing methods, we calculated the root mean squared errors (RMSE) between them, where the mean of the squared errors is taken over the surface. For the heat kernel regression, we used the bandwidth and eigenfunctions up to degree. For iterated kernel smoothing, we varied the number of iterations with correspondingly smaller bandwidth , which results in an effective bandwidth of . For diffusion smoothing, a sufficiently small step size was taken for 200 iterations resulting in bandwidth . The RMSE between the kernel regression and the iterated kernel smoothing reached up to 0.5901 (-coordinate) and did not decrease even when we increased the number of iterations (Figure 7). The RMSE between the kernel regression and diffusion smoothing was smaller than 0.0046 (-coordinates). Figure 7 (b) shows the squared differences between the two methods. For the iterated kernel smoothing, the difference is mainly localized in high curvature areas, where the Gaussian kernel used in the iterated kernel smoothing fails to approximate the heat kernel. This comparison clearly demonstrates the limitation of iterated heat kernel smoothing, which does not converge to heat diffusion. However, the heat kernel regression and diffusion smoothing gave almost identical results and there was no discernible difference (Figure 7 (c).

We also compared the performance of the three smoothing techniques at four different bandwidths . For the kernel regression, was used. For the iterated kernel smoothing and the diffusion smoothing, a fixed step size of was used with , 800, 2000, 4000 iterations. The diffusion smoothing and heat kernel smoothing gave visually identical results for bandwidths due to a sufficiently large number of iterations (Figure 8). However, the iterated kernel smoothing gave a different result. In this experiment, we replaced the original surface coordinates with smoothed ones for the final visualization. However, in the actual computation, we did not replace the original surface coordinates for the three methods. Iterated kernel smoothing compounded the discretization errors over iterations, so it did not converge to the kernel regression and diffusion smoothing. Diffusion smoothing and heat kernel smoothing share the same FEM discretization and converge to each other as the number of iterations increases.

Figure 8: Smoothed mandible surfaces using three different techniques. The x, y and z surface coordinates are treated as functional measurements on the original surface and smoothed. The proposed heat kernel smoothing is done with bandwidths, . Iterated kernel smoothing is done by iteratively approximating the heat kernel linearly with the Gaussian kernel (Chung et al., 2005). Diffusion smoothing directly solves the diffusion equation using the same FEM discretization (Chung and Taylor, 2004). Diffusion smoothing and heat kernel smoothing converge to each other as the bandwidth increases.

3.3 Simulation Studies

Since there is no known ground truth in the imaging data set we are using, it is uncertain how the proposed method will perform with real data. It is therefore necessary to perform simulation studies with ground truths. We performed two simulations with small and large SNR settings on a T-junction shaped surface (Figure 9), which was chosen because it was a surface with three different curvatures: convex, concave and almost flat regions. Note that surface smoothing methods perform differently under different curvatures. Three black signal regions of different sizes were taken as the ground truth at these regions and 60 independent functional measurements on the T-junction were simulated as , the absolute value of normal distribution with mean 0 and variance , at each mesh vertex. Value 1 was then added to the black regions in 30 of the measurements, which served as group 2, while the other 30 measurements were taken as group 1. Group 1 had distribution while group 2 had distribution in the signal regions. Larger variance corresponds to smaller SNR.

Figure 9: Simulation study I on a T-junction shaped surface where three black signal regions of different sizes are taken as the ground truth. 60 independent functional measurements on the T-junction were simulated as at each mesh vertex. We are only simulating positive numbers to better reflect the positive measurements used in the study. Value 1 was added to the black regions in 30 measurements, which served as group 2 while the other 30 measurements were taken as group 1. T-statistics are shown for these simulations (original) and three techniques with bandwidth . Heat kernel smoothing performed the best in detecting the ground truth.
Figure 10: Simulation study II on a T-junction shaped surface with the same ground truth as simulation study I (Figure 9). 60 independent functional measurements on the T-junction were simulated as at each mesh vertex. Value 1 was added to the black regions in 30 of the measurements, which served as group 2, while the remaining 30 measurements are taken as group 1. Due to the large SNR, the group means show visible group separations. All the methods detected the signal regions; however, the heat kernel smoothing and diffusion smoothing techniques were more sensitive at the large SNR.

In Study I, was used to simulate functional measurements with substantially smaller SNR. Figure 9 shows the simulation results. For iterated kernel and diffusion smoothing, we used bandwidth and 100 iterations. For smaller SNR, it is necessary to smooth with a larger bandwidth, which is determined empirically. For heat kernel smoothing, the same bandwidth and 1000 eigenfunctions were used. The same number of eigenfunctions was used throughout the study. For all three smoothing techniques, the bandwidth is the main parameter that determines performance. We then performed a two sample t-test with the RFT-based threshold of 4.90 to detect the group difference at level.

Neither the raw data nor iterated smoothing were able to correctly identify any signal region. However, heat kernel and diffusion smoothing correctly identified 94 and 91 of the signal regions respectively. In addition, heat kernel and diffusion smoothing incorrectly identified 0.26 and 0.26 of non-signal regions as signal. There are no visually discernible differences between the two methods as shown in Figure 10. The difference in performance is due to the discretization error associated with taking only 100 iterations, which disappears if we take smaller time steps. Alternatively, we can use better time-discretization schemes such as Padé-Chebyschev approximation (Patané, 2015), multi-step methods (Gottlieb et al., 2001) or higher-order Runge-Kutta schemes (Li and Alexiades, 2010).

In Study II, was used to simulate functional measurements with substantially larger SNR. Due to the large SNR, the group means showed visible group separations (Figure 10). For iterated kernel and diffusion smoothing, we used bandwidth and 100 iterations. For heat kernel smoothing, the same bandwidth and 1000 eigenfunctions were used. All the methods detected the signal regions; however, the heat kernel smoothing and diffusion smoothing techniques were more sensitive at large SNR. All the methods correctly identified the signal regions with 100 accuracy. There were no false discoveries in the raw data and iterated kernel smoothing methods. However, due to blurring effects, heat kernel and diffusion smoothing incorrectly identified 0.9 and 0.8 of non-signal regions as signal, which is negligible. For the large SNR setting, all the methods were reasonably able to detect the correct signal regions with minimal error.

In summary, in larger SNR, all three methods performed well. However, in substantially smaller SNR, the kernel regression performed best, closely followed by diffusion smoothing. Neither the raw data nor iterated kernel smoothing performed well in the low SNR setting.

4 Application: Mandible Growth Analysis

As an illustration of the proposed kernel regression technique, we analyzed mandible growth on a CT imaging data set consisting of 77 human subjects between the ages of 0 and 19 years. Subjects were divided into three age categories: 0 to 6 years (group I, 26 subjects), 7 to 12 years (group II, 20 subjects), and 13 to 19 years (group III, 31 subjects). The main biological question of interest was whether there were localized regions of growth between these different age groups. Mandible surface meshes for all subjects were constructed through the image acquisition and processing steps described in the previous section. For surface alignment, diffeomorphic surface registration was used to match mandible surfaces across subjects (Miller and Qiu, 2009; Vaillant et al., 2007; Qiu and Miller, 2008; Yang et al., 2011).

4.1 Diffeomorphic Surface Registration

Figure 11: Left: Mandible F155-12-08, which forms an initial template . All other mandibles were affine registered to F155-12-08. Middle: The superimposition of affine registered mandibles showing local misalignment. Diffeomorphic registration was then performed to warp misaligned affine transformed mandibles. Right: The average of deformation with respect to F155-12-08 provides the final population average template where statistical parametric maps were constructed.

We chose the mandible of a 12-year-old subject identified as F155-12-08, which served as the reference template in previous studies (Seo et al., 2010, 2011), as initial template and aligned the remaining 76 mandibles to the initial template affinely to remove the overall size variability. Some subjects may have larger mandibles than others, so it is necessary to remove the global size differences in localized shape modeling. From the affine transformed individual mandible surfaces , we performed an additional nonlinear surface registration to the template using the large deformation diffeomorphic metric mapping (LDDMM) framework (Miller and Qiu, 2009; Vaillant et al., 2007; Qiu and Miller, 2008; Yang et al., 2011).

In the LDDMM framework (Miller and Qiu, 2009; Vaillant et al., 2007; Qiu and Miller, 2008; Yang et al., 2011), the metric space is constructed as an orbit of surface under the group of diffeomorphic transformations , i.e., . The diffeomorphic transformations (one-to-one, smooth forward and inverse transformation) are introduced as transformations of the coordinates on the background space . The diffeomorphisms are constructed as a flow of ordinary differential equations (ODE), where follows


where denotes the identity map and are the associated velocity vector fields. The vector fields are constrained to be sufficiently smooth, so that (22) is integrable and generates diffeomorphic transformations over finite time. The smoothness is ensured by forcing to lie in a smooth vector field , which is modeled as a reproducing kernel Hilbert space with linear operator associated with norm (Dupuis et al., 1998). The group of diffeomorphisms is then the solutions of (22) with the vector fields satisfying .

Given the template surface and an individual surface , the geodesic , which lies in the manifold of diffeomorphisms and connects and , is defined as

Figure 12: Mandibles were grouped into three age cohorts: group I (ages 0 to 6 years), group II (ages 7 to 12 years) and group III (ages 13 to 19 years). Each row shows the mean group differences of the displacement: group II - group I (first row) and group III - group II (second row). The arrows are the mean displacement differences and colors indicate their lengths in mm. Longer arrows imply more mean displacement.

For our application, we employed the LDDMM approach to estimate the template among all subjects. The estimated template can be simply computed through averaging the initial velocity across all subjects (Zhong and Qiu, 2010), which is similar to the unbiased template estimation approach in Joshi et al. (2004). We then recomputed the displacement fields with respect to the initial template . We averaged the deformation fields from the initial template to individual subjects to obtain the final template . Figure 11 shows the initial and final templates. Figure 12 shows the mean displacement differences between groups I and II (top) and II and III (bottom). Each row shows the group differences of the displacement: group II - group I (first row) and group III - group II (second row). The arrows are the growth direction with arrow length being representative of mean displacement differences and colors indicating growth length in mm.

Figure 13: -statistic map showing the regions of significant growth as measured by mean displacement differences between the groups displayed in Figure 12. The kernel regression was used to smooth out surface measurements. The top row shows significant growth between groups I and II ; and bottom row between groups II and III. The thresholds 10.52 and 10.67 are considered significant at 0.01 level (corrected) for the top and bottom rows.
Figure 14: -statistic map showing the regions of significant growth as measured by mean displacement differences between the groups displayed in Figure 12. The iterated kernel smoothing with parameters and were used. The top row shows significant growth between groups I and II ; and bottom row between groups II and III. The thresholds 10.52 and 10.67 are considered significant at 0.01 level (corrected) for the top and bottom rows.

4.2 Statistical Analysis

We are interested in determining the significance of the mean displacement differences in Figure 12. Since the length measurement provides a much easier biological interpretation, we used the length of the displacement vector as a response variable. The RFT assumes the measurements to be a smooth Gaussian random field. Heat kernel smoothing on the length measurement will make the measurement smoother, more Gaussian and increase the SNR (Chung et al., 2005). Heat kernel smoothing is applied with bandwidth using eigenfunctions on the final template . The number of eigenfunctions used is more than sufficient to guarantee a relative error less than . The heat kernel smoothing of the displacement length is taken as the response variable. We constructed the random field testing the length difference between the age groups I and II, and II and III showing the regions of accelerated growth (Figure 13).

The comparison of groups I and II is based on an -field with 1 and 44 degrees of freedom. The comparison of groups II and III is based on an -field with 1 and 49 degrees of freedom. The multiple comparison corrected -statistics thresholds corresponding to and levels are respectively 8.00 and 10.52 (group II-I) and 8.00 and 10.67 (group III- II). In the -statistic map shown in Figure 13, black and red regions are considered as exhibiting growth spurts at 0.01 and 0.05 levels respectively. Our findings are consistent with previous findings of simultaneous forward and downward growth (Scott, 1976; Smartt Jr. et al., 2005; Walker and Kowalski, 1972; Lewis et al., 1982; Seo et al., 2011) and bilateral growth (Enlow and Hans, 1996).

We also performed the same statistical analysis to the iterated kernel smoothing and diffusion smoothing results. For the diffusion smoothing, 200 step sizes are used with the fixed time step , which results in the overall bandwidth . For the iterated kernel smoothing, bandwidth is split into small bandwidths. The diffusion smoothing results are similar to Figure 13 so the only iterated kernel smoothing result is shown in Figure 14. Since this is a high SNR setting, all three methods are expected to perform well and similarly. In the heat kernel regression, 25 of mesh vertices show -statistic value above 8.00 for the comparison of groups I and II (0.05 level). For the iterated kernel smoothing and diffusion smoothing, 24 and 24 of vertices are above 8.00. For the comparison of groups II and III, the numbers are 38, 36 and 36 respectively. The differences are not significant.

5 Conclusions

This study presents a novel heat kernel regression framework, where functional measurements are expanded analytically using the weighted Laplace-Beltrami eigenfunctions. The weighted eigenfunction expansion is related to isotropic heat diffusion and the diffusion wavelet transform. The method was validated and compared against exiting surface-based smoothing methods. Although the proposed kernel regression and diffusion smoothing share identical FEM discretization, the kernel regression is a parametric model, whereas diffusion smoothing is not. The flexibility of the parametric model enabled us to establish the mathematical equivalence of kernel regression, diffusion smoothing and diffusion wavelets.

The method was subsequently applied to characterize mandible growth. Based on the significant directions of growth identified in Figure 12 and 13, we quantified the regions, direction and extent of growth during the first two decades of life that contribute to the overall downward and forward growth of the mandible as described in the literature. To quantify mandibular growth using smaller age cohorts, we are currently securing additional CT images.


This work was supported by NIH Research Grant R01 DC6282 (MRI and CT Studies of the Developing Vocal Tract), from the National Institute of Deafness and other Communicative Disorders (NIDCD); core grant P-30 HD03352 to the Waisman Center from the National Institute of Child Health and Human Development (NICHHD); Clinical and Translational Science Award (CTSA) program, through the NIH National Center for Advancing Translational Sciences (NCATS) UL1TR000427, and Vilas Associate Award from University of Wisconsin-Madison. We thank Lindell R. Gentry, Mike S. Schimek, Brian J. Whyms, Reid B. Durtschi and Dongjun Chung at the University of Wisconsin-Madison for assistance with image acquisition and image segmentation. We also thank the anonymous reviewers and Yuan Wang and Jacqueline Houtman for comments on earlier versions of this paper.


  • Adler (1981) Adler, R., 1981. The Geometry of Random Fields. John Wiley Sons.
  • Andrade et al. (2001) Andrade, A., Kherif, F., Mangin, J., Worsley, K., Paradis, A., Simon, O., Dehaene, S., Le Bihan, D., Poline, J.B., 2001. Detection of fMRI activation using cortical surface mapping. Human Brain Mapping 12, 79–93.
  • Andresen et al. (2000) Andresen, P.R., Bookstein, F.L., Conradsen, K., Ersbøll, B.K., Marsh, J.L., Kreiborg, S., 2000. Surface-bounded growth modeling applied to human mandibles. IEEE Transactions on Medical Imaging 19, 1053–1063.
  • Antoine et al. (2010) Antoine, J.P., Roşca, D., Vandergheynst, P., 2010. Wavelet transform on manifolds: old and new approaches. Applied and Computational Harmonic Analysis 28, 189–202.
  • Belkin et al. (2006) Belkin, M., Niyogi, P., Sindhwani, V., 2006. Manifold regularization: A geometric framework for learning from labeled and unlabeled examples. The Journal of Machine Learning Research 7, 2399–2434.
  • Bernal-Rusiel et al. (2008) Bernal-Rusiel, J., Atienza, M., Cantero, J., 2008. Detection of focal changes in human cortical thickness: spherical wavelets versus gaussian smoothing. NeuroImage 41, 1278–1292.
  • Bronstein and Kokkinos (2010) Bronstein, M.M., Kokkinos, I., 2010. Scale-invariant heat kernel signatures for non-rigid shape recognition, in: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1704–1711.
  • Cachia et al. (2003a) Cachia, A., Mangin, J.F., Riviére, D., Kherif, F., Boddaert, N., Andrade, A., Papadopoulos-Orfanos, D., Poline, J.B., Bloch, I., Zilbovicius, M., Sonigo, P., Brunelle, F., Régis, J., 2003a. A primal sketch of the cortex mean curvature: a morphogenesis based approach to study the variability of the folding patterns. IEEE Transactions on Medical Imaging 22, 754–765.
  • Cachia et al. (2003b) Cachia, A., Mangin, J.F., Riviére, D., Papadopoulos-Orfanos, D., Kherif, F., Bloch, I., Régis, J., 2003b. A generic framework for parcellation of the cortical surface into gyri using geodesic Voronoï diagrams. Image Analysis 7, 403–416.
  • Cao and Worsley (2001) Cao, J., Worsley, K., 2001. Applications of random fields in human brain mapping. Spatial Statistics: Methodological Aspects and Applications 159, 170–182.
  • Chung (2001) Chung, M., 2001. Statistical Morphometry in Neuroanatomy. Ph.D. Thesis, McGill University. http://www.stat.wisc.edu/~mchung/papers/thesis.pdf.
  • Chung et al. (2008) Chung, M., Hartley, R., Dalton, K., Davidson, R., 2008. Encoding cortical surface by spherical harmonics. Statistica Sinica 18, 1269–1291.
  • Chung et al. (2005) Chung, M., Robbins, S., Evans, A., 2005. Unified statistical approach to cortical thickness analysis. Information Processing in Medical Imaging (IPMI), Lecture Notes in Computer Science 3565, 627–638.
  • Chung and Taylor (2004) Chung, M., Taylor, J., 2004. Diffusion smoothing on brain surface via finite element method, in: Proceedings of IEEE International Symposium on Biomedical Imaging (ISBI), pp. 432–435.
  • Chung et al. (2001) Chung, M., Worsley, K., Taylor, J., Ramsay, J., Robbins, S., Evans, A., 2001. Diffusion smoothing on the cortical surface. NeuroImage 13, S95.
  • Coifman and Maggioni (2006) Coifman, R., Maggioni, M., 2006. Diffusion wavelets. Applied and Computational Harmonic Analysis 21, 53–94.
  • Dupuis et al. (1998) Dupuis, P., Grenander, U., Miller, M.I., 1998. Variational problems on flows of diffeomorphisms for image matching. Quarterly of Applied Math. 56, 587–600.
  • Enlow and Hans (1996) Enlow, D.H., Hans, M.G., 1996. Essentials of Facial Growth. W.B. Saunders Company, Philadelphia.
  • Evans (1998) Evans, L., 1998. Partial differential equations. American Mathematical Society.
  • Gaser et al. (2006) Gaser, C., Luders, E., Thompson, P., Lee, A., Dutton, R., Geaga, J., Hayashi, K., Bellugi, U., Galaburda, A., Korenberg, J., Mills, D., Toga, A., Reiss, A., 2006. Increased local gyrification mapped in williams syndrome. NeuroImage 33, 46–54.
  • de Goes et al. (2008) de Goes, F., Goldenstein, S., Velho, L., 2008. A hierarchical segmentation of articulated bodies. Computer Graphics Forum 27, 1349–1356.
  • Gottlieb et al. (2001) Gottlieb, S., Shu, C.W., Tadmor, E., 2001. Strong stability-preserving high-order time discretization methods. SIAM review 43, 89–112.
  • Guskov and Wood (2001) Guskov, I., Wood, Z., 2001. Topological noise removal, in: Graphics Interface, pp. 19–26.
  • Hagler Jr. et al. (2006) Hagler Jr., D., Saygin, A., Sereno, M., 2006. Smoothing and cluster thresholding for cortical surface-based group analysis of fMRI data. NeuroImage 33, 1093–1103.
  • Hammond et al. (2011) Hammond, D., Vandergheynst, P., Gribonval, R., 2011. Wavelets on graphs via spectral graph theory. Applied and Computational Harmonic Analysis 30, 129–150.
  • Han et al. (2007) Han, J., Kim, J., Chung, C., Park, K., 2007. Evaluation of smoothing in an iterative lp-norm minimization algorithm for surface-based source localization of meg. Physics in Medicine and Biology 52, 4791–4803.
  • Han et al. (2006) Han, X., Jovicich, J., Salat, D., van der Kouwe, A., Quinn, B., Czanner, S., Busa, E., Pacheco, J., Albert, M., Killiany, R., et al., 2006. Reliability of MRI-derived measurements of human cerebral cortical thickness: The effects of field strength, scanner upgrade and manufacturer. NeuroImage 32, 180–194.
  • Hayasaka et al. (2004) Hayasaka, S., Phan, K.L., Liberzon, I., Worsley, K.J., Nichols, T.E., 2004. Nonstationary cluster-size inference with random field and permutation methods. Neuroimage 22, 676–687.
  • Hernandez et al. (2006) Hernandez, V., Roman, J.E., Tomas, A., Vidal, V., 2006. A Survey of Software for Sparse Eigenvalue Problems. Technical Report STR-6. Universidad Politécnica de Valencia. http://www.grycap.upv.es/slepc.
  • Jo et al. (2007) Jo, H., Lee, J.M., Kim, J.H., Shin, Y.W., Kim, I.Y., Kwon, J., Kim, S., 2007. Spatial accuracy of fmri activation influenced by volume- and surface-based spatial smoothing techniques. NeuroImage 34, 550–564.
  • Joshi et al. (2009) Joshi, A., Shattuck, D.W., Thompson, P.M., Leahy, R.M., 2009. A parameterization-based numerical method for isotropic and anisotropic diffusion smoothing on non-flat surfaces. IEEE Transactions on Image Processing 18, 1358–1365.
  • Joshi et al. (2004) Joshi, S.C., Davis, B., Jomier, M., Gerig, G., 2004. Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage 23, 151–160.
  • Kilner and Friston (2010) Kilner, J., Friston, K., 2010. Topological inference for eeg and meg. The Annals of Applied Statistics 4, 1272–1290.
  • Kim et al. (2011) Kim, S.G., Chung, M., Seo, S., Schaefer, S., van Reekum, C., Davidson, R., 2011. Heat kernel smoothing via Laplace-Beltrami eigenfunctions and its application to subcortical structure modeling, in: Pacific-Rim Symposium on Image and Video Technology (PSIVT). Lecture Notes in Computer Science (LNCS), pp. 36–47.
  • Kim et al. (2012) Kim, W., Pachauri, D., Hatt, C., Chung, M., Johnson, S., Singh, V., 2012. Wavelet based multi-scale shape features on arbitrary surfaces for cortical thickness discrimination, in: Advances in Neural Information Processing Systems, pp. 1250–1258.
  • Lai et al. (2010) Lai, Z., Hu, J., Liu, C., Taimouri, V., Pai, D., Zhu, J., Xu, J., Hua, J., 2010. Intra-patient supine-prone colon registration in CT colonography using shape spectrum, in: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2010, pp. 332–339.
  • Lehoucq et al. (1998) Lehoucq, R., Sorensen, D., Yang, C., 1998. ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM Publications, Philadelphia.
  • Lerch and Evans (2005) Lerch, J.P., Evans, A., 2005. Cortical thickness analysis examined through power analysis and a population simulation. NeuroImage 24, 163–173.
  • Lévy (2006) Lévy, B., 2006. Laplace-Beltrami eigenfunctions towards an algorithm that ”understands” geometry, in: IEEE International Conference on Shape Modeling and Applications, p. 13.
  • Lewis et al. (1982) Lewis, A., Roche, A., Wagner, B., 1982. Growth of the mandible during pubescence. The Angle Orthodontist 52, 325–342.
  • Li and Alexiades (2010) Li, C., Alexiades, V., 2010. Time stepping for the cable equation, Part 1: Serial performance. Proceedings of Neural, Parallel & Scientific Computations 4, 241–246.
  • Loubele et al. (2006) Loubele, M., Maes, F., Schutyser, F., Marchal, G., Jacobs, R., Suetens, P., 2006. Assessment of bone segmentation quality of cone-beam CT versus multislice spiral CT: a pilot study. Oral Surgery, Oral Medicine, Oral Pathology, Oral Radiology, and Endodontology 102, 225–234.
  • Luders et al. (2006a) Luders, E., Narr, K., Thompson, P., Rex, D., Woods, R., DeLuca, H., Jancke, L., Toga, A., 2006a. Gender effects on cortical thickness and the influence of scaling. Human Brain Mapping 27, 314–324.
  • Luders et al. (2006b) Luders, E., Thompson, P., Narr, K., Toga, A., Jancke, L., Gaser, C., 2006b. A curvature-based approach to estimate local gyrification on the cortical surface. NeuroImage 29, 1224–1230.
  • Mahadevan and Maggioni (2006) Mahadevan, S., Maggioni, M., 2006. Value function approximation with diffusion wavelets and laplacian eigenfunctions. Advances in neural information processing systems 18, 843.
  • Malladi and Ravve (2002) Malladi, R., Ravve, I., 2002. Fast difference schemes for edge enhancing Beltrami flow, in: Proceedings of Computer Vision-ECCV, Lecture Notes in Computer Science (LNCS), pp. 343–357.
  • Miller and Qiu (2009) Miller, M., Qiu, A., 2009. The emerging discipline of computational functional anatomy. Neuroimage 45, S16–S39.
  • Nain et al. (2007) Nain, D., Styner, M., Niethammer, M., Levitt, J., Shenton, M., Gerig, G., Bobick, A., Tannenbaum, A., 2007. Statistical shape analysis of brain structures using spherical wavelets, in: IEEE Symposium on Biomedical Imaging ISBI.
  • Niethammer et al. (2007) Niethammer, M., Reuter, M., Wolter, F., Bouix, S., Peinecke, N., Koo, M., Shenton, M., 2007. Global Medical Shape Analysis Using the Laplace-Beltrami Spectrum. Lecture Notes in Computer Science 4791, 850.
  • Nilsson et al. (2007) Nilsson, J., Sha, F., Jordan, M., 2007. Regression on manifolds using kernel dimension reduction, in: Proceedings of the 24th international conference on Machine learning, ACM. pp. 697–704.
  • Patané (2015) Patané, G., 2015. Volumetric heat kernel: Padé-Chebyshev approximation, convergence, and computation. Computers & Graphics 46, 64–71.
  • Perona and Malik (1990) Perona, P., Malik, J., 1990. Scale-space and edge detection using anisotropic diffusion. IEEE Trans. Pattern Analysis and Machine Intelligence 12, 629–639.
  • Qiu et al. (2006) Qiu, A., Bitouk, D., Miller, M., 2006. Smooth functional and structural maps on the neocortex via orthonormal bases of the laplace-beltrami operator. IEEE Transactions on Medical Imaging 25, 1296–1396.
  • Qiu and Miller (2008) Qiu, A., Miller, M., 2008. Multi-structure network shape analysis via normal surface momentum maps. NeuroImage 42, 1430–1438.
  • Reuter (2010) Reuter, M., 2010. Hierarchical shape segmentation and registration via topological features of Laplace-Beltrami eigenfunctions. International Journal of Computer Vision 89, 287–308.
  • Reuter et al. (2006) Reuter, M., Wolter, F.E., Peinecke, N., 2006. Laplace-Beltrami spectra as Ôshape-dnaÕ of surfaces and solids. Computer-Aided Design 38, 342–366.
  • Reuter et al. (2009) Reuter, M., Wolter, F.E., Shenton, M., Niethammer, M., 2009. Laplace-Beltrami eigenvalues and topological features of eigenfunctions for statistical shape analysis. Computer-Aided Design 41, 739–755.
  • Rosenberg (1997) Rosenberg, S., 1997. The Laplacian on a Riemannian Manifold. Cambridge University Press.
  • Rustamov (2007) Rustamov, R., 2007. Laplace-Beltrami eigenfunctions for deformation invariant shape representation, in: Proceedings of the fifth Eurographics symposium on Geometry processing, p. 233.
  • Schölkopf and Smola (2002) Schölkopf, B., Smola, A.J., 2002. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press.
  • Scott (1976) Scott, J.H., 1976. Dento-facial Development and Growth. Pergamon Press, Oxford.
  • Seo and Chung (2011) Seo, S., Chung, M., 2011. Laplace-Beltrami eigenfunction expansion of cortical manifolds, in: Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium on, IEEE. pp. 372–375.
  • Seo et al. (2011) Seo, S., Chung, M., Voperian, H., 2011. Mandible shape modeling using the second eigenfunction of the laplace-beltrami operator, in: Proceedings of SPIE Medical Imaging, p. 79620Z.
  • Seo et al. (2010) Seo, S., Chung, M., Vorperian, H., 2010. Heat kernel smothing using laplace-beltrami eigenfunctions, in: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2010, pp. 505–512.
  • Shawe-Taylor and Cristianini (2004) Shawe-Taylor, J., Cristianini, N., 2004. Kernel methods for pattern analysis. Cambridge University Press.
  • Shen et al. (2004) Shen, L., Ford, J., Makedon, F., Saykin, A., 2004. Surface-based approach for classification of 3d neuroanatomical structures. Intelligent Data Analysis 8, 519–542.
  • Shen et al. (2006) Shen, L., Saykin, A., Chung, M., Huang, H., Ford, J., Makedon, F., McHugh, T., Rhodes, C., 2006. Morphometric analysis of genetic variation in hippocampal shape in mild cognitive impairment: Role of an il-6 promoter polymorphism, in: Life Science Society Computational Systems Bioinformatics Conference.
  • Shi et al. (2009) Shi, Y., Dinov, I., Toga, A.W., 2009. Cortical shape analysis in the Laplace-Beltrami feature space, in: 12th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI 2009), pp. 208–215.
  • Shi et al. (2008) Shi, Y., Lai, R., Kern, K., Sicotte, N., Dinov, I., Toga, A.W., 2008. Harmonic surface mapping with Laplace-Beltrami eigenmaps, in: 11th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI 2008), pp. 147–154.
  • Smartt Jr. et al. (2005) Smartt Jr., J.M., Low, D.W., Bartlett, S.P., 2005. The pediatric mandible: I. A primer on growth and development. Plast. Reconstr. Surg. 116, 14e–23e.
  • Sochen et al. (1998) Sochen, N., Kimmel, R., Malladi, R., 1998. A general framework for low level vision. IEEE Transactions on Image Processing 7, 310–318.
  • Steinke and Hein (2008) Steinke, F., Hein, M., 2008. Non-parametric regression between manifolds. Advances in Neural Information Processing Systems .
  • Stratemann et al. (2010) Stratemann, S.A., Huang, J.C., Maki, K., Hatcher, D.C., Miller, A.J., 2010. Evaluating the mandible with cone-beam computed tomography. American Journal of Orthodontics and Dentofacial Orthopedics 137, S58–S70.
  • Styner et al. (2006) Styner, M., Oguz, I., Xu, S., Brechbuhler, C., Pantazis, D., Levitt, J., Shenton, M., Gerig, G., 2006. Framework for the statistical shape analysis of brain structures using spharm-pdm, in: Insight Journal, Special Edition on the Open Science Workshop at MICCAI.
  • Sun et al. (2009) Sun, J., Ovsjanikov, M., Guibas, L.J., 2009. A concise and provably informative multi-scale signature based on heat diffusion. Comput. Graph. Forum 28, 1383–1392.
  • Tang et al. (1999) Tang, B., Sapiro, G., Caselles, V., 1999. Direction diffusion, in: The Proceedings of the Seventh IEEE International Conference on Computer Vision, pp. 2:1245–1252.
  • Tasdizen et al. (2006) Tasdizen, T., Whitaker, R., Burchard, P., Osher, S., 2006. Geometric surface smoothing via anisotropic diffusion of normals, in: Geometric Modeling and Processing, pp. 687–693.
  • Taubin (2000) Taubin, G., 2000. Geometric Signal Processing on Polygonal Meshes, in: EUROGRAPHICS.
  • Taylor and Worsley (2007) Taylor, J., Worsley, K., 2007. Detecting sparse signals in random fields, with an application to brain mapping. Journal of the American Statistical Association 102, 913–928.
  • Vaillant et al. (2007) Vaillant, M., Qiu, A., Glaunès, J., Miller, M., 2007. Diffeomorphic metric surface mapping in subregion of the superior temporal gyrus. NeuroImage 34, 1149–1159.
  • Vallet and Lévy (2008) Vallet, B., Lévy, B., 2008. Spectral geometry processing with manifold harmonics. Computer Graphics Forum 27, 251–260.
  • Van Den Boomgaard and Van Balen (1992) Van Den Boomgaard, R., Van Balen, R., 1992. Methods for fast morphological image transforms using bitmapped binary images. CVGIP: Graphical Models and Image Processing 54, 252–258.
  • Walker and Kowalski (1972) Walker, G.F., Kowalski, C.J., 1972. On the growth of the mandible. American Journal of Physical Anthropology 36, 111–118.
  • Wang (1997) Wang, F.Y., 1997. Sharp explict lower bounds of heat kernels. Annals of Probability 24, 1995–2006.
  • Wang et al. (2010) Wang, Y., Zhang, J., Gutman, B., Chan, T., Becker, J., Aizenstein, H., Lopez, O., Tamburo, R., Toga, A., Thompson, P., 2010. Multivariate tensor-based morphometry on surfaces: Application to mapping ventricular abnormalities in HIV/AIDS. NeuroImage 49, 2141–2157.
  • Wardetzky (2008) Wardetzky, M., 2008. Convergence of the cotangent formula: An overview, in: Bobenko, A., Sullivan, J., Schröder, P., Ziegler, G. (Eds.), Discrete Differential Geometry. Birkhäuser Basel. volume 38 of Oberwolfach Seminars, pp. 275–286.
  • Wood et al. (2004) Wood, Z., Hoppe, H., Desbrun, M., Schröder, P., 2004. Removing excess topology from isosurfaces. ACM Transactions on Graphics (TOG) 23, 190–208.
  • Worlsey et al. (1995) Worlsey, K., Poline, J.B., Vandal, A., Friston, K., 1995. Test for distributed, non-focal brain activations. NeuroImage 2, 173–181.
  • Worsley (2003) Worsley, K., 2003. Detecting activation in fMRI data. Statistical Methods in Medical Research. 12, 401–418.
  • Worsley et al. (1999) Worsley, K., Andermann, M., Koulis, T., MacDonald, D., Evans, A., 1999. Detecting changes in nonisotropic images. Human Brain Mapping 8, 98–101.
  • Worsley et al. (1998) Worsley, K., Cao, J., Paus, T., Petrides, M., Evans, A., 1998. Applications of random field theory to functional connectivity. Human Brain Mapping 6, 364–7.
  • Worsley et al. (2004) Worsley, K., Taylor, J., Tomaiuolo, F., Lerch, J., 2004. Unified univariate and multivariate random field theory. NeuroImage 23, S189–195.
  • Yang et al. (2011) Yang, X., Goh, A., Qiu, A., 2011. Locally linear diffeomorphic metric embedding (LLDME) for surface-based anatomical shape modeling. NeuroImage 56, 149–161.
  • Yger and Rakotomamonjy (2011) Yger, F., Rakotomamonjy, A., 2011. Wavelet kernel learning. Pattern Recognition 44, 2614–2629.
  • Yotter et al. (2009) Yotter, R.A., Dahnke, R., Gaser, C., 2009. Topological correction of brain surface meshes using spherical harmonics, in: Medical Image Computing and Computer-Assisted Intervention (MICCAI), Springer. pp. 125–132.
  • Yushkevich et al. (2008) Yushkevich, P., Zhang, H., Simon, T., Gee, J., 2008. Structure-specific statistical mapping of white matter tracts. NeuroImage 41, 448–461.
  • Zhang et al. (2007) Zhang, H., van Kaick, O., Dyer, R., 2007. Spectral methods for mesh processing and analysis, in: EUROGRAPHICS, pp. 1–22.
  • Zhang et al. (2010) Zhang, H., van Kaick, O., Dyer, R., 2010. Spectral mesh processing. Computer Graphics Forum 29, 1865–1894.
  • Zhong and Qiu (2010) Zhong, J., Qiu, A., 2010. Multi-manifold diffeomorphic metric mapping for aligning cortical hemispheric surfaces. Neuroimage 49, 355–365.
  • Zhu et al. (2007) Zhu, H., Ibrahim, J., Tang, N., Rowe, D., Hao, X., Bansal, R., Peterson, B., 2007. A statistical analysis of brain morphology using wild bootstrapping. IEEE Transactions on Medical Imaging 26, 954–966.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description