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