Mandible Surface Modeling Using Laplace-Beltrami Eigenfunctions

(c) 2013 Moo K. Chung, Houri K. Vorperian   
Vocal Tract Development Laboratory
Waisman Center
University of Wisconsin-Madison
email: mkchung@wisc.edu

Description
January 2, 2015

The mandible surface data set for all 77 subjects and MATLAB codes are stored as a single zipped file mandible77v1.zip. The zipped file also contains this web. The data set has been published in references [2][4]. The surface registration has been done by Anqi Qiu of National University of Singapore. See [4] for details on diffeomorphic surface registration. The codes have been tested on Macbook Pro 2011 laptop with 16GB memory and MATLAB 7.9 (R2009b). If you use the MATLAB codes or the full data set, please reference [4]. Send email to mkchung@wisc.edu for any inequity about MATALAB codes and data. Here is the detailed instruction for loading and analyzing the data set. The script below is given as mandible77v1.m.


Basic Data Manipulation
January 26, 2013

The full data set is stored in mandible77subjects.mat.  The variable id stores gender and age information. For instance the first subject F203-01-03 indicates the subject is female subject number 203 at age 1 year and 3 months.


>> id

ans =
    'F203-01-03'
    'F204-01-06'
    'M217-01-06'

    ...

   
We code the group variable group to be 1 for male and 0 for female. Months are converted to years.

str=cell2mat(id);
group=str(:,1)
group=[group=='M']
year = str2num(str(:,6:7));
month=str2num(str(:,9:10));
age = (12*year+month)/12;


The statistical analysis will be performed on the average surface which serves as the template. The average surface template is constructed as follows.

template.faces=faces;
template.vertices=voxelsize(44)*mean(vertices,3);
figure; figure_patch(template,[0.74 0.71 0.61],0.7);
view([90 60]); camlight;
title('Template')


template consists of 43272 triangles and 21639 vertices. To speed up the computation, we reduced the size of mesh to 10% of the original size using reducepatch.

tt=reducepatch(template,0.1);
ind=mesh_commonvertex(template, tt);
p= vertices(ind,:,:);


77 subjects are grouped into three groups: age below 7 (group I), between 7 and 13 (group II), and above 13 (group III).

age1= age<7;
age2=(7<=age)&(age<13);
age3=13<=age;



We compute the displacement and its length in each group. For instance, p12 is the mean displacement between group I and group II while p12L is its length.

p1= p(:,:,age1);
p2= p(:,:,age2);
p3= p(:,:,age3);

p12=mean(p2,3)-mean(p1,3);
p23=mean(p3,3)-mean(p2,3);

p1L=squeeze(sqrt(sum(p1.^2,2)));
p2L=squeeze(sqrt(sum(p2.^2,2)));
p3L=squeeze(sqrt(sum(p3.^2,2)));

p12L=sqrt(sum(p12.^2,2));
p23L=sqrt(sum(p23.^2,2));


To visualize the group difference between group I and II, we display the mean vector difference p12 on top of its length  p12L as quivers. Unfortunately, the quivers are small so we artificially enlarged them by the factor of five. The age effect is then visualized in Figure 1.

s=voxelsize(44);

figure; figure_surf(tt, s*p12L);
figure_quiver3surface(tt, s*5*p12, 5);
caxis([0 10*s]); colorbar
title('Between age group I and II')
colormap('hot'); view([90 60]); camlight;

figure; figure_surf(tt, s*p23L);
figure_quiver3surface(tt, s*5*p23, 5);
caxis([0 10*s]); colorbar
title('Between age group II and III')
colormap('hot'); view([90 60]); camlight;



Figure 1.
Mandibles are binned into three age groups: group I (ages 0 and 6), group II (ages
7 and 12) and group III (ages 13 and 19). Each row shows the mean group differences of
the displacement: group II - group I (first row) and group III - group II (second row). The
arrows are the mean displacement differences and colors indicate their lengths in mm.


Laplace-Beltrami Eigenfunction Based Heat Kernel Smoothing
January 27, 2013

The Laplace-Beltrami eigenfunction based heat kernel smoothing is first introduced in Seo et al. [2]. The schematic of heat kernel smoothing is given in Figure 2. We compute the eigenvalues D and eigenfunctions V of the Laplace-Beltrami operator along the template by the finite element method. The eigenvalues D are sorted in increasing order. To display the 5th eigenfunction V(:,5), we run

[A, C]= FEM(tt);
[V, D]= eigs(C,A, 1000,'sa'); 

figure; figure_surf(tt,V(:,5))
view([90 60]); camlight;



We smooth out the length of displacement with the bandwidth 20 and
200 basis functions.


p1Lsmooth=lb_smooth(p1L,tt, 20, 200, V, D);
p2Lsmooth=lb_smooth(p2L,tt, 20, 200, V, D);
p3Lsmooth=lb_smooth(p3L,tt, 20, 200, V, D);


Figure 2. Schematic of heat kernel smoothing. Given a functional data on a surface, we
compute the eigenfunctions ψj and the Fourier coefficients βj. Then we combine all the
terms and reconstruct the functional signal back.


Statistical Analysis on surface
January 27, 2013

The two-sample t-test is then performed using the built-in function ttest2.m. The code below computes the t-statistic map for comparing the group I and II. The other case is similarly done. The resulting t-static map is shown in not shown here but its square should yield the F-statistic map in Figure 3.

tfield=[];
for i=1:length(p1Lsmooth)
[h p ci t]=ttest2(p2Lsmooth(i,:),p1Lsmooth(i,:),[],[],'equal');
tfield=[tfield t.tstat];
end;

figure; figure_surf(tt, tfield);
figure_quiver3surface(tt, s*5*p12, 4);
view([90 60]); colormap('hot'); colorbar
title('T-stat between age group I and II')
caxis([-6 6]); colorbar; camlight


We can also compute it slightly using Hotelling's T-square routine hotelT2.m. Note that the square of t-statistic results should be the F-statistic (Figure 3).

tfield=[];
for i=1:length(p1Lsmooth)
    h= hotelT2(p2Lsmooth(i,:)', p1Lsmooth(i,:)');
    tfield=[tfield  h.t];
end;

figure; figure_surf(tt, tfield);
figure_quiver3surface(tt, s*5*p12, 4);
view([90 60]);

caxis([0 16]); colorbar
camlight; colormap('hot')

c=colormap;
colormap(c(end:-1:1,:));
title('Hotelling T-square between age group I and II')



Figure 3. F-statistic map showing the regions of significant mean displacement difference
shown in Figure 1.


Multiple Comparisons on Surface
January 27, 2013

The multiple comparisons corrected p-value is computed using the random field theory. It requires computing the Minkowski functional mu0, mu1 and EC-density rho0 and rho1. Let us compute the corrected p-value of F-statistic value of 9, 12 and 16 for the bandwidth 20.

t=[9 12 16]
a=1; b=sum(age1) + sum(age2) -2;
sigma=20;

mu0=2;

ss.faces = tt.faces;
ss.vertices= s*tt.vertices;
area=GETarea(ss);
mu2=sum(area.faces)/2;

rho0 = 1 - fcdf(t,a,b);
rho2 = 1/(4*sigma^2*pi)*gamma((a+b-2)/2)/ gamma(a/2)/...
       gamma(b/2)*(a*t/b).^((a-2)/2).*(1+a*t/b).^...
       (-(a+b-2)/2).*((b-1)*a*t/b-(a-1));

>>p=mu0*rho0+rho2*mu2

p =
    0.0192    0.0060    0.0014



This whole procedure is implemented in a single function.

F.tstat=[9 12 16]
F.df= [a b];

>>p=stat_Fcorrected(F, sigma,ss)

p =
         tstat: [9 12 16]
            df: [1 44]
    pcorrected: [0.0192 0.0060 0.0014]


Between F-statistic value of 8 and 20, we plotted the corrected p-value (Figure 4). Since the p-value plot is monotonic, It is not that hard to identify the threshold corresponds to the specific p-values. For instance, the F-statistic value corresponding to the corrected p-value of 0.01 is 8 for the case of comparing the groups I and II.

F.tstat=[8:0.01:20]
stat=stat_Fcorrected(F,sigma,ss)
figure; plot(F.tstat, stat.pcorrected, '-k', 'LIneWidth', 2)
title('Corrected p-value')
xlabel('F-statistic')
ylabel('p-value')

ind=find(stat.pcorrected<0.01);
>>F.tstat(ind(1))

ans =
     8


Figure 4. The plot of the corrected p-value over F-statistic for testing the group
difference. The corrected p-value is computed using the random field theory.


References
January 2, 2015

  1. Chung, M.K., Taylor, J. 2004. Diffusion Smoothing on Brain Surface via Finite Element Method,  IEEE International Symposium on Biomedical Imaging (ISBI). 562.
  2. 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.
  3. Seo, S., Chung, M.K. 2011. Laplace-Beltrami eigenfunction expansion of cortical manifolds. IEEE International Symposium on Biomedical Imaging (ISBI).
  4. Chung, M.K., Qiu, A., Seo, S., Vorperian, H.K. 2014. Unified heat kernel regression for diffusion, kernel smoothing and wavelets on manifolds and its application to mandible growth modeling in CT images arXiv:1409.6498. Medical Image Analysis 22:63-76.