(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

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

[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);

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).

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).