January 8, 2008
Using the weighed
harmonic (SPHARM) representation   ,
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  or .
The complete and processed amygdala surface mesh data set used in Chung et al. (2010)  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 , 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.
Loading T1-weighted MRI (Analyze format) into
January 8, 2008
To load the left amygdala of subject 001 (Tracing_mask_left.zip) we run
>leftinfo = analyze75info(lefthdr(1,:));
above line reads the header information (*.hdr). To find out
the image dimension, we run
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 =
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
rightinfo = analyze75info(righthdr(i,:));
vol = analyze75read(rightinfo);
save volumetry.mat left_vol right_vol
%two sample t-test on amygala
[h,p] = ttest2(al,cl)
[h,p] = ttest2(ar,cr)
4. Amygdala Surface
January 20, 2008
flattening from an amygdala surface to a unit sphere is
needed for establishing a spherical harmonic parameterization on
the amygdala surface  . 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
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
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 , 30
iterations were used but for simple shape like amygdala, 5
iterations are probably sufficient.
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=LAPLACEcontour(stream,sphere, 1 - 2*alpha/5);
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,
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
, 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
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
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
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  and  for detail. Figure 5 shows various
Figure 5. Spherical harmonic representations of an amygdala surface.
series are known to introduce the Gibbs phenomeon
(ringing artifacts) . 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
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
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.
%put surface flattening + weighted spherical harmonic modeling here
We can load *.mat files iteratively and save them into a structrued array study1_left:
The Fourier coefficinets of the
10th-subject is given by
x: [43x85 double]
y: [43x85 double]
z: [43x85 double]
Average Amygdala Surface Template
January 30, 2008; March 29, 2008
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.
To construct average surface template with only control subjects, modify the above code to
sigma=0.01 introduce a bit of smoothness into the
Figure 7. The average left and right amygdala surfaces over 47 subjects (24 control and 23 autistic subjects).
T-square Test for Two Samples
January 30, 2008
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.
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
h is the
Hotelling's T-square statistic and p is the corresponding
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.
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 . 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
May 8, 2008
January 6, 2008; October 15, 2011