Heat Kernel Smoothing on Anatomical Manifolds

Description

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]. 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. With this motivation in mind, we constructed 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 cortical manifolds.

Figure 1. Heat kernel smoothing of real (top row) and simulated (bottom row) data. Simulation is descripted in [1].

MATLAB Codes (MATLAB 6.1. implementation)

One 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. If you are using the following MATLAB codes in your research, please reference [1] or [2] in the documentation.

hk_smooth.m heat kernel smoothing.

mni_getmesh.m load polygonal mesh data into MATLAB.

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

outersurface.obj sample our cortical surface mesh.

unitsphere.mat cortical surface flattened onto a sphere. North pole is marked for visual inspection.

The following cod will load the surface into MATLAB
[tri,coord,nbr,normal]=mni_getmesh('outersurface.obj');

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.

trisurf(tri, coord(1,:),coord(2,:),coord(3,:)); % displays the outersurface
load thickness.data;
output=hk_smooth(thickness',tri,coord,nbr,1,200); % heat kernel smoothing with bandwidth 1 and 200 iterations
trisurf(tri, coord(1,:),coord(2,:),coord(3,:),output);

Documentation

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

Created 02/05/2005. Update history 01/13/2006, 11/28/2006, 09/05/2007, 09/22/2007