Introduction
Gaussian kernel smoothing has been widely used in 3D brain images to increase the signaltonoise ratio. The Gaussian kernel weights observations according to their Euclidean distance. When the observations lie on a convoluted brain surface; however, it is more natural to assign the weight based on the geodesic distance along the surface [13]. On a 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. With this motivation in mind, we construct the kernel of a heat equation on manifolds that should be isotropic in the local conformal coordinates, and develop a framework for heat kernel smoothing on manifolds.The codes below have been tested under Matlab 7.5. If you are using the Matlab codes/sample data for your publication, please reference [1] or [2] in References.
Gaussian kernel smoothing
Heat
kernel smoothing generalizes Gaussian
kernel smoothing. Here we present a 2D
version of Gaussian kernel smoothing as
an example. The original binary image toykey.tif
is
contaminated with Gaussian white noise N(0,2^2).
signal= imread('toykey.tif');
signal=(double(signal)219)/36; %
makes the image into real numbers between 0
and 1.
figure; imagesc(signal);
colormap('bone'); colorbar
noise=
normrnd(0, 2, 596, 368);
f = signal +
noise;
figure;
imagesc(f); colormap('bone');
caxis([0 1]); colorbar
The noisy
image went though Gaussian kernel
smoothing with 10 and 20 FWHM to recover the
signal (Figure 1).
F = gaussblur(f,100);
figure;
imagesc(F); colormap('bone');
caxis([0 1]); colorbar
Heat kernel smoothing on brain surfaces
The reason that the code is short and simple is
that it has been implemented as iterative kernel
smoothing with very small bandwidth. For small
bandwidth, a heat kernel converges to a Gaussian
kernel. For details on heat kernel smoothing,
please read [1] or [2].
For MNI format, use the following lines to load a mesh into MATLAB. For this, you need to download the following codes:

hk_smooth.m heat kernel smoothing.

mni_getmesh.m load Montreal neurological institute (MNI) mesh data into MATLAB.

thickness.data sample cortical thickness corresponding to the outer surface below.

outersurface.obj sample our cortical surface mesh.
 [tri,coord,nbr,normal]=mni_getmesh('outersurface.obj');
 surf.vertices
= coord;
 surf.faces=tri;
nbr (40962x6) gives neighboring node indices. Each node has up to 6 neighboring nodes.
normal (3x40962) gives the normal vectors of the cortex.
tri (81920x3) gives three node indices that forms each of 81920 triangles.
For
FreeSurfer format, use the following code to
load a mesh into MATLAB. FreeSurfer mesh data into
MATLAB. The FreeSurfer package produces
cortical meshes with different mesh topology
for different subjects. The FreeSufer mesh
format starts the face indexing from zero
while the MATLAB can not handle zero
indexing so we need to inrease the indexing
by one.
mesh='lh.pial'
[surf.vertices, surf.faces] = read_surf(mesh);
surf.faces=surf.faces+1;
figure_wire(surf,'white','white');
To perform heat kernel smoothing using function hk_smooth.m, we generate 134506 uniform random numbers between 2 and 6mm. This simulated cortical thickness is smoothed out using heat kernel smoothing with bandwidth 1 and 20 iterations (Figure 2).
input=random('unif',2,6,[134506,1]);
figure_trimesh(surf,input,'rwb')
output=hk_smooth(input,surf,1,20);
figure_trimesh(surf,output,'rwb')
Heat kernel smoothing of multiple signals
As an application of heat kernel
smoothing, we show how to smooth a closed
surface topologically equivalent to a sphere.
hk_smooth2 is a modified version of hk_smooth
takes vector data in the first argument. Now we
take three coordinates in the form [x y z]
simultaneously.
sigma=1;
n_smooth=10;
ss.vertices=hk_smooth2(surf.vertices,surf,
sigma, n_smooth);
figure; figure_patch(surf,[0.74
0.71 0.61],0.7);
Full width at half maximum (FWHM)
The relationship between heat kernel smoothing
bandwidth and FWHM is as follows.
Kernel smoothing with Gaussian kernel K_tau =
c*exp( x^2/(2*tau^2)) corresponds to
FWHM of 2*sqrt(2*log 2) *tau. However,
in hk_smooth.m,
a slightly different Gaussian kernel form is used.
K=inline('exp(x/(4*sigma))/sum(exp(x/(4*sigma)))');
Validation against the ground truth
The performance is validated against the LaplaceBeltrami (LB) eigenfunction approach (Chung et al. 2015). In terms of performance, heat kernel smoothing is as good as the LBeigenfunction approach. You won't even be able to distinguish the resulting statistical map with naked eyes.
References
[1] 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:12561265[2] 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:627638. SpringerVerlag.
[3] Chung, M.K. 2004. Heat kernel smoothing and its application to cortical manifolds. University of WisconsinMadison, Department of Statistics, Technical Report 1090.