MIDUS Amygdala/Hippocampus Surface Data Set, SPHARM and the Random Field Theory Based Surface Data Analysis Pipeline

(c) 2013 Moo K. Chung 
mkchung@wisc.edu

Department of Biostatistics and Medical Informatics
Waisman Laboratory for Brain Imaging and Behavior
University of Wisconsin-Madison

Description
August 7, 2013

MIDUS amygdala/hippocampus surface mesh data processing pipeline is presented. Note that this is not a public domain data yet. If you wish to use the data set in your research for whatever forms and shapes, you need to contact Moo K. Chung at mkchung@wisc.edu. The codes have been tested under a Macbook pro with 16GB memory and MATLAB 7.5. If you are using the Matlab codes below for your publication or research, please reference one of [2-5], where the detailed description of the code is given. The MATLAB script below is given as midus-v1.m.

1) MIDUS surface mesh data
August 7, 2013


The surface mesh data is stored in midus.mat, which contains template surface mesh and age, total brain volume, gender and displacement vector fields for 69 subjects. For the detailed explanation of the data set and relevent neuroanatomical hypothesis of interest, please read [1-4]. To display the template surface (Figure 1), run the codes below. 

load midus.mat

amygl=template{1,1}
amygr=template{1,2}
hippol=template{2,1}
hippor=template{2,2}

figure;
hold on; figure_wire(amygl,'black','red');
hold on; figure_wire(amygr,'black','blue');
hold on; figure_wire(hippol,'black','yellow');
hold on; figure_wire(hippor,'black','green');
view([-120 15])

Figure 1. Hippocampus/amygdala surface mesh templates.


The individual surface is constructed by adding the displacement vector field from the template. For instance, the surfaces of the 10-th subject (Figure 2) is constructed and displayed as

n=10
amygl_ind.vertices = amygl.vertices + displacement{n,1,1};
amygl_ind.faces=amygl.faces;

amygr_ind.vertices = amygr.vertices + displacement{n,1,2};
amygr_ind.faces=amygr.faces;

hippol_ind.vertices = hippol.vertices + displacement{n,2,1};
hippol_ind.faces= hippol.faces;

hippor_ind.vertices = hippor.vertices + displacement{n,2,2};
hippor_ind.faces=hippor.faces;

figure;
hold on; figure_wire(amygl_ind,'black','red');
hold on; figure_wire(amygr_ind,'black','blue');
hold on; figure_wire(hippol_ind,'black','yellow');
hold on; figure_wire(hippor_ind,'black','green');
view([-120 15])


Figure 2. Hippocampus/amygdala surface mesh for the 10-th subject. All the subject share the identical mesh topology that correspond across subjects. 

2) MIDUS surface mesh data in ASCII text file format
August 7, 2013

In case you need to port the surface mesh data into other programs, here is how to export it to ASCII text files. Each variable were saved into seperate text files using save command:

temp=amygl.vertices
save('amygl.vertices','temp', '-ascii')

temp=amygl.faces
save('amygl.faces','temp', '-ascii')

amygl.vertices is the vertex coordinates of the left amygdala surface while amygl.faces is the vertex indices that form triangle faces. Displacement vector fields with respect to the template is stored in the following fashion:

for i=1:69
    temp=displacement{i,1,1};
    strname=strcat('amygl.disp', num2str(i));
    save(strname, 'temp', '-ascii')
end;

for i=1:69
    temp=displacement{i,1,2};
    strname=strcat('amygr.disp', num2str(i));
    save(strname, 'temp', '-ascii')
end;

This generates 69 displacement vector field text files such as amygl.disp1, amygl.disp2, ... The numbers correspond to the subject id. All variables are then stored as midus.text.zip. If you unzip it, it will create the directory called midus.text containing all the variables as separate text files.

3) What can go wrong with spherical harmonic (SPHARM) representation
August 10, 2013

In order to understand the limitation of SPHARM (spherical harmonic) representation, see [3]. SPHARM representation depends on the quality of the corresponding spherical meshes. If the mesh is not so regular, SPHARM may not correctly reconstruct the original surfaces. Here we use area-preserving surface flattening which iteratively regularize triangle areas according to A-matrix in the finite element method (FEM) [6,7]. For instance, left amygdala and hippocampus surfaces are flattened as follows (Figure 3, 4):

figure; subplot(1,2,1); figure_wire(amygl,'r','w');
sphal=mesh_areapreserve(amygl,10000);
subplot(1,2,2); figure_wire(sphal,'r','w');

figure; subplot(1,2,1); figure_wire(hippol,'y','w');
sphhl=mesh_areapreserve(hippol,10000);
subplot(1,2,2); figure_wire(sphhl,'y','w');

Figure 3. The left amygdala template flattened onto a unit sphere by area-preserving mapping. Due to simple spherical geometry of amygdala, it can be easily flattened.



Figure 4. The left hippocampus template flattened onto a unit sphere by area-preserving mapping. Due to the elongated hippocamus shape, it is not easy to flatten. At the south pole, we see irregularly high mesh densities. This will cause a serious problem in SPHARM. It is crucial to come up with a better surface flattening algorithm that works best for SPHARM.

Once we have the corresponding spherical meshes, we reconstruct the surface coordinates with respect to the spherical meshes. SPHARM formulation is based on the following publications that utilizes iterative residual fitting (IRF) algorithm [8,9]. For amygdala surfaces, due to the regularity of the corresponding spherical meshes, spherical SPHARM can reconstruct the original surfaces pretty well (Figure 5). However, there is a serious defect on reconstructing hippocampus surfaces (Figure 6).


figure; subplot(1,2,1); figure_wire(amygl,'r','w');
[SPHARMal, fourier]=SPHARMsmooth2(amygl,sphal,20,0);
subplot(1,2,2); figure_wire(SPHARMal,'r','w');

figure; subplot(1,2,1); figure_wire(hippol,'y','w');
[SPHARMhl, fourier]=SPHARMsmooth2(hippol,sphhl,20,0);
subplot(1,2,2); figure_wire(SPHARMhl,'y','w');


Figure 5. The left amygdala surface is reconstructed using SPHARM of degree 20. Due to the spherical geometry of the original amygdala surface, the corresponding spherical mesh is regular enough to reconstruct the amygdala surface.

Figure 6. The left hippocampus surface is reconstructed using SPHARM of degree 20. Due to the elongated shape of the hippocampus, it won't reconstruct the hippocampus shape properly.


4) Better spherical harmonic (SPHARM) representation with more uniform mesh resampling
August 15, 2013


The left hippocampus template consists of 4884 triangles and 2444 vertices. We have modified mesh_areapreserve.m in such a way that A-matrix is computed within each loop whenever a mesh is updated. Although it takes longer, it performs better as shown in Figure 6. However, it is still not good enough. 

figure; subplot(1,3,1); figure_wire(hippol,'y','w');
sphhl=mesh_areapreserve2(hippol,1000);
subplot(1,3,2); figure_wire(sphhl,'y','w');
[SPHARMhl, fourier]=SPHARMsmooth2(hippol,sphhl,20,0);
subplot(1,3,3); figure_wire(SPHARMhl,'y','w');

Figure 7. Left: left hippocampus template. Middle: corresponding spherical parametrization. Right: SPHARM with respect to the spherical parametrization. It performs better than Figure 6 although it still needs improvements.

To improve the reconstruction done in Figure 7, we resample the spherical mesh with a fine uniform spherical mesh. Based on triangle subdivision, sphere_tri.m generates uniform spherical mesh with 40962 vertices and 81920 triangles. Using sphere6, we subsequently resample nonregular spherical mesh sphhl. This results in more regularized spherical mesh S2hlres, where SPHARM will be constructed. See Figure 8 for SPHARM on the resampled spherical mesh that shows a better reconstruction.

sphere6= sphere_tri('ico',6,1);
figure; subplot(1,3,1); figure_wire(sphere6, 'y', 'w');

[S2hlres, ind] = sphere_resample(sphhl, sphere6);
hippol_res.faces=S2hlres.faces;
hippol_res.vertices=hippol.vertices(ind,:);
subplot(1,3,2);  figure_wire(hippol_res, 'y', 'w');

[hippolSP, fourier]=SPHARMsmooth2(hippol_res,S2hlres,40,0);
subplot(1,3,3); figure_wire(hippolSP,'y','w');


Figure 8. Left: spherical mesh sphere6 consisting of 81920 triangles. MNI and FreeSurfer also use the identical uniform spherical mesh. Middle: resampled left hippocampus surface with respect to sphere6. Right: degree 40 SPHARM of the resampled hippocampus. Compared SPHARM in Figure 6 and 7, uniform resampling drastically improves the reconstruction.  

The estimated spherical harmonic coefficients are stored in the variable fourier. Using the estimated coefficients fourier, we can resample hippocampus surface anyway we want. For instance, we can construct a smaller uniformly sampled spherical mesh sphere4 with 5120 triangles and resample the hippocampus surface with respect to sphere4 (Figure 9). Then the resampled left hippocampus will have 5120 triangles. Note that originally the left hippocampus had 4884 triangles. We have easily changed the mesh topology. However, the reconstructed hippocampus seems unnatural. This can be fixed by reconstructing SPHARM with the initial irregular spherical mesh sphhl that did not work previously. Figure 9 shows the reconstructed SPHARM that looks the best among all SPHARM reconstructions.

sphere4= sphere_tri('ico',4,1)

SPHARMsurf1=SPHARMrepresent2(sphere4,fourier,40,0)
figure; subplot(1,2,1); figure_wire(SPHARMsurf1,'y','w');

SPHARMsurf2=SPHARMrepresent2(sphhl,fourier,40,0)
subplot(1,2,2); figure_wire(SPHARMsurf2,'y','w');


Figure 9. Left: SPHARM reconstructed hippocampus with 5120 triangles. The mesh topology is different from the original surface. Right: SPHARM reconstructed hippocampus with 4884 triangles. The mesh topology is identical to the original surface.

If you want to compute the reconstruction accuracy of SPHARM in terms of the mean squared errors (MSE), you need to compute between the original surface (Figure 7 left) and the best SPHARM reconstruction (Figure 9 right). 

5) Reconstruction error plots using the mean squared errors (MSE)
September 14, 2013

The accuracy of reconstruction can be visualized using the mean squared errors (MSE). Here we computed MSE between the original left hippocampus and the reconstructed and resampled 40 degree SPHARM and displayed (Figure 10). If you want to change the color bar scale, use caxis to reduce the scale of color bar.

residual=hippol.vertices - SPHARMsurf2.vertices;
MSE=sqrt(sum(residual.^2, 2));

figure; figure_surf(hippol , MSE)
view([160 20])
colorbar; camlight; alpha(0.8)
title('Reconstruction Error')
caxis([0 2])

 
Figure 10. The mean square error (MSE) map between the original left hippocampus and the reconstructed and resampled 40 degree SPHARM.

6) Random field theory based statistical analysis on surface.
August 18, 2013


Here we present statistical analyses for testing the group and age effects on the left hippocampus surface. The statistical analyses are based on SurfStat package written by Keith J. Worsley. You must download it and install it with the correct Matlab path.

On the template surface, we need to test for gender and age effect of total 69 surfaces. This is equivalent to testing the corresponding 69 displacement vector fields. The displacement of the left hippocampus are stored in displacement{:,2,1}. SurfStat package inputs only a specific vector data format so the displacement is stored in disphl.Similarly, meshes need to be reformatted.

disphl=zeros(69,2444,3);
for i=1:69
    disphl(i,:,:)=displacement{i,2,1};
end

hippol.coord= hippol.vertices';
hippol.tri=hippol.faces';


Then testing for the gender effect on the displacement vector is as follows. It is based on the general multivariate linear model (GMLM) or multivariate general linear model (MGLM). See [10] for statistical details on MGLM.

Gender=term(gender);
Age=term(age);
Brain=term(brain);

slm1 = SurfStatLinMod( disphl, 1+Brain + Gender+Age, hippol);
slm2 = SurfStatT( slm1, gender);
figure; figure_surf( hippol ,slm2.t);
colorbar; view([-120 20]); camlight

T-statistic obtained in MGLM is visualized in Figure 11. To determine the regions of statistical significance, we need to threshold the t-statistic map at 0.05. T-stat threshold of 4.73 corresponds to the multiple comparisons corrected pvalue of 0.05. However, we did not detect any region that is above 4.73.

thresholdhl=randomfield_threshold(slm, 0.05)


Figure 11. T-stat showing the lack of gender effect on the left hippocampus.


To testing age effect, we set up a similar MGLM. The T-statistic map is shown in Figure 12, where most of the tail regions of the hippocampus show significant age related anatomical change.

slm3 = SurfStatT( slm1, age);
figure; figure_surf( hippol ,slm3.t);
colorbar; view([-120 20]); camlight


Figure 12. T-stat showing the strong age effect on the left hippocampus. Any region above 4.73 is considered as significant at 0.05 level (corrected).

Publication on MIDUS dataset
August 7, 2013

  1. Chung, M.K. et al. 2013. Unified kernel regression on manifolds. under review.
  2. Chung, M.K. 2013. Statistical and Computational Methods in Brain Image Analysis.  Chapman & Hall/CRC. This book has chapters on various surface shape modeling techniques and corresponding MATLAB codes used on MIDUS data.
  3. Hosseinbor, A.P., Chung, M.K., Schaefer, S.M., van Reekum, C., Peschke-Schmitz, L., Sutterer, M., Alexander, A.L., Davidson, R.J. 2013. 4D Hyperspherical harmonic (HyperSPHARM) representation of  multiple disconnected brain subcortical structures, 16th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI). Lecture Notes in Computer Science (LNCS) accepted.
  4. Kim, S.-G., Chung, M.K., Schaefer, S.M., van Reekum, C., Davidson, R.J. 2012. Sparse shape representation using the Laplace-Beltrami eigenfunctions and its application to modeling subcortical structures. IEEE Computer Society Workshop on Mathematical Methods in Biomedical Image Analysis (MMBIA). 25-32. MATLAB
  5. Kim, S.-G., Chung, M.K., Seo, S., Schaefer, S.M., van Reekum, C., Davidson, R.J. 2011. Heat kernel smoothing via Laplace-Beltrami eigenfunctions and its application to subcortical structure modeling. Pacific-Rim Symposium on Image and Video Technology (PSIVT). Lecture Notes in Computer Science (LNCS). 7087:36-47. MATLAB
  6. Chung, M.K. 2001. Statistical morphometry in neuroanatomy, PhD Thesis, McGill University (PhD defense slides)
  7. Chung, M.K., Taylor, J. 2004. Diffusion smoothing on brain surface via finite element method,  IEEE International Symposium on Biomedical Imaging (ISBI). 562
  8. Chung, M.K. Hartley, R., Dalton, K.M., Davidson, R.J. 2008. Encoding cortical surface by spherical harmonics.  Satistica Sinica 18:1269-1291. MATLAB
  9. Chung, M.K., Dalton, K.M., Shen, L., L., Evans, A.C., Davidson, R.J. 2007. Weighted Fourier series representation and its application to quantifying the amount of gray matter.. IEEE Transactions on Medical Imaging 26:566-581. MATLAB
  10. Chung, M.K., Worsley, K.J., Nacewicz, B.M., Dalton, K.M., Davidson, R.J. 2010. General multivariate linear modeling of surface shapes using SurfStat. NeuroImage. 53:491-505.