Kernel Smoothing on Arbitrary Manifolds via
(c) Moo K. Chung (wrote
the original code)
Department of Biostastics and Medical Infomatics
Waisman Laboratory for Brain Imaging and Behavior
University of Wisconsin-Madison
Florian Bachmann (drastically speed up the original code by
Institut für Geophysik und Geoinformatik
TU Bergakademie Freiberg
Seong Ho Seo
(corrected bugs in the original code)
Department of Brain
and Cognitive Sciences
Seoul National University
April 24, 2010; July 16, 2011
The concept of heat kernel
smoothing along an arbitary manifold has been first
introduced in 2005 . The Gaussian kernel weights
observations according to their Euclidean distance. When
the observations lie on a convoluted brain surface and
arbitary manifolds; however, it is more natural to assign
the weight based on the geodesic distance along the
manifolds. On the curved manifold, a straight line between
two points is not the shortest distance so one may
incorrectly assign less weights to closer observations.
Therefore, smoothing data residing on manifolds requires
constructing a kernel that is isotropic along the geodesic
The original implemenation
given in  was fomulated as iterated kernel
convolutions that approximate Gaussian kernel smoothing in
the tangent space locally. The MATLAB code is given in http://www.stat.wisc.edu/~mchung/softwares/hk/hk.html.
To correct for the confounding numerical error over each
iteration, we introduced a new smoothing framework that
uses the eigenfunctions of the Laplace Beltrami operator.
This implementation of heat kernel smoothing probably
solves an isosotropic heat diffusion on the manifolds mos
accruately without the problem of divergence. Our new
method performs surface data smoothing by
constructing the series expansion of the eigenfunctions of
the Laplace-Beltrami operator . This new analytic
framework improves upon the previous iterated kernel
smoothing fomulation   with improved numerical
accuracy and stability. The
codes have been tested on an old Mac computer (intel processor) with 4GB
memory and MATLAB 7.5, and 2.66Ghz Quad-Core
Intel Xeon Mac computer with 32GB and MATLAB 7.9. If you are using the Matlab codes/sample
data below for your publication, please reference 
or , which introdued the method given here. This is
the only publically available MATLAB code out there.
This is a sample
data we will use as an illustration. The left and right
hippocampus surfaces are saved as triangular mesh formats and
displayed in Figure 1.
The left and right hippocampi is of the format:
vertices: [2338x3 double]
faces: [4672x3 double]
vertices: [2312x3 double]
faces: [4620x3 double]
Figure 1. Left hippocampus
showing possible mesh noise.
of Laplace-Beltrami Operator
To compute the eigenfunctions of the Laplace-Beltrami
operator, we need to discretize using the finite
element method (FEM). The FEM descretization of the
Laplace-Beltrami operator was originally given in my PhD thesis
in 2001 . Our cotan implmentation follows notation and
implmentation from ,  and Qiu et at. (2006, IEEE
Transactions on Medical Imaging), which basically followed the
derivation given in  and . The MATLAB function FEM produces sparse
matrices A and C from the left hippocampus surface. If you are using the code,
please reference .
2010; July 1, 2011
[A, C] =FEM(hippoleft);
The generalized eigenvalue problem is then computed by
using the sparse general eigenvalue problem routine in MATLAB.
D] = eigs(C,A,999,'sm');
It takes 104 seconds for computing 999 eigenfunctions in old MacBook
pro laptop while takes 45seconds in Quad-Core Intel Xeon Mac
computer. For a larger surface meshes, you need huge memory. In
MATLAB 7.9 with 2.66Ghz Quad-Core Intel Xeon Mac computer with 32GB
memory, we can construct 1000 eigenfunctions for MNI cortical
surface mesh formats (40962 vertices) in about 8 mininutes. You only
need to construct the eigenfunctions for the cortical template so
the computation burden is not severe. We haven't tried FreeSurfer
format which will likely require more memory than 32GB.
April 22, 2010;
July 16, 2011
of heat kernel smoothing in
manifolds was originally given in  . The original
implemenation was based on the iterated linear approximation
of the heat kernel using a Gaussian kernel in the tangent
space, which componds error when the number of iterations
increase. FreeSurfer package is also based on a similar
iterated Gaussian kernel smoothing. Our new
Laplace-Beltraim eigenfunction approach uses the
eigenfunctions of the Laplace-Beltrami operator in
representing the heat kernel as a series expansion involing
the eigenfunctions  .
left hippocampus surface, with sigma = 0.2 and 2000
eigenfunctions, we run
= lb_smooth(,hippoleft, 0.2,
500, V, D);
first argument is where input signal should go but since the
surface coordinates themselves are signal, we do not need to put
Figure 2. Heat kernel smoothing of left hippocampus
(Figure 1) with sigma=0.2 and 2000 eigenfunctions.
We can also smooth signal on the hippocampus surface. Supppose
signal is simply the y-coordinate function plus Gaussian white
noise N(0,10). Then we put the signal in the first argument:
signal = signal + normrnd(0,10,2338,1);
smoothed = lb_smooth(signal,hippoleft, 0.2, 500, V, D);
April 22, 2010
The validation of our proposed analytic framework is done on a
unit sphere using the spherical harmonics. The technical detail is
hown in . Here we show that our the eigenvalues obtained from
the FEM-discretization converges to l*(l+1) for l=0,1, 2,
..... on the unit sphere as the mesh resolution increases.
Figure 4. Icosahedron which ahs 20 triangles, 12 vertices and 30
edges. The triangle subdivision of the icosahedron increases
the number of triangles by a factor of 4. At each subdivision, the
vertices are projected onto the sphere. At the subsequent level of
refinement, we have (42, 80), (162, 320), (642, 1280), (2562,
5120), (10242, 20480), (40962, 81920) faces and vertices.
The following code generates the spherical mesh with 10242
vertices. The first 30 eigenvalues are then computed. sphere_tri
additional codes sphere_project
and mesh_refine_tri4 to run.
As expected, we obtin 0, 1*2, 2*3, 3*4, .... within numerical
kernel Smoothing in 2D images: analytic construction of
November 22, 2011
will show you how to construct the scale-space representation of
2D images analytically via heat kernel smoothing. The
Laplace-Beltrami eigenfunctions in the Euclidean space is simply
given in terms of the product of sine and cosine functions. We
first estimate the Fourier coefficients corresponding to the basis
functions. Then we only need to change the bandwith of haet kernel
in the computation to have another scale-space. This is an
extremely powerful framework not often used since it used to be
rectangle of size of size 300 x 300, eigenfunction_box will
construct the eigenfunction of degree l in the box size L (Figure
figure; imagesc(psi); colorbar;
Figure 5. The first six eigenfunctions in a rectangular domain. The
first eigenfunction is a constant.
Let us perform eigenfunction expansion of an image using the
constucted eigenfunctions. The deform.tif image (Figure 6 left)
shows the part of T1 weighted MRI showing the deformation field
(blue arrow) around the hipocampus (enlongated yellow structure).
If we perform the Fourier series expansion using the 10000
eigenfunctions, we obtain the Figure 6 (middle). With 90000 basis
expansion, the Fourier series expansion is almost like the original
l=[j k]; % degree of basis
coeff=sum(sum(psi.*A))/300^2; %Fourier coefficient
fourier(j+1,k+1) = coeff;
expansion=expansion + coeff*psi; %Fourier expansion
Left: original deform.tif image. Middle: 10000 basis expansion.
Left: 90000 basis expansion.
Once we obtained the Fourier coefficient, heat kernel smoothing is
done by simply changing the bandwith.
lambda= (j*pi/L(1))^2 + (k*pi/L(2))^2; %eigenvalue
coeff = fourier(j+1,k+1);
expansion=expansion + exp(-sigma*lambda)*coeff*psi; %heat kernel
Figure 7. Heat kernel smoothing with 40000 basis and different
April 22, 2010
- Chung, M.K. 2001. Statistical
Morphometry in Neuroanatomy, PhD Thesis, McGill
- Chung, M.K., Taylor, J. 2004. Diffusion
via Finite Element Method, IEEE International
Symposium on Biomedical Imaging (ISBI). 562.
- Chung, M.K., Robbins,S., Dalton, K.M.,
Davidson, Alexander, A.L., R.J., Evans, A.C. 2005. Cortical
via heat kernel smoothing. NeuroImage 25:1256-1265.
- Chung, M.K., Robbins, S., Evans, A.C. 2005. Unified
thickness analysis. Information Processing in Medical
Imaging (IPMI). Lecture Notes in Computer Science (LNCS)
- Chung, M.K., Dalton, K.M., Shen, L., L., Evans,
A.C., Davidson, R.J. 2007. Weighted
its application to quantifying the amount of gray matter.
IEEE Transactions on Medical Imaging. 26:566-581.
- Seo, S., Chung, M.K., Vorperian, H. K. 2010. Heat
eigenfunctions. 13th International Conference on Medical
Image Computing and Computer Assisted Intervention (MICCAI).
Lecture Notes in Computer Science (LNCS). 6363:505-512.
- Seo, S., Chung, M.K. 2011. Laplace-Beltrami
expansion of cortical manifolds. IEEE International
Symposium on Biomedical Imaging (ISBI).