

\documentclass{llncs}
\usepackage{llncsdoc}
\usepackage{epsfig,epsf,psfrag}
\usepackage{amsmath,amsfonts}
\usepackage{amsmath}
\usepackage{amssymb}


%\documentclass[10pt,letterpaper]{article}
%\usepackage[dvips]{graphicx,graphics,color}
%\usepackage{times,amsmath,amsfonts}

\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}

\newtheorem{theorem}{\bf Theorem}
\newtheorem{algorithm}{\bf Algorithm}
\newtheorem{example}{\bf Example}

\usepackage[usenames,dvips]{color}
%\usepackage{color}

\title{Tiling Manifolds with Orthonormal Basis
\vspace{0.3cm}\\
\textcolor{BrickRed}{\large \tt University of Wisconsin, Madison\\ 
Department of Biostatistics and Medical Informatics\\
Technical Report 203}\\
\vspace{0.5cm}\\
{\large \tt \large Published in the Proceedings of the MICCAI 2008 Workshop\\ on 
Mathematical Foundations of Computational Anatomy}
%}
\author{Moo K. Chung\inst{1}\inst{2}, Anqi Qiu\inst{4}, Brendon, M. Nacewicz\inst{2},\\ 
%Kim M. Dalton\inst{2}, 
Seth Pollak\inst{3}, Richard J. Davidson\inst{2}\inst{3}}%}\\\
%\skiplinehalf
\institute{$^1$Department of Biostatistics and Medical Informatics, \\
$^2$Waisman Laboratory for Brain Imaging and Behavior,\\ 
$^3$ Department of Psychology and Psychiatry\\
University of Wisconsin, Madison, USA\\
$^4$Division of Bioengineering, Faculty of Engineering\\
National University of Singapore, Singapore\\
\email{mkchung@wisc.edu}}
\authorrunning{Chung et al.}}
%\authorinfo{E-mail: mkchung@wisc.edu, http://www.stat.wisc.edu/$\sim$mchung}

%GreenYellow, Yellow, Goldenrod, Dandelion, Apricot, Peach, Melon, YellowOrange, Orange, BurntOrange, Bittersweet, RedOrange, Mahogany, Maroon, BrickRed, Red, OrangeRed, RubineRed, WildStrawberry, Salmon, CarnationPink, Magenta, VioletRed, Rhodamine, Mulberry, RedViolet, Fuchsia, Lavender, Thistle, Orchid, DarkOrchid, Purple, Plum, Violet, RoyalPurple, BlueViolet, Periwinkle, CadetBlue, CornflowerBlue, MidnightBlue, NavyBlue, RoyalBlue, Blue, Cerulean, Cyan, ProcessBlue, SkyBlue, Turquoise, TealBlue, Aquamarine, BlueGreen, Emerald, JungleGreen, SeaGreen, Green, ForestGreen, PineGreen, LimeGreen, YellowGreen, SpringGreen, OliveGreen, RawSienna, Sepia, Brown, Tan, Gray
 

\begin{document}
  \maketitle 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
\begin{abstract}
One main obstacle in building a sophisticated parametric model along an arbitary anatomical manifold is the lack of an easily available orthonormal basis. Although there are at least two numerical techniques available for constructing an orhonormal basis such as the Laplacian eigenfunction approach and the Gram-Smidth orthogonalization, they are computationally not so trivial and costly. We present a relatively simpler method for constructing an orthonormal basis for an arbitrary anatomical manifold. On a unit sphere, a natural orthonormal basis is the spherical harmonics which can be easily computed.  Assuming the manifold is topologically equivalent to the sphere, we can establish a smooth mapping $\zeta $ from the manifold to the sphere. Such mapping can be obtained from various surface flattening techniques. If we project the spherical harmonics to the manifold, they are no longer orthonormal. However, we claim that there exists an orthonormal basis  that is the function of spherical harmonics and the spherical mapping $\zeta$. 

The detailed step by step procedures for the construction is given along with the numerical validation using amygdala surfaces as an illustration. As an application, we propose the pullback representation that reconstructs surfaces using the orthonormal basis obtained from an average template. The pullback representation introduces less inter-subject variability and thus requires far less number of coefficients than the traditional spherical harmonic representation. The source code used in the study is freely available at \\
\textcolor{BrickRed}{\tt http://www.stat.wisc.edu/$\sim$mchung/research/amygdala}.

\end{abstract}


%%%%%%%%% TITLE

%{\em Spherical Harmonics, Cortical Surface, Spectral Method, Orthonormal Basis}

\section{Introduction}

We present a novel orthonormal basis construction method for an arbitrary anatomical surface that is topologically equivalent to a sphere. The method avoids the well known Gram-Smidth orthogonalization procedure \cite{gorski.1994}, which is inefficient for high resolution polygonal meshes. In order to perform the Gram-Smidth orthogonalization as described in \cite{gorski.1994}, for a surface mesh with $n$ vertices, we need to perform the Choleski decomposition as well as the inversion of matrix of size $n \times n$. For a cortical mesh generated with FreeSurfer \cite{fischl.2000}, $n$ can easily reach up to 200000. 

On the other hand, Qiu et al. \cite{qiu.2006} constructed an orthonormal basis as the eigenfunctions of the Laplace-Beltrami operator in a bounded regions of interest (ROI) on a cortical surface (Figure \ref{fig:inner.sphere}). The finite element method (FEM) is used to numerically construct the orthonormal basis by solving a system of large linear equations. The weakness of the FEM approach is the computational burden of inverting a matrix of size $n \times n$. One advantage of the eigenfunction approach is that since the eigenfunctions and eigenvalues are directly related to the Laplace-Beltrami operator, it is trivial and geometrically intuitive to construct the heat kernel analytically and perform a various heat kernel smoothing based modeling \cite{chung.ni.2005}. 

We propose a completely different method that avoids the computational bottleneck by using a conceptually different machinery. We assume an arbitrary anatomical surface to be topologically equivalent to a sphere. Then using a smooth mapping $\zeta$ obtained from a surface flattening technique, we project the spherical harmonics to the anatomical surface. Obviously the projected spherical harmonics will no longer be orthonormal. However, if we  correct the metric distortion introduced from the surface flattening, we may able to  make the projected spherical harmonics orthonormal somehow. This is the basic idea behind our new proposed method. For the surface flattening, we present a new method that treats the mapping $\zeta$ as the geodesic path of the heat equilibrium state.

 
As an application of the proposed technique, we present the novel {\em pullback representation} for parameterizing anatomical boundaries that outperforms the traditional spherical harmonic (SPHARM) representation \cite{brechbuhler} \cite{chung.tmi.2007} \cite{gu.2004} \cite{shen.2004}  \cite{styner.2006}. We claim that our proposed representation has far less intersubject variability in the estimated parameters than SPHARM and converges faster to the true boundary with less number of basis. 



\section{Methods}

It is assumed that the anatomical boundary $\mathcal{M}$ is a smooth 2-dimensional Riemannian manifold parameterized by two parameters.  %\cite{chung.tmi.2007a}. 
The one-to-one mapping $\zeta$ from point ${\bf p}=(p_1,p_2,p_3)' \in \mathcal{M}$ to  ${\bf u}=(u_1,u_2,u_3)' \in S^2$, a unit sphere,  can be obtained from various surface flattening techniques such as conformal mapping  \cite{angenent.1999} \cite{gu.2004} \cite{hurdal.2004}, quasi-isometric mapping \cite{timsari.2000}, area preserving mapping \cite{brechbuhler} \cite{shen.2004} \cite{styner.2006} and the deformable surface algorithm \cite{macdonald.2000}. Since the conformal mapping tend to introduce huge area distortion, most spherical harmonic literature tend to use area preserving mapping \cite{brechbuhler} \cite{shen.2004} \cite{styner.2006}. 


In this paper, we present a new flattening technique via the geodesic trajectory of the equilibrium sate of heat diffusion. The proposed flattening technique is numerically simpler than any other available methods and does not require optimizing a cost function. The methodology is illustrated using the 47 amygdala binary segmentation obtained from the 3-Tesla magnetic resonance images (MRI). 


High resolution anatomical MRI were obtained using a 3-Tesla GE SIGNA scanner with a quadrature head coil. Details on image acquisition parameters are given in \cite{nacewicz.2006}. MRIs are reoriented to the pathological plane  for optimal segmentation and comparison with an atlas. This global alignment guarantee that amygdala are approximately aligned in the same orientation.

Manual amygdala segmentation was done by a trained expert and the reliability of the manual segmentation was validated by two raters on 10 amygdale resulting in interclass correlation of 0.95 and the intersection over the union of 0.84 \cite{nacewicz.2006}. Afterwards a marching cubes algorithm was used to extract the boundary of the binary segmentation as a triangle mesh with approximately 2000-3000 vertices. The amygdala surface is then mapped onto a sphere using the new flattening algorithm. 

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{stream.eps}
\end{center}
 \caption{The diffusion equation with a heat source (amygdala) and a heat sink (enclosing sphere) corresponds. After sufficient amount of diffusion, the heat equilibrium state is reached. By tracing the geodesic path from the heat source to the heat sink using the geodesic contours, we obtain a smooth mapping $\zeta$.}\label{fig:stream}
\end{figure}


\begin{figure}[t]
\begin{center}
%\begin{array}{cc}
\includegraphics[width=0.9\linewidth]{amyg_flattening.eps} 
%&\\(a) & (b) \\
%\end{array}
\end{center}
 \caption{Amygala surface flattening is done by tracing the geodesic path of the heat equilibrium state. The numbers corresponds to the different the geodesic contours. For simple shapes like amygde, 5 to 10 contours are sufficient for tracing the geodesic path.}\label{fig:flattening}
\end{figure}



\subsection{Diffusion-Based Surface Flattening}
 Given an  amygdala binary segmentation $\mathcal{M}_a$, we put a larger sphere $\mathcal{M}_s$ that encloses the amygala (Figure \ref{fig:stream} left). The amygdala and sphere serve as Dirichlet boundary conditions for solving the Laplace equation. The amygdala is assigned the value 1 while the enclosing sphere is assigned the value -1, i.e. 
 \bqn f(\mathcal{M}_a,\sigma)=1, f(\mathcal{M}_s,\sigma)=-1\label{eq:boundary}\eqn for all $\sigma$. The amygdala and the sphere serve as a heat source and a heat sink respectively. Then we solve an isotropic diffusion 
\bqn \frac{\partial f}{\partial \sigma} = \Delta f \label{eq:diffusion} \eqn
within the empty space bounded by the amygdala and the sphere. $\Delta$ is the 3D Laplacian. After enough diffusion, the system reaches the heat equilibrium state where the additional diffusion does not make any difference in the heat distribution (Figure  \ref{fig:stream} middle). Once we obtained the equilibrium state, we trace the geodesic path from the heat source to the heat sink for every mesh vertices. The trajectory of the geodesic path provides a smooth mapping from the amygdala surface to the sphere. The geodesic path can be easily traced by constructing geodesic contours that correspond to the level set of the equilibrium state (Figure \ref{fig:stream} right}). Then the geodesic path is constructed by finding the shortest distance from one contour to the next and  iteratively connecting the path together. Figure \ref{fig:flattening} shows the process of flattening using five contours corresponding to the temperature 0.6, 0.2, -0.2, -0.6, -1.0. 

Although we did not apply our flattening technique to other anatomical objects, the proposed method can be applied to more complex object than the amygdala. At the equilibrium state, we no longer has change in heat change over time, i.e. $\frac{\partial f}{\partial \sigma} =0$, so we have the Laplace equation $$\Delta f = 0$$ with the same boundary condition. The Laplace equation has been previously used to trace the distance between outer and inner cortical surfaces and to compute cortical thickness \cite{jones.2000}  \cite{lerch.2005} \cite{yezzi.2003}. Since the solution to the Laplace equation with the boundary condition (\ref{eq:boundary}) is unique even for highly convoluted and folded structures, the geodesic path will be uniquely defined.  

 
\subsection{Orhonormal basis in two sphere $S^2$}


Suppose a unit sphere $S^2$ is represented as a high resolution triangle mesh consisting of the vertex set $\mathcal{V}(S^2)$. We have used an almost uniformly sampled mesh with 2562 vertices and 5120 faces. Let us parameterize coordinates ${\bf u} \in S^2$ with parameters $\theta,\varphi$:
$$(u_1,u_2,u_3)=(\sin \theta \cos \varphi, \sin \theta \sin \varphi,
\cos \theta),$$ where $(\theta, \varphi) \in \mathcal{N}
=[0,\pi]\otimes [0,2\pi)$. The polar angle $\theta$ is the angle from the north pole and $\varphi$ is the azimuthal angle. The orthonormal basis on the unit sphere is given by the
eigenfunctions of $$\Delta f + \lambda f =0,$$ where $\Delta$
is the spherical Laplacian. The eigenfunction $Y_{lm}$ corresponding to the eigenvalue $l(l+1)$ is called the spherical harmonic of degree $l$ and order $m$ \cite{courant.vol2}. 
With respect to the inner product 
\bqn \langle f, g \rangle_{S^2} = \int_{S^2} f({\bf u})g({\bf u}) \;d\mu({\bf u}), \label{eq:inner}\eqn with measure $d\mu({\bf u}) = \sin \theta d\theta d\varphi$, $Y_{lm}$ form the orthonormal basis in $L^2(S^2)$, the space of square integrable functions on $S^2$, i.e. \bqn \langle Y_{lm}, Y_{l'm'} \rangle_{S^2} = \delta_{ll'}\delta_{mm'}. \label{eq:inner-sphere}\eqn %The norm is defined as $\| f \| = \langle f, f \rangle^{1/2}$. 

The inner product can be numerically computed as the Riemann sum over mesh vertices as
\bqn \langle Y_{lm}, Y_{l'm'} \rangle_{S^2} \approx \sum_{{\bf u}_j \in \mathcal{V}(S^2)} Y_{lm}({\bf u}_j)Y_{l'm'}({\bf u}_j) D_{S^2}({\bf u}_j),\label{eq:inner-sphere} \eqn
where $D_{S^2}({\bf u}_j)$ is the discrete approximation of  $d\mu({\bf u})$. Let $T_{{\bf u}_j}^1,T_{{\bf u}_j}^2, \cdots, T_{{\bf u}_j}^{j_m}$ be the area of triangles containing the vertex ${\bf u}_j$. Then we estimate $D_{S^2}({\bf u}_j)$ as 
\bqn  D_{S^2}({\bf u}_j) = \frac{1}{3} \sum_{k=1}^{j_m} T_{{\bf u}_j}^{k}. \label{eq:area-element}\eqn
%The summation is the total area of all triangles that contains vertex ${\bf u}_j$. 
The discrete approximation (\ref{eq:area-element})  defines the area of triangles at a mesh vertex.  The factor $1/3$ is chosen in such a way that 
$$\sum_{{\bf u}_j \in \mathcal{V}(S^2)} D_{S^2}({\bf u}_j) = 12.5514 = 4\cdot 3.1378,$$ analogous to the relationship
$$\int_{S^2} \; d\mu({\bf p}) = 4\pi.$$
The discrepancy between the integral and its discrete counter part is due to the mesh resolution and it should become smaller as the mesh resolution increases.


Based on the proposed discretization scheme, we have computed the inner product (\ref{eq:inner-sphere}) for all degrees $0 \leq l,l' \leq 20$. Figure \ref{fig:inner.sphere} (left) shows the inner products for every possible pairs. Since for up to the $k$-th degree, there are total $(k+1)^2$ basis functions, we have total $441^2$ possible inner product pairs, which is displayed as a matrix. For the diagonal terms, we obtained $0.9988 \pm 0.0017$ while for the off-diagonal terms, we have obtained $0.0000 \pm 0.0005$ indicating our basis and the discretization scheme is orthonormal with two decimal accuracy. 

    
\subsection{Orthonormal basis on manifold $\mathcal{M}$}


%Previously, in the spherical harmonic representation \cite{gerig.2001, gu.2004, shen.2004}, the surface coordinates ${\bf p} \in \mathcal{M}$ is parameterize as \bqn {\bf p} (\theta,\varphi) =\zeta^{-1} \circ {\bf u} = \sum_{l=0}^k \sum_{m=-l}^l {\bf p}_{lm}Y_{lm}(\theta,\varphi) \label{eq:coordinate.m},\eqn
%where ${\bf p}_{lm} = \langle {\bf p}, Y_{lm} \rangle$ is the vector of Fourier coefficient.

For $f, g \in L^2(\mathcal{M})$, the orthonormality is defined with respect to the inner product 
$$ \langle f, g \rangle_{\mathcal{M}} = \int_{\mathcal{M}} f({\bf p})g({\bf p}) \;d\mu({\bf p}).$$
Using the spherical harmonics in $S^2$, it is possible to construct an orthonormal basis in $\mathcal{M}$ numerically without the computational burden of solving the large matrix inversion associated with the eigenfunction method or the Gram-Smidth orthogonalization. 
Since the spherical harmonics are orthonormal in $S^2$ and, the manifolds $S^2$ and $\mathcal{M}$ can be deformed to each other by the mapping $\zeta$, one would guess that the orthonormal basis in $\mathcal{M}$ can be obtained somehow using the spherical harmonics.  Surprisingly this guess is not wrong as we will show in this section.

For $f \in L^2(S^2)$, let us define the {\em pullback operation} $^*$ as
$$\zeta^* f = f \circ \zeta.$$
While $f$ is defined on $S^2$, the pullbacked function $\zeta^*f$ is defined on $\mathcal{M}$. The schematic of the pull back operation is given in Figure \ref{fig:variance} (a). Then even though we do not have orthonormality on the pullbacked spherical harmonics, i.e.,
$$\langle \zeta^*Y_{lm}, \zeta^*Y_{l'm'} \rangle_{\mathcal{M}} \neq \delta_{ll'}\delta_{mm'},$$ 
we can make them orthonormal by using the Jacobian determinant of the mapping $\zeta$ somehow.

Consider the Jacobian $J_{\zeta}$ of the mapping $\zeta: {\bf p} \in \mathcal{M} \to {\bf u} \in S^2$ defined as
$$J_{\zeta} = \frac{\partial  {\bf u}(\theta,\varphi)}{\partial {\bf p}(\theta,\varphi)}.$$
For functions $f, g \in L^2(S^2)$,  we have the following change of variable relationship:
\bqn \langle f,g \rangle_{S^2} = \int_{\mathcal{M}}  \zeta^* f({\bf p}) \zeta^* g({\bf p}) |\det J_{\zeta}| \;d\mu({\bf p}). \label{eq:inner_S2}\eqn
Similarly we have the inverse relationship given as
\bqn \langle \zeta^*f,\zeta^*g \rangle_{\mathcal{M}} = \int_{S^2}  f({\bf u}) g({\bf u}) |\det J_{\zeta^{-1}}| \;d \mu({\bf u}). \label{eq:inner.M}\eqn
%$$\int_{S^2} f \;d\mu({\bf u}) = \int_{\mathcal{M}} |\det J_{\zeta}| \zeta^* f \;d\mu({\bf p}),$$
%$$\int_{\mathcal{M}} g \;d\mu({\bf p}) = \int_{S^2} |\det J_{\zeta^{-1}}| \zeta^{-1}^* g \; d\mu({\bf u}).$$
By letting $f=Y_{lm}$ and $g=Y_{l'm'}$ in (\ref{eq:inner_S2}), we obtain 
\bqn \delta_{ll'}\delta_{mm'} &=& \int_{\mathcal{M}}  \zeta^* Y_{lm}  \zeta^*Y_{l'm'}  |\det J_{\zeta}|  \;d\mu({\bf p})\label{eq:normality.pullback}\eqn
%\\
%&=& \langle  |\det J_{\zeta}|^{1/2}  \zeta^* Y_{lm}, |\det J_{\zeta}|^{1/2}  \zeta^* Y_{l'm'}  \rangle_{\mathcal{M}}.\eq
Equation (\ref{eq:normality.pullback}) demonstrates that functions 
\bqn Z_{lm} = |\det J_{\zeta}|^{1/2}  \zeta^*Y_{lm}\label{eq:Zlm1}\eqn are orthonormal in $\mathcal{M}$. We will refer $l$ as degree and $m$ as order of the basis function. Using the Riesz-Fischer theorem \cite{kolmogorov}, it is not hard to show that $Z_{lm}$ form a complete basis in $L^2(\mathcal{M})$. 
%In triangle meshes, the average area of adjacent triangles to a given vertex serves as the discrete estimation of the area element. 
%The discrete estimation of the area elements would break down if the triangle mesh resolution is crude. However we are using using high resolution mesh with 2562 vertices for fairly small amygdala and a unit sphere. 

\begin{figure*}[!tb]
\begin{center}
\begin{array}{cc}
\includegraphics[width=0.5\linewidth]{inner.anqi.eps} &
\includegraphics[width=0.5\linewidth]{anqi.eps} \end{array}
\end{center}
 \caption{Left: inner products of eigenfunctions of the Laplacian  for every pairs \cite{qiu.2006}. The pairs are rearranged from low to high degree. Right: representative eigenfunctions $\Psi_j$ on the left amygdala template surface obtained by solving $\Delta \Psi_j + \lambda_j \Psi_j =0$.}\label{fig:inner.sphere}
\end{figure}

\begin{figure*}[tb]
\begin{center}
\begin{array}{cc}
\includegraphics[width=0.5\linewidth]{inner.sphere.eps} &
\includegraphics[width=0.5\linewidth]{Zbais.eps} \end{array}
\end{center}
 \caption{Left: inner products of spherical harmonics computed using formula (\ref{eq:inner}) for every pairs. The pairs are rearranged from low to high degree and order. There are total $(20+1)^2=441$ possible pairs for up to degree 20. Right: representative orthonormal basis $Z_{lm}$ on the left amygdala template surface.}\label{fig:inner.sphere2}
\end{figure}

\subsection{Numerical Implementation} 

Although the expression (\ref{eq:Zlm1}) provides a nice analytical form for an orthonormal basis for an arbitrary manifold $\mathcal{M}$, it is not practical. If one want to use the basis (\ref{eq:Zlm1}), the Jacobian determinant needs to be numerically estimated somehow. We present a new discrete estimation technique for the surface Jacobian determinant that avoids estimating unstable spatial derivative estimation.  

The Jacobian determinant $J_{\zeta}$ of the mapping $\zeta$ can be expressed in terms of the Riemannian metric tensors associated with the manifolds $S^2$ and $\mathcal{M}$. Consider determinants $ \det g_{S^2}$ and $\det g_{\mathcal{M}}$ of the  Riemannian metric tensors associated with the parameterizations ${\bf u}(\theta,\varphi)$ and ${\bf p}(\theta,\varphi)$ respectively. Note that the integral of the area elements $ \sqrt{\det g_{S^2}}$ and $\sqrt{\det g_{\mathcal{M}}}$ with respect to the parameter space $\mathcal{N}$ gives the total area of the manifolds, i.e.
$$ \int_{\mathcal{N}} \sqrt{\det g_{S^2}} \;d\mu(\theta,\varphi) = 4\pi,\;
\int_{\mathcal{N}} \sqrt{\det g_{\mathcal{M}}} \;d\mu(\theta,\varphi) = \mu(\mathcal{M}).$$
Then we have the relationship 
$$ |\det J_{\zeta^{-1}}| =\frac{\sqrt{\det g_{\mathcal{M}}}}{\sqrt{\det g_{S^2}}},\; |\det J_{\zeta}| =\frac{\sqrt{\det g_{S^2}}}{\sqrt{\det g_{\mathcal{M}}}}.$$
Note that the Jacobian determinant $\det J_{\zeta}$ measures the amount of contraction or expansion in the mapping $\zeta$ from $\mathcal{M}$ to $S^2$. So it is intuitive to have this quantity to be expressed as the ratio of the area elements. Consequently the discrete estimation of the Jacobian determinant at mesh vertex ${\bf u}_j = \zeta({\bf p}_j)$ is obtained as
$$|\det J_{\zeta}| \approx \frac{D_{S^2}({\bf u}_j)}{D_{\mathcal{M}}({\bf p}_j)}.$$
Then our orthonormal basis is given by
\bqn Z_{lm}({\bf p}_j) = \sqrt{ \frac{D_{S^2}(\zeta({\bf p}_j))}{D_{\mathcal{M}}({\bf p}_j)} } \zeta^* Y_{lm}({\bf p}_j).\eqn
The numerical accuracy can be determined by computing the inner product
\bq \langle Z_{lm}, Z_{l'm'}  \rangle_{\mathcal{M}}  &\approx&  \sum_{{\bf p}_j \in \mathcal{V}(\mathcal{M})} Z_{lm}({\bf p}_j) Z_{lm}({\bf p}_j) {D_{\mathcal{M}}({\bf p}_j).\\
&=& \sum_{{\bf p}_j \in \mathcal{V}(\mathcal{M})} \zeta^*Y_{lm} ({\bf p}_j) \zeta^*Y_{l'm'} ({\bf p}_j)D_{S^2}( \zeta({\bf p}_j))\\
&=& \sum_{{\bf u}_j \in \mathcal{V}(S^2)} Y_{lm} ({\bf u}_j) Y_{l'm'} ({\bf u}_j)D_{S^2}({\bf u}_j)\\
&=&  \langle Y_{lm}, Y_{l'm'}  \rangle_{S^2} \eq
Since this is tautology, the order of the numerical accuracy in $Z_{lm}$ is identical to that of spherical harmonics given in the previous section. There is no need for additional validation other than given in the previous section. Hence we conclude that our basis is in fact orthonormal within two decimal accuracy. Figure \ref{fig:inner.sphere2} shows the result of our numerical procedure applied to the average amygdala surface template. The template surface is constructed by averaging the surface using the spherical harmonic correspondence given in \cite{chung.tmi.2007}. 

We have also constructed an orthonormal basis on a cortical surface with more than 40000 mesh vertices (Figure \ref{fig:cortex-zbasis}).  The diagonal elements in the inner product matrix are $0.9999 \pm 0.0001$ indicating that our basis is orthonormal within three decimal accuracy. As the mesh resolution increases, we expect to have increased accuracy. The proposed orthonormal basis construction methods avoids inverting matrix of size larger than $40000 \times 40000$ associated with the eigenfunction approach and the Gram-Smidth orthogonalization process. 


\begin{figure*}[!tb]
\begin{center}
\includegraphics[width=1\linewidth]{cortex-Zbasis.eps}&
\end{center}
 \caption{Orthonormal basis $Z_{lm}$ on a cortical surface. The basis is projected on a sphere to show how the nonuniformity of the Jacobian determinant is effecting the spherical harmonics $Y_{lm}$. The color scale is thresholded at $\pm 0.003$ for better visualization.}\label{fig:cortex-zbasis}
\end{figure}

Although the pattern of tiling in the eigenfunction approach (Figure \ref{fig:inner.sphere}) and the pullback based method (Figure \ref{fig:inner.sphere2}) looks different, it can be shown that they are actually linearly dependent. 

%Suppose $E_1$ is the vectorized version of the eigenfunctions \cite{qiu.2006} $\Psi_j$ over all mesh vertices, i.e.  
%$$E_1=(\Psi_j({\bf p}_1), \cdots, \Psi_j({\bf p}_n))'.$$
%Similarly we let $E_2$ be the vectorized version of the basis $Z_{lm}$ over all mesh vertices, i.e.
%$$E_2=(Y_{lm}({\bf p}_1), \cdots, Y_{lm}({\bf p}_n))'.$$
%Then the orthonormality condition is expressed as $E_1' D_1 E_1 =I $ with the area element matrix $D_1$ as defined in \cite{qiu.2006}. Similarly we have $E_2'D_2E_2 = I$ with the diagonal matrix $D_2 = diag(D_{\mathcal{M}}({\bf p}_1), \cdots, D_{\mathcal{M}}({\bf p}_n))$.
%
%With respect to the area element matrix $D_1$, $E_1'D_1 E_1 =I$, the identity matrix. Similarly with respect to our own area element matrix $D_2 = $
%
%
%Then there exists a matrix $A$ such that $Q=AM$. Note that we have $Q'KQ=I$ for the matrix of area element. Then we have $I=M'A'KAM$. Hence $A$ matrix should satisfies $A'KA=I$. $LQ=cQ$ implies $LAM=cAM$. Applying $A'K$ on the both sides, we have $LM=cM$. Hence we get the same eigenvalues.  





\section{Application: Pullback Representation}


As an application of the proposed orthonormal basis construction, we present a new variance reducing  Fourier Series representation that outperforms the traditional spherical harmonic representation \cite{brechbuhler} \cite{chung.tmi.2007} \cite{gu.2004} \cite{shen.2004}  \cite{styner.2006}. We will call this method as the {\em pullback representation}.

The spherical harmonic  (SPHARM) representation models the surface coordinates with respect to a unit sphere as
\bqn {\bf p}(\theta,\varphi) = \sum_{l=0}^k \sum_{m=-l}^l {\bf p}_{lm}^0Y_{lm}(\theta,\varphi) \label{eq:SPH} \eqn
where ${\bf p}_{lm}^0 =\langle {\bf p},Y_{lm} \rangle_{S^2}$ are spherical harmonic coefficients, which can be viewed as random variables. The coefficients are estimated using the {\em iterative residual fitting algorithm} \cite{chung.tmi.2007} that breaks a larger least squares problem into smaller ones in an iterative fashion. The MATLAB code for performing the iterative residual fitting algorithm for arbitrary surface mesh is given in
\textcolor{BrickRed}{\tt \small http://www.stat.wisc.edu}\\
\textcolor{BrickRed}{\tt \small /$\sim$mchung/softwares/weighted-SPHARM/weighted-SPHARM.hmtl}. Note that all MRIs were reoriented to the pathological plane guaranteeing an approximate global alignment before the surface flattening to increase the robustness of the coefficient estimation. 

The shortcoming of the spherical harmonic representation is that the reconstruction is respect to a unit sphere that is not geometrically related to the original anatomical surface. On the other hand, the pullback representation will reconstruct the surface with respect to the average template surface reducing substantial amount of variability compared to SPHARM. 

In the pullback representation, we represent the surface coordinates with respect to the template surface $\mathcal{M}$ as
\bqn {\bf p}(\theta,\varphi) = \sum_{l=0}^k \sum_{m=-l}^l {\bf p}_{lm}^1Z_{lm}(\theta,\varphi) \label{eq:WFS} \eqn
with ${\bf p}^1_{lm} = \langle {\bf p},Z_{lm} \rangle_{\mathcal{M}}$. 
Then we claim that the pullback representation has smaller variance in the estimated coefficients so that
\bqn \mbox{Var}({\bf p}_{lm}^1) \leq \mbox{Var}({\bf p}_{lm}^0). \label{eq:ineqaulity}\eqn
The equality in (\ref{eq:ineqaulity}) is obtained when the template $\mathcal{M}$ becomes the unit sphere, in which case the spherical mapping $\zeta$ collapses to the identity, and the inner products coincide.  We have computed the sample standard deviation of Fourier coefficients for 47 subjects using the both representations. In average, the SPHARM contains $88\%$ more intersubject variability compared to the pullback representation (Figure \ref{fig:variance} right). This implies that SPHARM is an inefficient representation and requires more number of basis to represent surfaces compared to the pullback method. 

\begin{figure*}[!tb]
\begin{center}
\begin{array}{cc}
\includegraphics[width=0.4\linewidth]{function-diagrm.eps}&
\includegraphics[width=0.55\linewidth]{WFSvariance.eps}
%a)&(b)
\end{array}
\end{center}
 \caption{Left: schematic showing how the pullback operation $^*$ is working. Point ${\bf p} \in \mathcal{M}$ is mapped to ${\bf u} \in S^2$ via our new flattening technique. As an illustration $f = Y_{3,2} + 0.6Y_{2,1}$ is plotted on $S^2$. The function $f$ is pulled back onto $\mathcal{M}$ by $\zeta$. Right: sample standard deviation of Fourier coefficients of for 47 subjects plotted over the index of basis. In average, the traditional SPHARM representation (black) has $88\%$ more variability than the pull back method (red).}\label{fig:variance}
\end{figure}


Although the pullback method is more efficient, the both representations (\ref{eq:SPH}) and (\ref{eq:WFS}) converge to each other as $k$ goes to infinity. We have computed the squared Euclidean distance between two representations numerically (Figure \ref{fig:pullback}). In average, the difference is 0.0569 mm for 20 degree representation negligible for 1mm resolution MR. Figure  \ref{fig:pullback} also visually demonstrate that the pullback representation converges to the true manifold faster than SPHARM again showing the inefficiency of the SPHARM representation.

\begin{figure*}[!tb]
\begin{center}
%\begin{array}{cc}
%\includegraphics[width=0.5\linewidth]{WFSvariance.eps} &
\includegraphics[width=0.9\linewidth]{pullback.eps} 
%\end{array}
\end{center}
 \caption{Comparison of SPHARM and the pullback representations for degree 5 to 25. Red colored numbers are the average Euclidean distance between two representations in mm.}\label{fig:pullback}
\end{figure}



\section{Conclusion}

We have introduced a computationally efficient way to construct an approximate orthonormal basis on an arbitrary manifold by pulling back the spherical harmonics to the manifold and accounting for the metric distortion using the Jacobian determinant. The proposed technique is very general so that it can be applicable to other types of anatomical manifolds. The constructed basis on an amygdala is used to show the new pullback representation that reconstruct the manifold as linear combination of the basis functions. 

%\bibliography{reference.JUN.2008}
\bibliographystyle{plain}

\begin{thebibliography}{10}

\bibitem{angenent.1999}
S.~Angenent, S.~Hacker, A.~Tannenbaum, and R.~Kikinis.
\newblock On the laplace-beltrami operator and brain surface flattening.
\newblock {\em IEEE Transactions on Medical Imaging}, 18:700--711, 1999.

\bibitem{brechbuhler}
C.~Brechbuhler, G.~Gerig, and O.~Kubler.
\newblock Parametrization of closed surfaces for 3d shape description.
\newblock {\em Computer Vision and Image Understanding}, 61:154--170, 1995.

\bibitem{chung.tmi.2007}
M.K. Chung, L.~Shen Dalton, K.M., A.C. Evans, and R.J. Davidson.
\newblock Weighted fourier representation and its application to quantifying
  the amount of gray matter.
\newblock {\em IEEE Transactions on Medical Imaging}, 26:566--581, 2007.

\bibitem{chung.ni.2005}
M.K. Chung, S.~Robbins, Davidson R.J. Alexander~A.L. Dalton, K.M., and A.C.
  Evans.
\newblock Cortical thickness analysis in autism with heat kernel smoothing.
\newblock {\em NeuroImage}, 25:1256--1265, 2005.

\bibitem{courant.vol2}
R.~Courant and D.~Hilbert.
\newblock {\em Methods of Mathematical Physics: Volume II}.
\newblock Interscience, New York, english edition, 1953.

\bibitem{fischl.2000}
B.~Fischl and A.M. Dale.
\newblock Measuring the thickness of the human cerebral cortex from magnetic
  resonance images.
\newblock {\em PNAS}, 97:11050--11055, 2000.

\bibitem{gorski.1994}
K.M. Gorski.
\newblock On determining the spectrum of primordial inhomogeneity from the cobe
  dmr sky maps: I. method.
\newblock {\em Astrophysical Journal, 430}, 430:L85, 1994.

\bibitem{gu.2004}
X.~Gu, Y.L. Wang, T.F. Chan, T.M. Thompson, and S.T. Yau.
\newblock Genus zero surface conformal mapping and its application to brain
  surface mapping.
\newblock {\em IEEE Transactions on Medical Imaging}, 23:1--10, 2004.

\bibitem{hurdal.2004}
M.~K. Hurdal and K.~Stephenson.
\newblock Cortical cartography using the discrete conformal approach of circle
  packings.
\newblock {\em NeuroImage}, 23:S119�S128, 2004.

\bibitem{jones.2000}
S.E. Jones, B.R. Buchbinder, and I.~Aharon.
\newblock Three-dimensional mapping of cortical thickness using laplace's
  equation.
\newblock {\em Human Brain Mapping}, 11:12--32, 2000.

\bibitem{kolmogorov}
A.N. Kolmogorov and S.V. Fomin.
\newblock Introductory real analysis.
\newblock Dover Publications, Inc, New York, 1970.

\bibitem{lerch.2005}
J.~P. Lerch and A.C. Evans.
\newblock Cortical thickness analysis examined through power analysis and a
  population simulation.
\newblock {\em NeuroImage}, 24:163--173, 2005.

\bibitem{macdonald.2000}
J.D. MacDonald, N.~Kabani, D.~Avis, and A.C. Evans.
\newblock Automated 3-{D} extraction of inner and outer surfaces of cerebral
  cortex from mri.
\newblock {\em NeuroImage}, 12:340--356, 2000.

\bibitem{nacewicz.2006}
B.M. Nacewicz, K.M. Dalton, T.~Johnstone, M.T. Long, E.M. McAuliff, T.R. Oakes,
  A.L Alexander, and R.J. Davidson.
\newblock Amygdala volume and nonverbal social impairment in adolescent and
  adult males with autism.
\newblock {\em Arch. Gen. Psychiatry}, 63:1417--1428, 2006.

\bibitem{qiu.2006}
A.~Qiu, D.~Bitouk, and M.I. Miller.
\newblock Smooth functional and structural maps on the neocortex via
  orthonormal bases of the laplace-beltrami operator.
\newblock {\em IEEE Transactions on Medical Imaging}, 25:1296--1396, 2006.

\bibitem{shen.2004}
L.~Shen, J.~Ford, F.~Makedon, and A.~Saykin.
\newblock surface-based approach for classification of 3d neuroanatomical
  structures.
\newblock {\em Intelligent Data Analysis}, 8:519--542, 2004.

\bibitem{styner.2006}
M.~Styner, I.~Oguz, S.~Xu, C.~Brechbuhler, D.~Pantazis, J.~Levitt, M.~Shenton,
  and G.~Gerig.
\newblock Framework for the statistical shape analysis of brain structures
  using spharm-pdm.
\newblock In {\em Insight Journal, Special Edition on the Open Science Workshop
  at MICCAI}, 2006.

\bibitem{timsari.2000}
B.~Timsari and R.~Leahy.
\newblock An optimization method for creating semi-isometric flat maps of the
  cerebral cortex.
\newblock In {\em The Proceedings of SPIE, Medical Imaging}, 2000.

\bibitem{yezzi.2003}
A.~Yezzi and J.L. Prince.
\newblock An eulerian pde approach for computing tissue thickness.
\newblock {\em IEEE Transactions on Medical Imaging}, 22:1332--1339, 2003.

\end{thebibliography}

\end{document}



%\subsection{Mesh Resampling}
%Using the {\em spherical harmonic correspondence} [blind], the polygonal mesh of every 
%amygdala surfaces will be resampled in such a way that every surface will have identical mesh topology with equal number of vertices and edges. Although it may not clear the motivation for how the mesh resampling can be used to construct spherical harmonic induced orthonormal basis. 


The main challenge in surface-based brain image analyses is the surface-specific normalization process that is needed to nonlinearly register a population of cortical surfaces to a template for a statistical analysis. We present a new  spectral method for normalizing cortical surfaces. By representing the Cartesian coordinates of cortical surfaces as a weighted linear combination of spherical harmonics (SPHARM), we obtain the spectral representation of cortical surfaces. Then the representations are matched in the least squares fashion to obtain the optimal surface displacement vector field. For validation purpose, we developed a new method for simulating cortical surfaces. The simulated cortical surface serves as the ground truth for validating our methodology.
As an application, we apply our technique in quantifying abnormal cortical anatomy in a group of high functioning autistic subjects. 
\end{abstract}


\section{Introduction}

The cortical surface normalization to a template is necessary for quantifying shape variations in a clinical population \cite{chung.ni.2005, thompson.1996}. Cortical surface normalization was mainly done by minimizing
an objective function that measures the global fit of two surfaces
while maximizing the smoothness of the deformation in such a way
that the gyral patterns are matched smoothly \cite{chung.ni.2005, thompson.1996}. These warping algorithms parameterize cortical surfaces and the warping function itself on the unit sphere for computational convenience. The resulting optimal surface displacement vector field is given implicitly without any parametric from. We present a new spectral method for 


\section{Cortical Thickness and Surface Warping}

Thickness is measured using two cortical surfaces: the white-matter/gray-matter
boundary and the gray-matter/CSF boundary. A single surface is used
to capture both hemispheres of the brain. Each surface is estimated
from the MR data using a deformable triangulated mesh. At each node
on the inner surface, a distance to the outer surface is computed
and taken to be a measure of cortical thickness. Finally, a spatial
normalization step is performed to facilitate inter-individual thickness
comparison. The spatial normalization used is a non-linear 2D warping
of each inner brain surface to a template inner brain surface.


For each brain image, the inner and outer surfaces are obtained by
iteratively deforming a triangulated mesh until it conforms to the
desired boundary. Such a \emph{deformable surface} approach has the
advantage that the surface topology can be fixed to be spherical and
the deformation process can maintain a non-intersecting surface at
all times \cite{macdonald00:cortical-surfaces}, obviating the need
for manual editing. The mesh starts as an ellipsoid located outside
the brain and is shrunk to obtain the inner cortical surface. Finally,
the inner surface is expanded, with constraints, to obtain the outer
cortical surface. 


There is a natural correspondence between nodes on the inner and outer
surface meshes, since the outer surface is obtained by deforming the
inner surface. The distance between corresponding nodes, dubbed $T_{link}$
in \cite{macdonald00:cortical-surfaces}, is the measurement used
here as cortical thickness. The thickness values of each individual
are mapped on to the individual's inner cortical surface.


A correspondence is required between inner surfaces of different individuals
in order to compare thickness values across subjects. We achieve this
by a 2D spatial normalization, i.e. nonlinearly warping the surface
of each individual to a template surface. 




In the WFS representation,
the surface normalization is straightforward and does not require
any sort of optimizations explicitly. 

The MRIs were spatially normalized into the Montreal neurological institute (MNI)
stereotaxic space using a global affine transform \cite{collins.1994}, which gives the 3D alignment of brain. Then an additional surface alignment was done during the cortical surface extraction process via the anatomic segmentation using the proximities (ASP) algorithm \cite{macdonald.2000}. The algorithm generates 40,962 vertices and 81,920 triangles with the identical mesh topology for all subjects. The vertices indexed identically on two cortical meshes will have a very close anatomic homology  and this defines the surface alignment \cite{chung.ni.2003} \cite{macdonald.2000}.  This provides the same spherical parameterization at identically indexed vertices across different cortical surfaces. Our approach avoids  the surface alignment done by coinciding the first order ellipsoid meridian and equator in the SPHARM-correspondence approach \cite{gerig.2001} \cite{martin.2006}.  Surface meshes obtained from other segmentation techniques such as FreeSurfer \cite{fischl.2000} may require the SPHARM-correpsondence approach. After this initial surface alignment, the more refined WFS-correspondence is established. 




