Description
January 8, 2008
Using the weighed
spherical
harmonic (SPHARM) representation [1] [2] [3][6],
amydala surfaces will be parameterized and discriminated in
a two group comparision setting. The implementation is based
on MATLAB 7.5 and a Mac-Intel computer (MacBook
Pro). If you are
using any of the codes or imaging data set given here,
please reference [1] or [6].
The complete and processed amygdala
surface mesh data set used in Chung et al. (2010) [6] is stored
as chung.2010.NI.mat. It
contains the surface displacement vector fields on the template
surface, group identifier and age and the template surface. The subset of
data published in [6], where the gaze fixation duration is
avaialble is stored as a seperate file amygdala.volume.mat. It consists
of the manual amygdala
segmention of MRI of 10 control
and 12 autistic subjects that also have gaze fixation
duration.
1.
Loading T1-weighted MRI (Analyze format) into
MATLAB
January 8,
2008
lefthdr =
/Users/chung/Desktop/amygdala/study/subject001/Tracing_mask_left.hdr
/Users/chung/Desktop/amygdala/study/subject002/Tracing_mask_left.hdr
/Users/chung/Desktop/amygdala/study/subject005/Tracing_mask_left.hdr
.....
To load the left amygdala of subject 001 (Tracing_mask_left.zip) we run
>leftinfo = analyze75info(lefthdr(1,:));
The
above line reads the header information (*.hdr). To find out
the image dimension, we run
>>
leftinfo.Dimensions
ans =
191 236
171 1
From Figure 1, we see that the
image dimension is Saggital 191 x Coronal 236 x Axial 171. To read
*.img, we run
>vol = analyze75read(leftinfo);
Using the marching cubes algorithm, we segment the amydala surface as a triangle mesh:
>surf = isosurface(vol)
surf =
Figure 1.
The coordinates directly correspond to the voxel positions of
Tracing_mask_left.img. Right figure is generated using MRIcro.
2. Data Set Published in
Chung et al. (2010), NeuroImage.
October 15, 2011
3. Amygdala Volumetry
Febuary 3,
2008
The traditional amygdala
volumetry can also be done. Make sure that vol only
contains values 1 and 0 otherwise, the following code will not
work.
for i=1:size(id,1)
rightinfo = analyze75info(righthdr(i,:));
vol = analyze75read(rightinfo);
leftvol(i)=sum(sum(sum(vol)));
rightvol(i)=sum(sum(sum(vol)));
end;
save volumetry.mat left_vol right_vol
load volumetry
al= left_vol(find(group))
ar=right_vol(find(group))
cl= left_vol(find(~group))
cr=right_vol(find(~group))
[mean(al) std(al)]
[mean(ar) std(ar)]
[mean(cl) std(cl)]
[mean(cr) std(cr)]
plot(al,ar,'.r')
hold on
plot(cl,cr,'.b')
legend('Autism','Control')
%two sample t-test on amygala
volume
[h,p]
= ttest2(al,cl)
[h,p]
= ttest2(ar,cr)
4. Amygdala Surface
Flattening
January 20,
2008
The surface
flattening from an amygdala surface to a unit sphere is
needed for establishing a spherical harmonic parameterization on
the amygdala surface [1] [6]. We have developed a novel
flattening techniue using heat diffusion. We first enclose the
amygdala binary mask with a larger sphere (Figure 2 left). Take the
amygadala as a source (value +1) and the sphere as a sink (value
-1).
[amyg,sphere,amygsphere]=CREATEenclosedamyg(vol,surf);
figure;imagesc(squeeze(amygsphere(:,25,:)));colorbar
With the heat sink and source, we perform isotropic
heat diffusion for long time. After suffiient amount of
time, we reach a stady state, which is equivalent to solving the
Laplace
quation. The following code will run approximately 1
min. per surface in a qaud-core computer. The
result is given in Figure 2.
The last arguemnt indicates the number of iterations used to
obtain the steady state heat equation. In publication [6], 30
iterations were used but for simple shape like amygdala, 5
iterations are probably sufficient.
stream=LAPLACE3Dsmooth(amygsphere,amyg,-sphere,5);
figure;imagesc(squeeze(stream(:,25,:)));colorbar
Figure 2. Left:
The amygdala is asigned the value 1 and the sphere enclosing
the amydala is asigned the value -1. Right: After solving
isotropic heat diffusion for long time, we reach a statedy
state, which can be used to generate a mapping from the
amygdala surface to the sphere by taking the geodesic path
from value 1 to -1.
The stady state map (Figure 2) is used to generate a mapping from the amygdala surface to the sphere. This is done by computing the geodesic path from the source to the sink. We first compute the levet sets of the stady state corresponding to f(p) = c for -1 <= c <=1 and propagate the amygdala boundary (c=1.0 in Figure 3) to the next level set (c=0.6) by computing the geodesic path. The propagation from one level set to the next is done by LAPLACEcontour. The flattening process is expected to produce discretization error. For each level set, we regularize mesh by shrinking down larger than expected area of triangles using REGULARIZEarea. 20% of the larger than expected triangles are reduced in size by setting the parameter to 0.8 in REGULARIZEarea.
sphere=isosurface(amyg);
for
alpha=1:5
sphere=LAPLACEcontour(stream,sphere,
1
- 2*alpha/5);
sphere=REGULARIZEarea(sphere,
0.8);
figure_wire(sphere,'yellow')
end;
Figure 3. Amygdala surface flattening process
using the geodesic path of of the stady state map. The
parameters correspond to the level set f(p)=1.0, 0.6, ..., -0.6,
-1.0.
5. Weighted Spherical
Harmonic Representation (SPHARM)
January 22,
2008
Once we establish the mapping
from the amygdala surface onto a sphere, we can determine the spherical
angles \theta and \varphi. The domain of \theta
and \varphi in this study follows the convention established in
[1][2][3][6], which is different from the result directly
obtained from the built-in MATLAB function cart2sph.m.
Continuing directly from the above for-end loop, we run
[theta varphi]=EULERangles(surf);
surf=isosurface(amyg);
figure_trimesh(surf,theta)
figure_trimesh(surf,varphi)
figure_trimesh
projects Euler angles \theta and \varphi to the amygdala surface
(Figure 4). The
spherical angles are necessary to establish the weighted
spherical harmonic representation on a unit sphere. The
additional MATLAB implementation details on the weighted
spherical harmonic representation can be found in this link [click
here].
Figure 4. The spherical
angles projected onto an amygdala surface. North pole is taken
where \theta=0 (blue peak in the left figure).
The weighted spherical harmonic
representation of degree 60 and bandwidth sigma=0 is given by
running the code
[surf_smooth,
fourier]=SPHARMsmooth2(surf,sphere,60,0);
figure_wire(surf_smooth,'yellow')
where surf is the amygdala
surface and sphere is the spherical mesh
obtained by running LAPLACEcontour
and REGULARIZEarea. When the bandwidth is
zero, i.e. sigma=0, the weighted spherical harmonic
representation becomes the traditional spherical harmonic
representation. See refernce [2] and [3] for detail. Figure 5 shows various
degree representation.
Figure 5. Spherical harmonic representations of
an amygdala surface.
The Fourier
series are known to introduce the Gibbs phenomeon
(ringing artifacts) [1]. The Gibbs phoenomeon occurs near
rapidgly changing or discontinous measurments. Since the
spherical harmonics are continuous and smooth basis functions,
it is not possible to represent discontinuity of data using
spherical harmonics.
Figure 6. The Gibbs
phenomeon is particuarly visible in the SPHARM representation
with k=42 and sigma=0. The ringing artifact is reduced if we
introduce a bit of smoothing sigma=0.001.
6. Saving and Loading
Multiple Surfaces
January
30, 2008
The following line of codes will process a
surface corresponding to subject id(i,:) and save the results
into a seperate *.mat file.
for i=1:size(id,1)
%put surface flattening + weighted
spherical harmonic modeling here
id_write=strcat('study1.right.',id(i,:),'.mat');
save(id_write,'fourier');
end
We can load *.mat files iteratively and save them into a structrued array study1_left:
for i=1:size(id,1)
file_name=strcat('study1.left.',id1(i,:),'.mat');
load(file_name);
study_left(i,:)=fourier;
end;
The Fourier coefficinets of the
10th-subject is given by
>
study_left(10,:)
ans =
x: [43x85 double]
y: [43x85 double]
z: [43x85 double]
7. Constructing
Average Amygdala Surface Template
January
30, 2008; March 29, 2008
The average
amygdala surfaces are constructed by averaging the collection
of Fourier coefficients up to degree 15. unitsphere2562.mat
is a spherical mesh with 2562 vertices. This is the mesh
resolution we will resample the average amygdala surface.
load unitsphere2562.mat
left_average=SPHARMaverage(study_left,sphere,
15,0.01);
figure_patch(left_average,'yellow');
right_average=SPHARMaverage(study_right,sphere,
15,0.01);
figure_patch(right_average,'yellow');
To construct average surface template with only control subjects, modify the above code to
left_average=SPHARMaverage(study_left(~group),sphere,
15,0.01);
The parameter
sigma=0.01 introduce a bit of smoothness into the
representation [1].
Figure 7. The average left
and right amygdala surfaces over 47 subjects (24 control and 23
autistic subjects).
8. Hotelling's
T-square Test for Two Samples
January 30,
2008
From
study_left and study_right that contain Fourier coefficients, we
construct the weighted spherical harmonics with degree 15
and bandwidth sigma 0.01. The resulting coordinates are
stored as matrices left_surf and right_surf.
left_surf=zeros(2562,3,47);
right_surf=zeros(2562,3,47);
for
i=1:47
temp=SPHARMrepresent2(sphere,study_left(i),15,0.01);
left_surf(:,:,i)=temp.vertices;
temp=SPHARMrepresent2(sphere,study_right(i),15,0.01);
right_surf(:,:,i)=temp.vertices;
end;
left_c=left_surf(:,:,find(~group));
left_a=left_surf(:,:,find(group));
right_c=right_surf(:,:,find(~group));
right_a=right_surf(:,:,find(group));
left_c
contails all surface coordinates for control subjects while
left_a contains all surface coordinates for autistic subjects. Then the Hotelling's T-square statistic is
used to measure the discrepancy between autistic and control
amygdala surfaces.
[h
p]=hotelT2(left_c,left_a);
h is the
Hotelling's T-square statistic and p is the corresponding
p-value (uncorrected).
figure_hotelT2(left_c,left_a,0.1,sphere);
Figure 8. The p-value of testing the
significance of group difference projected on to the mean
control surface. The white arrows indicate where to move the the
mean control surface to match it to the mean autistic surface. The
Hotelling's T-square did not detect significant difference (at
the corrected p-value of 0.05) although there is a huge cluster
on the right amygdala.
9. Multivariate
General Linear Model (MGLM) Using SurfStat
April 25,
2008
The problem with the Hotelling's T-square approach is the lack
of control for other covariates such as the global brain size
and age. So we need to covriate these nuisance variables
into a statistical model. This can be done if we use the multivarite
general linear models (MGLM) framework, which is implmented in Keith
Worsley's SurfStat
package [6]. In order to peform the
analysis given here, you need to download the package,
install it and set up the proper path.
Suppose we load age, brain size and group variables
(0=control, 1=autism) as
10. Brain and Behavior
Correlation
May 8, 2008
References
January 6, 2008;
October 15, 2011