(c) 2013 Moo K. Chung, Houri K. Vorperian

Vocal Tract Development Laboratory

Waisman Center

University of Wisconsin-Madison

**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

- Chung, M.K., Taylor, J. 2004. Diffusion
Smoothing
on
Brain
Surface
via Finite Element Method, IEEE International
Symposium on Biomedical Imaging (ISBI). 562.

- 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.
- Seo, S., Chung, M.K. 2011. Laplace-Beltrami eigenfunction expansion of cortical manifolds. IEEE International Symposium on Biomedical Imaging (ISBI).
- 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.