Weigthed Spherical Harmonic Representation

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

Department of Biostastics and Medical Infomatics
Waisman Laboratory for Brain Imaging and Behavior
University of Wisconsin-Madison

Description
January 6, 2008

The weighed spherical harmonic (SPHARM) representation is a unified surface data smoothing, surface parameterization and surface registration technique formulated in a unified Hilberst 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] or [2] in Referenes. The codes have been tested under Matlab versions 6.1, 6.5 and 7.5.


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.

directory='d:\harmonics\basis\'
 SPHARMconstruct(directory,85);
Download  a sample outer cortical surface outersurface.obj and read it with mni_getmesh.m. [tri,coord,nbr,normal]=mni_getmesh('outersurface.obj');
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

This implementation is based on the FreeSurfer mesh format. The FreeSurfer package produces cortical meshes with different mesh topology for different subjects. The implmentation can constructs up to the 85th degree representation. It has been tesed on a Mac computer (intel processor) with 2GB memory and MATLAB 7.5.

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

[surf_smooth, fourier]=SPHARMsmooth2(surf,sphere,40,0);
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 coeffients
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,coord1,nbr,normal]=mni_getmesh('outersurface1.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,:,:)).

SPHARPcorrespondence
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
January 6, 2008

  1. 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.
  2. 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. Special Issue of IEEE Transactions on Medical Imaging, on Computational Neuroanatomy. 26:566-581.
  3. Chung, M.K. Dalton, K.M., Davidson, R.J. 2007. Encoding neuroanatomical information using weighted spherical harmonic representation. IEEE Statistical Signal Processing (SSP) Workshop
  4. 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).
  5. 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).
  6. 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).