Heat
Kernel Smoothing on Arbitrary Manifolds via
Laplace-Beltrami Eigenfunctions
(c) Moo K. Chung (wrote
the original code)
email://mkchung@wisc.edu
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
removing loops)
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
Description
April 24, 2010; July 16, 2011
The concept of heat kernel
smoothing along an arbitary manifold has been first
introduced in 2005 [3][4]. 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
curves.
The original implemenation
given in [3][4] 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 [6][7]. This new analytic
framework improves upon the previous iterated kernel
smoothing fomulation [3] [4] 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 [6]
or [7], which introdued the method given here. This is
the only publically available MATLAB code out there.
Hippocampus
Surface Data
April 22,
2010
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.
load hippocampus.mat;
The left and right hippocampi is of the format:
hippoleft =
vertices: [2338x3 double]
faces: [4672x3 double]
hipporight =
vertices: [2312x3 double]
faces: [4620x3 double]
figure;figure_patch(hippoleft,[0.7
0.7 0.6],0.5)
Figure 1. Left hippocampus
showing possible mesh noise.
Eigenfunctions
of Laplace-Beltrami Operator
April 22,
2010; July 1, 2011
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 [1]. Our cotan implmentation follows notation and
implmentation from [1], [2] and Qiu et at. (2006, IEEE
Transactions on Medical Imaging), which basically followed the
derivation given in [1] and [2]. The MATLAB function FEM produces sparse
matrices A and C from the left hippocampus surface. If you are using the code,
please reference [2].
[A, C] =FEM(hippoleft);
The generalized eigenvalue problem is then computed by
using the sparse general eigenvalue problem routine in MATLAB.
[V,
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.
Heat Kernel
Smoothing
April 22, 2010;
July 16, 2011
The concept
of heat kernel smoothing in
manifolds was originally given in [3] [4]. 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 [6] [7].
To smooth
left hippocampus surface, with sigma = 0.2 and 2000
eigenfunctions, we run
hippolefts
= lb_smooth([],hippoleft, 0.2,
500, V, D);
The
first argument is where input signal should go but since the
surface coordinates themselves are signal, we do not need to put
input signal.
figure; figure_patch(hippolefts,[0.7
0.7 0.6],0.5);
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=hippoleft.vertices(:,2)
signal = signal + normrnd(0,10,2338,1);
figure_trimesh(hippolefts,signal,'rwb')
smoothed = lb_smooth(signal,hippoleft, 0.2, 500, V, D);
figure_trimesh(hippolefts,smoothed,'rwb');
Validation
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 [7]. 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
requires two
additional codes sphere_project
and mesh_refine_tri4 to run.
sphere
= sphere_tri('ico',5,1);
figure;figure_wire(sphere, 'black',
'white');
[A,
C]= FEM(sphere);
[V,
D]= eigs(C,A,30,'sm');
sort(diag(D))
As expected, we obtin 0, 1*2, 2*3, 3*4, .... within numerical
accraucy:
-0.0000
2.0007
2.0007
2.0007
6.0044
6.0044
6.0044
6.0044
6.0044
....
Heat
kernel Smoothing in 2D images: analytic construction of
scale-spaces
November 22, 2011
We
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
computationally demanding.
In 2D
rectangle of size of size 300 x 300, eigenfunction_box will
construct the eigenfunction of degree l in the box size L (Figure
5).
L=[300
300];
for
j=0:5
l=[j j];
psi=eigenfunction_box(l,L);
figure; imagesc(psi); colorbar;
end;
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).
A =
imread('deform.tif');
A=double(A);
imagesc(A);colorbar
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
image.
expansion=zeros(300,300);
fourier=zeros(100,100);
L=[300
300]
for
j=0:99
for k=0:99
l=[j k]; % degree of basis
psi=eigenfunction_box(l,L);
coeff=sum(sum(psi.*A))/300^2; %Fourier coefficient
fourier(j+1,k+1) = coeff;
expansion=expansion + coeff*psi; %Fourier expansion
end;
end;
figure;
imagesc(expansion); colorbar
Figure 6.
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.
sigma=100;
expansion=zeros(300,300);
degree=200;
for
j=0:deg-1
for k=0:deg-1
l=[j k];
psi=eigenfunction_box(l,L);
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
expansion
end;
end;
Figure 7. Heat kernel smoothing with 40000 basis and different
bandwidths.
References
April 22, 2010
- Chung, M.K. 2001. Statistical
Morphometry in Neuroanatomy, PhD Thesis, McGill
University.
- Chung, M.K., Taylor, J. 2004. Diffusion
Smoothing
on
Brain
Surface
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
thickness
analysis
in
autism
via heat kernel smoothing. NeuroImage 25:1256-1265.
- Chung, M.K., Robbins, S., Evans, A.C. 2005. Unified
statistical
approach
to
cortical
thickness analysis. Information Processing in Medical
Imaging (IPMI). Lecture Notes in Computer Science (LNCS)
3565:627-638.
- Chung, M.K., Dalton, K.M., Shen, L., L., Evans,
A.C., Davidson, R.J. 2007. Weighted
Fourier
series
representation
and
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
kernel
smoothing
using
Laplace-Beltrami
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
eigenfunction
expansion of cortical manifolds. IEEE International
Symposium on Biomedical Imaging (ISBI).