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 [1] or [2] in References.

177. 
     

Figure 1. The original binary image was contaminated with Gaussian white noise N(0,2^2). Gaussian kernel smoothing with 10 and 20mm FWHM were performed to recover the original image.

 

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= imread('toy-key.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;

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

177. 
     

Figure 2. Heat kernel smoothing of real (top) and simulated (bottom) data with different number of iterations. The simulation detail is given in Chung et al (2005).

 

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

[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:1256-1265

[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:627-638. Springer-Verlag.

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