### Introduction

Gaussian kernel smoothing has been widely used in 3D brain images to increase the signal-to-noise 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 [1-3]. 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  or  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 toy-key.tif is contaminated with Gaussian white noise N(0,2^2).

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

For MNI format, use the following lines to load a mesh into MATLAB. For this, you need to download the following codes:

[tri,coord,nbr,normal]=mni_getmesh('outersurface.obj');
surf.vertices = coord;
surf.faces=tri;

coord (3x40962) gives the x,y,z coordinates of nodes.
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.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)))');

This corresponds to FWHM of 4*sqrt(log 2* sigma). If n_smooth number of iterations with bandwidth sigma is used, FWHM is 4*sqrt( log 2 *n_smooth *sigma).

### Validation against the ground truth

The performance is validated against the Laplace-Beltrami (LB) eigenfunction approach (Chung et al. 2015). In terms of performance, heat kernel smoothing is as good as the LB-eigenfunction approach. You won't even be able to distinguish the resulting statistical map with naked eyes.

### References

 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. Springer-Verlag.

 Chung, M.K. 2004. Heat kernel smoothing and its application to cortical manifolds. University of Wisconsin-Madison, Department of Statistics, Technical Report 1090.

Update history

2005/02/05/2005 Created
2010/05/26 Smooth FreeSurfer mesh format.
2011/08/16 Gaussian kernel smoothing 2015/10/26/2015 Smooth multiple signals
2018/01/05 Cegree computation error fixed in hk_smooth2.m. Website updated