Heat Kernel Smoothing on Arbitrary Manifolds via Laplace-Beltrami  Eigenfunctions 

(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

This is a sample data we will use as an illustration. The left and right hippocampus surfaces are saved as triangular mesh formats and displayed in Figure 1.
 

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

To compute the eigenfunctions of the Laplace-Beltrami operator, we need to discretize using the finite element method (FEM). The FEM descretization of the Laplace-Beltrami operator was originally given in my PhD thesis in 2001 [1]. Our cotan implmentation follows notation and implmentation from [1], [2] and Qiu et at. (2006, IEEE Transactions on Medical Imaging), which basically followed the derivation given in [1] and [2]. The MATLAB function FEM produces sparse matrices A and C from the left hippocampus surface. If you are using the code, please reference [2].
 
[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);


Figure 2. Heat kernel smoothing of left hippocampus (Figure 1) with sigma=0.2 and 2000 eigenfunctions.

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

L=[300 300];

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

  1. Chung, M.K. 2001. Statistical Morphometry in Neuroanatomy, PhD Thesis, McGill University.
  2. Chung, M.K., Taylor, J. 2004. Diffusion Smoothing on Brain Surface via Finite Element Method,  IEEE International Symposium on Biomedical Imaging (ISBI). 562.
  3. 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.
  4. 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.
  5. 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.
  6. 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.
  7. Seo, S., Chung, M.K. 2011. Laplace-Beltrami eigenfunction expansion of cortical manifolds. IEEE International Symposium on Biomedical Imaging (ISBI).