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

Department of Biostatistics and Medical Informatics

Waisman Laboratory for Brain Imaging and Behavior

University of Wisconsin-Madison

Description

January 6, 2008

The weighed
spherical
harmonic (weighted-SPHARM) representation is a
unified surface data smoothing, surface parametrization and
surface registration technique formulated in a unified
Hilbert space framework. The weighted-SPHARM expresses
surface data as a weighted linear combination of spherical
harmonics. The weighted-SPHARM generalizes the traditional
SPHARM representation as a special case. If
you are using the Matlab codes/sample data below for
your publication, please reference [1], [2] or [4] in
the References. The codes have been tested under Matlab
versions 6.1, 6.5 and 7.5.

Autism cortical surface mesh data

April 15, 2013

The data set used in [1],[2][4] is stored Matlab MAT file brainimaging.waisman.wisc.edu/~chung/BIA/download/matlab.v1/autism.surface.mat. The data set consists of triangular meshes of brain surfaces. The MAT file contains the coordinates of the mesh vertices for 16 autistic (autisminner, autismouter) and 11 control (controlinner and controlouter) subjects. The variable tri is the collection of triangles that form each triangle element. To visualize the inner and outer surfaces of the first autistic subject (

surf.vertices=squeeze(autisminner(1,:,:))';

surf.faces=tri;

figure; figure_wire(surf, 'k', 'w');

surf.vertices=squeeze(autismouter(1,:,:))';

surf.faces=tri;

figure; figure_wire(surf, 'k', 'w');

To read the data set and perform weighted-SPHARM, run chapter11-SPHARM.m.The code is explained in detail below. The script requires additional function calls, which are found in brainimaging.waisman.wisc.edu/~chung/BIA/download/matlab.v1

Constructing
Weighted-SPHARM representation (MNI mesh format)

September 19,
2006

This implementation is based on the Montreal Neurological Institute (MNI) mesh format. Follow the step by step procedures below. It can construct up to the 85th degree representation and store them into a hard drive. It requires 2.4GB of space and 6 minutes in Pentium 4 computer. Download unitsphere.mat. It contains coordinates and topological information of the mesh of a unit sphere. The weighted-SPHARM representation is built upon this mesh. Download Y_l.m. It computes spherical harmonics Y_lm [1]. This is needed for SPHARMconstruct.m.

SPHARMconstruct(directory,85);

Download a sample outer cortical surface outersurface.obj and read it with mni_getmesh.m.

trisurf(tri,coord(1,:),coord(2,:),coord(3,:))

Run SPHARMsmooth.m to get the weighted-SPHARM representation. In order to smooth with up to 78th degree weighted-SPHARM with bandwidth 0.0001 (Figure 1), run

L=78

sigma=0.0001

[coord_smooth,fourier_coeff]= SPHARMsmooth(coord,L,sigma);

trisurf(tri,coord_smooth(1,:),coord_smooth(2,:),coord_smooth(3,:));

Figure 1.
Weighted-SPHARM representation of an outer cortical surface
with different bandwidths. The p-value plot automatically
determine the optimal degree for the representaion in a
forward model selection framework. See [2] for detail.

Constructing
Weighted-SPHARM representation (FreeSurfer mesh format)

June 20, 2008

Note: The FreeSufer format starts the face indexing from zero while the MATLAB can not handle zero indexing so we need to inrease the indexing by one. The codes below read the left pial surface (lh.pial) and the corresponding sphere mesh (lh.sphere) generated from FreeSurfer. If you don't have the corresponding sphere mesh (mainly obtained from a surface flattening algorithm), you can't construct the weighted-SPHARM for the given surface.

dir='lh.pial'

[vertices, faces] = read_surf(dir);

surf.vertices=vertices;

surf.faces=faces+1;

figure_wire(surf,'white');

dir='lh.sphere'

[vertices2, faces2] = read_surf(dir);

sphere.vertices=vertices2;

sphere.faces=faces2+1;

figure_wire(sphere,'white');

Figure 1-a. A left pial surface and the corresponding spherical mesh obtained from FreeSurfer.

Using the spherical mesh
lh.sphere, we can establish spherical angles \theta and
\varphi. The domains of \theta and \varphi follows the
convention established in [1, 2], which is different from the
MATLAB convension used in cart2sph.m. Our convension is more
widely used in mathematics literature. We check if we spherical
angles are properly computed:

[theta varphi]=EULERangles(sphere);

figure_trimesh(surf,theta,'rwb');

figure_trimesh(surf,varphi,'rwb');

figure_trimesh
projects the spherical angles \theta and \varphi to the cortical
surface. The spherical
angles are necessary to establish the weighted spherical harmonic
representation on a unit sphere. The weighted
spherical harmonic representation of degree 40 and bandwidth
sigma=0 is given by the code

figure_wire(surf_smooth,'white');

It
takes about 243, 624 and 1265 seconds for constructing degree 40,
60 and 80 representations respectively (Figure 1-b). For the choice of optimal degree, see
[2].

Figure
1-b. Weighted spherical harmonic representation of
various degree. Degree 40-5o would be sufficient for most
applications.

Extracting Spherical harmonic coefficients

January 6, 2008

The spherical harmonic coefficients
can be extracted for each weighted-SPHARM representation
possibly for data reduction and classification purposes. For
each representation, we have Fourier coefficients for x,y
and z coordinates. The above SPHARMsmooth.m generates a structured array called
fourier_coeff. If we are using the 42-th degree
representation for 12 subjects, we have

> fourier_coeff

fourier_coeff
=

x:
[43x85x12 double]

y:
[43x85x12 double]

z:
[43x85x12 double]

The first dimension is the degree between 0
and 42; the second dimesion is the order between -42 and 42;
the last dimension is the subject number. See Figure 2
for the display of the Fourier coefficients for subject 1.
For the l-th
degree, there are total 2l+1
degrees so not all elements in the matrix
contains Fourier coefficients. Figure 2 shows the upper
triangle elements in the matrix is padded with zeros (Note:
MATLAB convention for array indexing is that the downward
direction increases the array index). To make the Fourier coefficients up to degree
k into a single
feature vector, run the following code:

>k=42

>fourier_vector=SPHARMvectorize(fourier_coeff,k)

fourier_vector=

x:
[1849x12 double]

y:
[1849x12 double]

z:
[1849x12 double]

Note that there are total 1 + (2+1) + (2*2 +1
) + ... + (2*42+1) = (42+1)^2 = 1849 Fourier coefficients.

Figure 2.
Fourier coefficients for x-coordinates. The order
(horizontal) is arranged from -42 to 42. The upper triangle
elements in the matrix is padded with zeros.

Figure 3. Plot of
vectorized Fourier coefficients of x-coodinates for 16
autistic subject (red) and 12 control subjects (blue).

Reconstruction
from Spherical harmonic coefficients

January 9, 2008

Using the Fourier coefficients matrix
fourier_coeff, surface coordinates x,y,z can be
reconstructed.

>load unitsphere.mat

>surf=SPHARMrepresent2(sphere,
fourier_coeff,42);

To display the surface,

figure_trimesh(surf,surf_x,'rwb');

Figure 4.
Display of the 42-th degree weighted-SPHARM representation
of a cortical surface.

Surface
registration via SPHARM-correspondence

July 5, 2007

The weighted-SPHARM
can be used to establish surface
correspondence between subjects. The concept is
introduced in [2]. Download two sample outer surfaces: outsurface1.obj
and outsurface2.obj

>[tri,coord2,nbr,normal]=mni_getmesh('outersurface2.obj');

>coord(1,:,:)=coord1;

>coord(2,:,:)=coord2;

>[coord_smooth,coeff]=SPHARMsmooth(coord,78,0.0001)

>coeff1.x=coeff.x(:,:,1);

>coeff1.y=coeff.y(:,:,1);

>coeff1.z=coeff.z(:,:,1);

>coeff2.x=coeff.x(:,:,2);

>coeff2.y=coeff.y(:,:,2);

>coeff2.z=coeff.z(:,:,2);

Run SPHARMcorrespondence.m to obtain displacement vector field from the weighted-SPHARM of outsurface1.obj to the weighted-SPHARM of outsurface2.obj. The displacement is defined on resampled weighted-SPHARM surfaces.

>disp=SPHARMcorrespondence(coeff1,coeff2,sigma);

Note that the deformation vector field coord_smooth(1,:,:))+disp should exactly match coord_smooth(2,:,:)).

Figure 5. The trajectory
of deformation from the weighted-SPHARM
of outsurface1.obj
(1st surface) to the
weighted-SPHARM
of outsurface2.obj
(5th surface).

Gibbs Phenomenon
(Ringing Artifacts)

January 23,
2009

The Fourier
series are known to introduce the Gibbs phenomeon
(ringing artifacts) [1] [2]. The Gibbs phoenomeon occurs in
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.

Reduction of Partial Volume Effect

January 23, 2009

The discretization process in ditigal images are expected to produce the partial volme effect. By increasing the bandwidth \sigma of the weighted-SPHARM, we can reduce the partial volume effect. Figure 7 is generated by rounding off (x-35)^2/10^2 + (y-35)^2/10^2 + (z-35)^2/5^2 <1 in a 3D array of size 70 x 70 x 70. Obviously we are seeing the effect of the discretization.

Figure 7. The partial volume effect caused by the discretization process can be reduced by the proper choice of bandwidth \sigma.

References
on weighted-SPHARM

January 6, 2008

- Chung, M.K. Hartley, R., Dalton, K.M., Davidson, R.J. 2008. Encoding
cortical surface by spherical harmonics. Satistica Sinica
18:1269-1291

- Chung, M.K., Dalton, K.M., Davidson, R.J. 2008. Tensor-based
cortical
surface morphometry via weighed spherical harmonic
representation.
*IEEE Transactions on Medical Imaging*27:1143-1151. - Chung, M.K. Nacewicz, B.M., Wang, S., Dalton, K.M., Pollak, S., Davidson, R.J. 2008. Amygdala surface modeling with weighted spherical harmonics. 4th International Workshop on Medical Imaging and Augmented Reality (MIAR). Lecture Notes in Computer Science (LNCS) 5128:177-184
- 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. - Chung, M.K. Dalton, K.M., Davidson, R.J. 2007. Encoding
neuroanatomical
information using weighted spherical harmonic representation.
*IEEE Statistical Signal Processing (SSP) Workshop* - Shen, L., Chung, M.K. 2006. Large-scale
modeling
of parametric surfaces using spherical harmonics.
*Third International Symposium on 3D Data Processing, Visualization and Transmission (3DPVT)*. - Chung, M.K., Shen, L., Dalton, K.M.,
Davidson, D.J. 2006. Multi-scale
voxel-based
morphometry via weighed spherical harmonic
representation.
*International Workshop on Medical Imaging and Augmented Reality (MIAR)*. Lecture Notes in Computer Science (LNCS). 4091:36-43. - Chung, M.K., Robbins, S., Dalton,
K.M., Wang, S., Evans, A.C., Davidson, R.J. 2006. Tensor-based
cortical
morphometry via weighted spherical harmonic
representation.
*IEEE Computer Society Workshop on Mathematical Methods in Biomedical Image Analysis (MMBIA)*.

Last updated April 15, 2013