
\documentclass[journal]{IEEEtran}
%\documentclass{article}


\usepackage{graphicx}  % Written by David Carlisle and Sebastian Rahtz
%\usepackage{epsfig,epsf,psfrag}
%\usepackage[dvips]{graphicx,graphics,color}
\usepackage{psfrag}
%\ifx\pdfoutput\undefined
%\usepackage{graphicx}
%\else
%\usepackage[pdftex]{graphicx}
%\fi



%\usepackage{url}       % Written by Donald Arseneau
                        % Provides better support for handling and breaking
                        % URLs. url.sty is already installed on most LaTeX
                        % systems. The latest version can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/other/misc/
                        % Read the url.sty source comments for usage information.

\usepackage[]{amsmath,amssymb}
%\usepackage{cite}
%this produces bold math symbols
%ex. \boldsymbol{\pi}

\newcommand{\bm}[1]{\boldsymbol{#1}}
\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}


\newtheorem{theorem}{\bf Theorem}
\newtheorem{algorithm}{\bf Algorithm}
\begin{document}
%
% paper title
\title{WFS Representation and Its Application to Quantifying the Amount of Gray Matter}
%
%
% author names and IEEE memberships
% note positions of commas and nonbreaking spaces ( ~ ) LaTeX will not break
% a structure at a ~ so this keeps an author's name from being broken across
% two lines.
% use \thanks{} to gain access to the first footnote area
% a separate \thanks must be used for each paragraph as LaTeX2e's \thanks
% was not built to handle multiple paragraphs
\author{Moo K.~Chung$^*$,~
        ~Kim~M.~Dalton,~
        ~Li~Shen,~
        ~Alan~C.~Evans,~
        and~Richard~J.~Davidson
\thanks{Manuscript received July 20, 2006; revised September 30, 2006  \em{Asterisk indicates corresponding author.}}        
\thanks{M.K. Chung is with the Department of Statistics, Biostatistics and Medical Informatics, and the  Waisman Laboratory for Brain Imaging and Behavior, University of Wisconsin, Madison, WI 53706 (e-mail: mchung@stat.wisc.edu).}
\thanks{K. M. Dalton and R. J. Davidson are with the Waisman Laboratory for Brain Imaging and Behavior, University of Wisconsin, Madison.}
\thanks{L. Shen is with the Department of Computer and Information Science, University of Massachusetts, Darthmouth.}
\thanks{A.C. Evans is with the McConnell Brain Imaging Centre, Montreal Neurological Institute, McGill University, Montreal.}}
 
% The paper headers
\markboth{IEEE Transactions on Medical Imaging}{Shell \MakeLowercase{\textit{et al.}}:
}

% use only for invited papers
\specialpapernotice{(Invited Paper)}

% make the title area
\maketitle


\begin{abstract}
We present a novel weighted Fourier Series (WFS) representation for cortical surfaces. The WFS representation is a data smoothing technique that provides the explicit smooth functional estimation of unknown cortical boundary as a linear combination of basis functions. The basic properties of the representation are investigated in connection with a self-adjoint partial differential equation and the traditional spherical harmonic (SPHARM) representation. To reduce steep computational requirements, a new iterative residual fitting (IRF) algorithm is developed. Its computational and numerical implementation issues are discussed in detail. The computer codes are also available at http://www.stat.wisc.edu/$\sim$mchung/softwares/weighted-SPHARM/weighted-SPHARM.html.

As an illustration, the WFS is applied in quantifying the amount of gray matter  in a group of high functioning autistic subjects. Within the WFS framework, cortical thickness and gray matter density are computed and compared.
\end{abstract}


\begin{keywords}
Spherical harmonics, cortical thickness, gray matter density, SPHARM, iterative residual fitting, diffusion smoothing \end{keywords}

\IEEEpeerreviewmaketitle



\section{Introduction}

 In this paper, we present a new morphometric framework called the {\em weighed Fourier series} (WFS)
representation. The WFS is both a global hierarchical parameterization and an explicit data smoothing technique formulated as a solution to a self-adjoint partial differential equation (PDE). WFS generalizes the traditional spherical harmonic (SPHARM) representation \cite{gerig.2001} \cite{shen.2004} with additional exponential weights. The exponentially decaying weights make the representation converges faster and reduce the ringing artifacts (Gibbs phenomenon) \cite{gelb} significantly. Unlike SPHARM, WFS can be reformulated as  heat kernel smoothing \cite{chung.ni.2005} when the self-adjoint operator becomes the Laplace-Beltrami operator. Then as similar to heat kernel smoothing, the random field theory \cite{worsley.1994} \cite{worsley.HBM.1996} \cite{worsley.2005} can be used for statistical inference on localizing abnormal shape variations in a clinical population. Many basic theoretical properties of WFS  and its numerical implementation issues are presented in great detail. The WFS representation requires estimating 18,723 unknown Fourier coefficients on a high resolution spherical mesh sampled at more than 40,000 vertices. This requires a specialized linear solver with fairly steep computational resources. To address this issue, we have developed a new estimation technique called the {\em iterative residual fitting} (IRF) algorithm \cite{shen.2006}. The computational burden has been reduced substantially by decomposing the subspace spanned by spherical harmonics into smaller subspaces, and iteratively performing the least squares estimation on the smaller subspaces.  The correctness of the algorithm is proven and its accuracy is numerically evaluated. 

As an illustration of the new technique, we have applied it in localizing abnormal amounts of gray mater in a group of high functioning autistic subjects. The cerebral cortex has a highly convoluted geometry and it is
likely that the local difference in gray matter concentration can
characterize a clinical population. Within the WFS framework, gray matter density and cortical thickness are also computed.  Then statistical parametric maps (SPM) are constructed to localize the regions of abnormal gray matter. These two measurements are also compared to determine if increased cortical thickness corresponds to increased gray matter locally. 


The main contributions of the paper are the development of the underlying theory of the WFS representation and its numerical implementation using the iterative residual fitting (IRF) algorithm. We also made the computer code freely available to public. In the following subsections, we will review the literature that is directly related to our methodology and address what our specific contributions are in the context of the previous literature. 

\subsection{Spherical Harmonic Representation}

The SPHARM representation \cite{brechbuhler} has been applied to subcortical structures such as the hippocampus and the amygdala \cite{gerig.2001} \cite{gu.2004} \cite{kelemen} \cite{shen.2004}. In particular, Gerig et al. used the mean squared distance (MSD) of the SPHARM coefficients in quantifying ventricle surface shape in a twin study \cite{gerig.2001}. Shen et al. used the principal component analysis technique on the SPHARM coefficients of schizophrenic hippocampal surfaces in reducing the data dimension \cite{shen.2004}. Recently it has begun to be applied to more complex cortical surfaces \cite{gu.2004} \cite{shen.2006}. Gu et al. presented SPHARM as a surface compression
technique, where the main geometric feasures are encoded in the low
degree spherical harmonics, while the noises are in the high degree
spherical harmonics \cite{gu.2004}.

In SPHARM, the spherical harmonic are used in  constructing the Fourier series expansion of the mapping from cortical surfaces to a unit sphere. So SPHARM is more of an interpolation technique than a smoothing technique, and thus it will have the ringing artifacts \cite{gelb}. On the other hand, WFS is a kernel smoothing technique given as a solution to a self-adjoint PDE. The solution to the PDE is expanded in basis functions. In a similar spirit, Bulow used the spherical harmonics in isotropic heat diffusion via the Fourier transform on a unit sphere as a form of hierarchical surface representation \cite{bulow.2004}. WFS offers many advantages over the previous PDE-based smoothing techniques \cite{andrade.2001} \cite{chung.ni.2003}. The PDE-based smoothing methods tend to suffer numerical instability \cite{andrade.2001} \cite{cachia.2003} \cite{chung.ni.2003} while WFS has no such problem.  Since the traditional PDE-based smoothing gives an implicit numerical solution, setting up a statistical model is not straightforward. However, WFS provides an explicit series expansion so it is easy to apply a wide variety of statistical modeling techniques such as the GLM \cite{friston.2002} , principal component analysis (PCA) \cite{shen.2004} and functional-PCA \cite{muller.2005} \cite{ramsay.1997}. 


SPHARM will be shown to be the special case of WFS.  In the SPHARM representation, all measurements are assigned equal weights and the coefficients of the series expansion is estimated in the least squares fashion. In WFS, closer measurements are weighted more and the coefficients of the series expansion is estimated in the weighted least squares fashion. So WFS is more suitable than SPHARM when the realization of the cortical boundaries, as triangle meshes, are noisy and possibly discontinuous. In most SPHARM literature, the degree of the Fourier series expansion has been arbitrarily determined and the problem of the optimal degree has not been addressed. Our WFS formulation addresses the determination of the optimal degree in a unified statistical modeling framework. The WFS-based global parametrization is computationally expensive compared to the local quadratic polynomial fitting \cite{brechbuhler} \cite{chung.ni.2003} \cite{dale.1999} \cite{joshi.1995} \cite{quicken} while providing more accuracy and flexibility for hierarchical representation.


\subsection{Cortical Thickness}

The cerebral cortex has the topology of a 2-dimensional convoluted sheet. Most of the features that distinguish these cortical regions can only be measured relative to that local orientation of the cortical surface \cite{dale.1999}. The CSF and gray matter interface is referred to as the {\em outer surface} (pial surface) while the gray and white matter interface is referred to as the {\em inner surface} \cite{chung.ni.2003} \cite{macdonald.2000}. Then the
distance between the outer and inner surfaces is defined as the {\em
cortical thickness} and it has been widely used as
an anatomical index for quantifying the amount of gray matter in the
brain \cite{chung.ni.2005} \cite{chung.ni.2003} \cite{fischl.2000}.


Unlike 3D whole brain volume based gray matter density, 1D cortical thickness measures have the advantage of providing a direct quantification of cortical geometry. It is likely that different clinical populations will exhibit different cortical thickness. By analyzing cortical thickness, brain shape differences can be quantified locally \cite{chung.ni.2003} \cite{fischl.1999} \cite{macdonald.2000} \cite{miller.2000}. The cortical surfaces are usually
segmented as triangle meshes that are constructed from deformable
surface algorithms \cite{davatzikos.1995} \cite{fischl.2000}
\cite{macdonald.2000}. Then the cortical
thickness is mainly defined and estimated as the shortest distance
between vertices of the two triangle meshes \cite{fischl.2000}
\cite{macdonald.2000}.  The mesh construction and discrete thickness
computation procedures introduce substantial noise in the thickness
measure \cite{chung.ni.2005} (Figure \ref{fig:thick}). So it is
necessary to increase the signal-to-noise ratio (SNR) and smoothness
of data for the subsequent random field based statistical analysis. For smoothing
cortical data, diffusion equation based methods have been used
\cite{andrade.2001} \cite{cachia.2003} \cite{chung.ni.2005} \cite{chung.ni.2003}. The shortcoming of these approaches is the need for numerically solving the diffusion equation possibly via the
finite element technique. This is an additional computational step on top of the cortical thickness estimation. In this paper, we present a more direct approach that smoothes and parameterizes the coordinates of a mesh directly via WFS such that the resulting thickness measures are already smooth. In WFS, the cortical surfaces are estimated as a weighted linear combination of smooth basis functions so that most algebraic operations on WFS will also be smooth. 

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{SSE.eps}
\end{center}\caption{Plots of RMSE over degree for bandwidths $t=0.01$, $0.001$, $0.0001$,
$0$. The smoothed outer surfaces are at the 85-th degree showing the bandwidth
controls the amount of smoothing in representing the cortex. The bandwidth $t=0$
corresponds to the traditional SPHARM. As $t \to 0$, the WFS converges to SPHARM.
\label{fig:SSE}}
\end{figure}

\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{IRF2.eps}
\end{center}\caption{Automatic degree selection in the IRF algorithm. For each bandwidth
$t$, the optimal degree is automatically selected by checking if
adding an additional degree will be statistical significant. The
outer cortical surfaces are the results of the optimal selection
procedure. The optimal degrees are $k=18 (t=0.01), k=42 (t=0.001),
k=52 (t=0.0005), k=78 (t=0.0001)$. For our study, $t=0.0001$ and the corresponding degree $k=78$ is used
through the paper. \label{fig:IRF}}
\end{figure*}

\section{Cauchy problem as a smoothing process}
Consider $\mathcal{M} \in \mathbb{R}^d$ to be a compact
differentiable manifold. Let $L^2(\mathcal{M})$ be the space of square integrable functions
in $\mathcal{M}$ with inner product \bqn \langle g_1, g_2 \rangle =
\int_{\mathcal{M}} g_1(p)g_2(p) \;d\mu(p),\label{eq:inner}\eqn where
$\mu$ is the Lebegue measure such that $\mu(\mathcal{M})$ is the
total volume of $\mathcal{M}$. The norm $\| \cdot \|$ is
defined as $$\|g\| = \langle g,g \rangle^{1/2}.$$
The linear partial differential operator $\mathcal{L}$ is {\em
self-adjoint} if
$$\langle g_1, \mathcal{L}g_2 \rangle = \langle
\mathcal{L}g_1, g_2 \rangle$$ for all $g_1, g_2 \in
L^2(\mathcal{M})$. Then the eigenvalues $\lambda_j$ and
eigenfunctions $\psi_j$ of the operator $\mathcal{L}$ are obtained
by solving \bqn \mathcal{L} \psi_j = \lambda_j
\psi_j.\label{eq:eigen}\eqn Without the loss of generality, we can order
eigenvalues $$0 < \lambda_0 \leq \lambda_1 \leq \lambda_2 \leq
\cdots$$ and the eigenfunctions to be orthonormal with respect to
the inner product (\ref{eq:inner}). Consider a Cauchy problem of the following form: \bqn \partial_t g +
\mathcal{L}g =0 \label{eq:cauchy1}, g(p,t=0) =f(p).\eqn
The initial functional data $f(p)$ can be further
stochastically modeled as
\bqn f(p) = \nu(p) + \epsilon(p),\label{eq:stochastic1}\eqn where
$\epsilon$ is a stochastic noise modeled as a mean zero Gaussian
random field and $\nu$ is the unknown signal to be estiamted. The PDE
(\ref{eq:cauchy1}) diffuses noisy initial data $f$ over time and estimate
the unknown signal $\nu$ as a solution. The time $t$ controls the
amount of smoothing and will be termed as the {\em bandwidth}. The unique solution to
equation (\ref{eq:cauchy1}) is given by the
following theorem.\\



\begin{theorem}\label{theorem:cauchy} For the self-adjoint linear differential operator
$\mathcal{L}$, the unique solution of the Cauchy problem (\ref{eq:cauchy1})  is given by
\bqn g(p,t) = \sum_{j=0}^{\infty} e^{-\lambda_j t} \langle f, \psi_j
\rangle \psi_j(p).\label{eq:cauchy.solution}\eqn
\end{theorem}
\begin{proof} For each fixed $t$,
$g$ has expansion \bqn g(p,t) = \sum_{j=0}^{\infty} c_j(t) \psi_j(p)
\label{eq:expansion}.\eqn Substitute equation (\ref{eq:expansion})
into (\ref{eq:cauchy1}). Then we obtain \bqn \partial_t c_j(t) +
\lambda_j c_j(t)=0.\label{eq:ODE}\eqn The solution of equation
(\ref{eq:ODE}) is given by $c_j(t) = b_je^{-\lambda_j t}$. So we
have solution $$g(p,t) = \sum_{j=0}^{\infty} b_j e^{-\lambda_j t}
\psi_j(p).$$ At $t=0$, we have
$$g(p,0) = \sum_{j=0}^{\infty} b_j \psi_j(p) = f(p).$$
The coefficients $b_j$ must be the Fourier coefficients $\langle
f,\psi_j \rangle$.\end{proof}

The implication of Theorem \ref{theorem:cauchy} is obvious. The
solution decreases exponentially as time $t$ increases and smoothes
out high spatial frequency noise much faster than low
frequency noise. This is the basis of many
of PDE-based image smoothing methods. PDE involving self-adjoint linear
partial differential operators such as the Laplace-Beltrami operator
or iterated Laplacian have been widely used in medical image
analysis as a way to smooth either scalar or vector data along
anatomical boundaries \cite{andrade.2001} \cite{bulow.2004} \cite{cachia.2003} \cite{chung.ni.2003}. These methods directly solve the PDE using standard numerical techniques such as the finite
difference method or the finite element method. The problem
with directly solving PDEs is the numerical instability and the
complexity of setting up the numerical scheme. WFS 
differs from these previous methods in such a way that we only need
to estimate the Fourier coefficients in a hierarchical fashion to solve the PDE.



\subsection{Weighted Fourier Series}
We will investigate the properties of the finite expansion of
(\ref{eq:cauchy.solution}) denoted by
$$\mathcal{F}_t^{k}[f](p)=\sum_{j=0}^{k} e^{-\lambda_j t} \langle f, \psi_j \rangle
\psi_j(p).$$ This expansion will be called as the {\em weighted Fourier
Series} (WFS). By rearranging the inner product, the WFS can be rewritten as kernel smoothing:
 \bqn \mathcal{F}_t^k[f](p) &=& \sum_{j=0}^{k}  e^{-\lambda_j t} \psi_j(p)\int_{S^2}
f(q)\psi_j(q)\;d\mu(q)\\
&=& \int_{S^2} f(q) K_t^k(p,q)
\;d\mu(q)\label{eq:kernel.smoothing}\eqn with symmetric positive
definite kernel $K_t^{k}$ given by \bqn K_t^{k}(p,q) =
\sum_{j=0}^{k} e^{-\lambda_j t}
\psi_j(p)\psi_j(q).\label{eq:kernel.tk}\eqn The subscript $t$ is
introduced to show the dependence of the kernel on time $t$.  This shows that the
solution of the Cauchy problem (\ref{eq:cauchy1}) can be interpreted as kernel
smoothing.

When the differential operator $\mathcal{L}=\Delta$, the
Laplace-Beltrami operator, the Cauchy problem (\ref{eq:cauchy1})
becomes an isotropic diffusion equation. For this particular case,
$K_t^{\infty}$ is called the {\em heat kernel} with bandwidth $t$ \cite{chung.isbi.2006}  \cite{chung.ni.2005}. For an arbitrary cortical manifold, the basis functions \psi_j$ can be computed and the exact shape of heat kernel can be determined numerically. Although it can be done by setting up a huge finite element method \cite{qiu.2006},  this is not a trivial numerical computation. A simpler approach is to use the first order approximation of the heat kernel for small bandwidth and iteratively apply it up to the desired bandwidth \cite{chung.ni.2005}.


\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{heat.kernel.shape.eps}
\end{center}
\caption{Shape of heat kernel with different bandwidth $t=0.01$,
$0.05$, $0.1$, $0.5$. The horizontal axis is the $\theta$ value from
the north pole $(\theta =0)$ to the south pole ($\theta = \pi$). The weighting scheme 
used in WFS follows the shape of the heat kernel weighting closer measurements heavily. 
 \label{fig:heat.kernel.shape}}
\end{figure}

WFS can be reformulated as a kernel regression problem \cite{fan}.
 At each fixed point $p$, we estimate unknown signal $\nu$ (\ref{eq:stochastic1}) with smooth function $h \in
 L^2(\mathcal{M})$ by minimizing the integral of the weighted
squared distance between $f$ and $h$: \bqn \min_{h \in
L^2(\mathcal{M})}\int_{\mathcal{M}} K_t(p,q)\big|f(q) - h(p)\big|^2
\; d\mu(q).\label{eq:min2}\eqn
%Under the stochastic model (\ref{eq:stochastic1}), we need to minimize the expectation
%$$  \min_{h \in
%L^2(\mathcal{M})}\mathbb{E} \int_{\mathcal{M}} K_t(p,q)\big|f(q) - h(p)\big|^2
%\; d\mu(q).$$
The minimizer of (\ref{eq:min2}) is
given by the following theorem.\\
\begin{theorem}\label{theorem:L2}\bq & & \frac{\mathcal{F}_t^k[f](p)}{\int_{\mathcal{M}} K_t^k(p,q)\;d\mu(q)}\\
& =& \arg\min_{h \in L^2(\mathcal{M})} \int_{\mathcal{M}}
K_t^k(p,q)\big|f(q) - h(p)\big|^2 \; d\mu(q).\eq
\end{theorem}
\begin{proof} Since the integral is quadratic in $h$, the minimum exists and
obtained when \bq &&\frac{\partial}{\partial h} \int_{\mathcal{M}}
K_t^k(p,q)|f(q) - h(p)|^2 \;d\mu(q) \\&=& -2\int_{\mathcal{M}}
K_t^k(p,q)[f(q) - h(p)] \;d\mu(q)=0.\eq Solving the equation, we
obtain the result.
\end{proof}
Theorem \ref{theorem:L2} shows WFS is the solution of a weighted least squares minimization problem.

When $\mathcal{L}$ is the Laplace-Beltrami operator with $k=\infty$, the heat kernel $K_t^{\infty}$ is a probability distribution in $\mathcal{M}$, i.e.
$$\int_{\mathcal{M}} K_t^{\infty}(p,q) \;d\mu(p) =1.$$
For this special case, Theorem \ref{theorem:L2} simplifies to
$$\mathcal{F}_t^{\infty}[f](p)= \arg\min_{h \in L^2(\mathcal{M})} \int_{\mathcal{M}}
K_t^{\infty}(p,q)\big|f(q) - h(p)\big|^2 \; d\mu(q).$$

In minimizing the weighted least squares in Theorem \ref{theorem:L2}, it is possible to restrict the function space $L^2(\mathcal{M})$ to a finite subspace that is more useful in numerical implementation. Let
$$\mathcal{H}_l =\{ \sum_{j=0}^l \beta_j \psi_j(p): \beta_j \in
\mathbb{R}\}$$ be the subspace spanned by basis
$\psi_0,\cdots,\psi_l$. Then  we have the following theorem.

\begin{theorem} \label{theorem:poly} If $\int_{\mathcal{M}} K_t^k(p,q) \;d\mu(q)=1$ and
$l \leq k$, then\bq \mathcal{F}_t^l[f](p)=\arg\min_{h \in
\mathcal{H}_l}\int_{\mathcal{M}} K_t^k(p,q)|f(q) - h(p)|^2 \;
d\mu(q).\eq
\end{theorem}

\begin{proof} Let $h(p)= \sum_{j=0}^l \beta_j \psi_j(p) \in \mathcal{H}_l.$ The
integral is written as $$ I(\beta_0, \cdots,
\beta_{l})=\int_{\mathcal{M}} K_t^k(p,q)\Big|f(q) - \sum_{j=0}^l
\beta_{j}\psi_j(p)\Big|^2\;d\mu(q).\label{eq:functional}$$ Since the
functional $I$ is quadratic in $\beta_0,\cdots,\beta_l$, the minimum
exists and it is obtained when $\frac{\partial I}{\partial
\beta_{j'}}=0$ for all $j'$. By differentiating $I$ and rearranging
terms, we obtain \bq &&\sum_{j=0}^k e^{-\lambda_j
t}\psi_j(p)\psi_{j'}(p) \int_{\mathcal{M}}
f(q)\psi_j(q)\;d\mu(q)\\
&=& \sum_{j=0}^l \beta_j \psi_j(p) \psi_{j'}(p).\eq Now integrate
the equations respect to measure $\mu(p)$ and obtain
$$e^{-\lambda_{j'}t}\langle f,\psi_{j'}\rangle = \beta_{j'}.$$
\end{proof}
If $K_t^k$ is a probability distribution, this theorem holds. For any other symmetric positive definite kernel, it can be made to be a probability distribution by renormalizing it. So Theorem \ref{theorem:poly} can be applicable in wide variety of kernels.


\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{FWHM.eps}
\end{center}\caption{Plot of FWHM (vertical) over bandwidth $t$ (horizontal) for both heat kernel and Gaussian kernel. The FWHM has to be numerically estimated in the case of the heat kernel. The numerically computed FWHM is used in the random field theory based multiple comparison correction. \label{fig:FWHM}}
\end{figure}

\subsection{Isotropic Diffusion on Unit Sphere}
Let us apply the WFS theory to a unit sphere denoted by $S^2$.  Since algebraic surfaces provide basis functions in a close form, it is not necessary to construct numerical basis \cite{qiu.2006}. The WFS in $S^2$ is given by the solution of isotropic diffusion. The spherical parametrization of $S^2$ is given
by the polar angle $\theta$ and the azimuthal angel $\varphi$: \bqn
p=(\sin \theta \cos \varphi, \sin \theta \sin \varphi, \cos
\theta)\label{eq:laplace2}\eqn with $p=(\theta, \varphi) \in
[0,\pi]\otimes [0,2\pi)$. The spherical Laplacian
$\Delta$ corresponding to the parametrization
(\ref{eq:laplace2}) is given by
$$\Delta = \frac{1}{\sin \theta}\frac{\partial }{\partial
\theta}\Big(\sin \theta \frac{\partial}{\partial
 \theta}\Big) +\frac{1}{\sin^2 \theta} \frac{\partial^2 }{\partial^2
 \varphi}.$$
There are $2l+1$ eigenfunctions $Y_{lm} ( -l \leq m \leq l)$,
corresponding to the same eigenvalue $\lambda_l =l(l+1)$ satisfying
 $$\Delta Y_{lm} = \lambda_l Y_{lm}.$$
$Y_{lm}$ is called the {\em spherical harmonic} of degree $l$ and
order $m$ \cite{courant, wahba.1990}. It is given explicitly as
 \bq Y_{lm} = \left\{\begin{array}{cc}
c_{lm}P^{|m|}_l(\cos\theta)\sin
(|m|\varphi), &-l \leq m\leq -1, \\
 \frac{c_{lm}}{\sqrt{2}}P_l^0(\cos\theta),& m=0,\\
c_{lm} P^{|m|
}_l(\cos\theta)\cos (|m|\varphi),& 1 \leq m\leq l,\\
\end{array}\eq
where $c_{lm}=\sqrt{\frac{2l+1}{2\pi}\frac{(l-|m|)!}{(l+|m|)!}}$ and
$P^m_l$ is the associated Legendre polynomials of order $m$.

Unlike many previous imaging
literatures on spherical harmonics that used the complex-valued spherical harmonics \cite{bulow.2004} \cite{gerig.2001} \cite{gu.2004} \cite{shen.2004}, only real-valued spherical harmonics  are used throughout the paper for convenience in setting up a
real-valued stochastic model.

For $f,h \in L^2(S^2)$, we define the inner product as
$$\langle f, h \rangle = \int_{\varphi=0}^{2\pi} \int_{\theta=0}^{\pi} f(p)
h(p) \;d\mu(p),$$ where Lebesgue measure $d\mu(p)=\sin \theta
d\theta d\varphi$. Then with respect to the inner product, the
spherical harmonics satisfies the orthonormal condition
$$ \int_{S^2}
Y_{ij}(p)Y_{lm}(p)\;d\mu(p) = \delta_{il}\delta_{jm},$$ where
$\delta_{il}$ is the Kroneker's delta. The kernel $K_t^k$ is given by
\bqn K_t^k(p,q) = \sum_{l=0}^{k} \sum_{m=-l}^{l} e^{-l(l+1)t}Y_{lm}(p)Y_{lm}(q). \label{eq:heat.kernel.sphere}\eqn
The associated WFS is given by
$$\mathcal{F}_t^k[f](p) = \sum_{l=0}^k \sum_{m=-l}^l
e^{-l(l+1)t}f_{lm}Y_{lm}(p)$$
with Fourier coefficient $f_{lm} = \langle f, Y_{lm} \rangle$. We will call this form of WFS as the
{\em weighted-SPHARM}. The special case $\mathcal{F}_0^k[f]$ is the traditional SPHARM representation used in representing the
Cartesian coordinates of anatomical boundaries \cite{gerig.2001} \cite{gu.2004} \cite{shen.2004}. Consider
subspace
$$\mathcal{H}_k = \{ \sum_{l=0}^k \sum_{m=-l}^l\beta_{lm} Y_{lm}:
\beta_{lm} \in \mathbb{R}\} \subset L^2(S^2),$ which is spanned by up
to the $k$-th degree spherical harmonics. Then the SPHARM satisfy
the least squares minimization problem different from Theorem
\ref{theorem:L2} and Theorem \ref{theorem:poly}.

\begin{theorem}\label{theorem:LSE}
\bqn \mathcal{F}_0^{k}[f](p) = \arg\min_{h \in \mathcal{H}_k} \|f -
h\|^2. \label{eq:LSE}\eqn
\end{theorem}



\section{Numerical Implementation}

In constructing the WFS representation, all we need is estimating Fourier coefficients $f_{lm} = \langle f, Y_{lm}\rangle$. There are three major techniques for computing the Fourier
coefficients. The first method numerically integrates the Fourier
coefficients over a high resolution triangle mesh \cite{chung.isbi.2006}. Although
this approach is the simplest to implement numerically and more accurate, the computation is extremely slow, due to the brute force nature of the technique. The second method is based on the fast Fourier
transform (FFT)  \cite{bulow.2004} \cite{gu.2004}. The drawback of FFT is the need for a
predefined regular grid system so if the mesh topology is different
for each subject as in the case of FreeSurfer \cite{fischl.2000}, a time
consuming interpolation is needed. The third method is based on
solving a system of linear equations \cite{gerig.2001} \cite{shen.2006} \cite{shen.2004} that minimize the least squares problem in Theorem
\ref{theorem:LSE}. This is the most widely used numerical technique
in the SPHARM literature. However, the direct application of the
least squares estimation is not desirable when the size of the
linear equation is extremely large.

Let $$h = \sum_{l=0}^k \sum_{m=-l}^l \beta_{lm}Y_{lm}(p) \in
\mathcal{H}_k.$$ Given $n$ nodes $p_1,\cdots, p_n$ in $S^2$ mesh, the
discretization of (\ref{eq:LSE}) is given by \bqn \|f-h\|^2 \approx
\sum_{i=1}^n \Big[f(p_i)-\sum_{l=0}^{k} \sum_{m=-l }^{l} \beta_{lm}
Y_{lm}(p_i)\Big]^2.\label{eq:discrete}\eqn The minimum of
(\ref{eq:discrete}) is obtained when \bqn f(p_i) = \sum_{l=0}^{k}
\sum_{m=-l }^{l} \beta_{lm} Y_{lm}(p_i)\label{equation:normal}\eqn
all $1 \leq i \leq n$.
The equation (\ref{equation:normal}) is referred as the {\em normal
equation} in statistical literatures. The normal equation is usually
solved via a matrix inversion. Let $$\bm{ f} =
(f(p_1),\cdots,f(p_n))'$$ and $$\bm{\beta}_l
=(\beta_{l,-l},\cdots,\beta_{l,l})'.$$ Also let
$${\bf Y}_l=\left[\begin{array}{ccc}
  Y_{l,-l}(p_1) & \cdots & Y_{l,l}(p_1) \\
  \vdots & \ddots & \vdots \\
  Y_{l,-l}(p_n) & \cdots & Y_{l,l}(p_n) \\
\end{array}\right]$$ be a $n \times (2l+1)$ submatrix consisting of
the $l$-th degree spherical harmonics evaluated at each node $p_i$.
Then (\ref{equation:normal}) can be rewritten in the following
matrix form: \bqn \bm{ f} = {\bf
Y}\bm{\beta}.\label{equation.normal2}\eqn with the {\em design
matrix} ${\bf Y} = [{\bf Y}_0, {\bf Y}_1,\cdots, {\bf Y}_k]$ and
unknown parameter vector $\bm{\beta} = (\bm{\beta}_0',\cdots,
\bm{\beta}_k')'$. The linear system is solved via \bqn \bm{\beta}=
({\bf Y}'{\bf Y})^{-1}{\bf Y}'\bm{ f}.\label{eq:normal3}\eqn The
problem with this widely used formulation is that the size of the
matrix ${\bf Y}$ is $n \times (k+1)^2$, which becomes fairly
large and may not fit in typical computer memory. So it becomes unpractical to perform matrix operation
(\ref{eq:normal3}) directly. This is true for many cortical surface
extraction tools such as FreeSurfer \cite{fischl.2000}  that
produces no less than $n> 100,000$ nodes for each hemisphere. This
computational bottleneck can be overcome by breaking the least
squares problem in the subspace $\mathcal{H}_k$ into smaller
subspaces using the {\em iterative residual fitting} (IRF) algorithm
\cite{shen.2006}. The IRF algorithm was first introduced in
\cite{shen.2006} although the the correctness of the algorithm was not given. In this paper, we present Theorem \ref{theorem:IRF} that proves the correctness of the IRF for the first time. The IRF algorithm can be also used in estimating SPHARM coefficients by letting the bandwidth $t=0$ in the algorithm.



\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.8\linewidth]{gibbs.eps}
\end{center}
\caption{Ringing artifacts (Gibbs phenomenon) for SPHARM and WFS. The traditional SPHARM and WFS representations are performed on the discontinuous measurements, which are defined as 1 in region  $ \frac{1}{8} < \theta < \frac{1}{4}$ and 0 in other regions. The first row shows severe ringing artifacts in the traditional SPHARM representation. The second row shows  the negligible ringing effect in the WFS representation. Columns are $k=18, 42, 52, 78$ degree representations. \label{gibbs}}
\end{figure*}


\subsection{Iterative residual fitting (IRF)
algorithm}Decompose the subspace $\mathcal{H}_k$ into smaller
subspaces as the direct sum: $$\mathcal{H}_k = \mathcal{I}_0 \oplus
\mathcal{I}_1 \cdots \oplus \mathcal{I}_k,$$ where subspace
$$\mathcal{I}_l = \{ \sum_{m=-l}^l \beta_{lm} Y_{lm}(p): \beta_{lm}
\in \mathbb{R} \}$$ is spanned by the $l$-th degree spherical
harmonics only. Then the IRF algorithm estimates the Fourier coefficients
$\bm{\beta}_l$ in each subspace $\mathcal{I}_l$ iteratively from
increasing the degree from $0$ to $k$. Suppose we estimated the
coefficients $\bm{\beta}_0, \cdots, \bm{\beta}_{l-1}$ up to degree
$l-1$ somehow. Then the residual vector $\bm{r}_{l-1}$ based on this estimation is given by \bqn
\bm{r}_{l-1} =\bm{ f} - \sum_{j=0}^{l-1} e^{-j(j+1)t}{\bf
Y}_{j}\bm{\beta}_{j}.\label{eq:residual}\eqn The components of the
residual vector $\bm{r}_{l-1}$ are identical so we denote all of them as
$r_{l-1}$. At the next degree $l$, we estimate the coefficients
$\bm{\beta}_l$ by minimizing the difference between the residual
$r_{l-1}$ and $\sum_{m=-l}^l \beta_{lm}Y_{lm} \in \mathcal{I}_l$.
This is formally stated as the following theorem:

\begin{theorem}\label{theorem:IRF}
\bqn \sum_{m=-l}^{l} f_{lm} Y_{lm} = \arg \min_{h \in
\mathcal{I}_{l}} \|r_{l-1} -h\|^2. \label{eq:r}\eqn
\end{theorem}
\begin{proof} Let $h = \sum_{m=-l}^l \beta_{lm}Y_{lm}(p) \in \mathcal{I}_l.$ The
squared norm is denoted as
$$I(\beta_{l,-l},\cdots,\beta_{l,l})= \int_{S^2}
\Big[\bm{r}_{l-1}(p) - \sum_{m=-l}^l \beta_{lm}Y_{lm}(p)\Big]^2
\;d\mu(p).$$ Since $I$ is quadratic in $\beta_{lm}$'s, the minimum
is obtained when $\frac{\partial I}{\partial \beta_{lm'}} =0$ for
all $m'$. By differentiating $I$ with respect to $\beta_{lm'}$ and
letting it equal to zero, we obtain \bqn &&\int_{S^2}
Y_{lm'}(p)\Big[\bm{r}_{l-1}(p) - \sum_{m=-l}^l
\beta_{lm}Y_{lm}(p) \Big]\;d\mu(p)\nonumber\\
&=& \int_{S^2} Y_{lm'}(p)\bm{r}_{l-1}(p)
\;d\mu(p)-\beta_{lm'}.\label{eq:1stterm}\eqn From
(\ref{eq:residual}), we have
$$r_{l-1} = f - \sum_{j=0}^{l-1} \sum_{m=-j}^j
e^{-j(j+1)t}\beta_{jm}Y_{jm},$$ which is a linear combination of
spherical harmonics up to $(l-1)$-th degree so it is orthonormal to
$Y_{lm}$. Then the first term in (\ref{eq:1stterm}) simplifies to
$$\int_{S^2} Y_{lm'}(p)\bm{r}_{l-1}(p) \;d\mu(p) = \langle f, Y_{lm'}
\rangle$$
\end{proof}
Theorem \ref{theorem:IRF}  proves that the correctness of the IRF algorithm. Then the discretization and the optimization is based on the normal equation approach (\ref{eq:LSE}):
$$\bm{\beta}_l = ({\bf Y}_l{\bf Y}_l)^{-1}{\bf Y}_l'\bm{r}_{l-1}.$$
Summarizing the results, the IRF algorithm is given below.

\begin{algorithm}\label{algorithm:IRF} Iterative Residual Fitting (IRF)
\begin{enumerate}
\item[1.] Let $l=0$.
\item[2.] $\bm{\beta}_{0} \leftarrow ({\bf Y}_0 {\bf Y}_0)^{-1}{\bf
Y}_0' \bm{ f}$.
\item[3.] $l \leftarrow l+1$.
\item[4.]$\bm{r}_{l-1} \leftarrow \bm{ f} - \sum_{j=0}^{l-1}e^{-j(j+1)t}{\bf Y}_{j}\bm{\beta}_{j}.$

\item[5.] $\bm{\beta}_l \leftarrow ({\bf Y}_l {\bf
Y}_l)^{-1}{\bf Y}_l'\bm{ r}_{l-1}.$

\item[6.] If $l < k$, go to step 3.

\end{enumerate}
\end{algorithm}




\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
degree $l$ & order $m$ & bandwidth $t$ & mean error & $f_{lm}$\\
\hline
18 & 17 & 0 &           0.0077 &  0.9979\\
18 & 17 & 0.0001 &    0.0078 &   0.9979\\
18 &  17 & 0.0005 &    0.0083&   0.9981\\
18 &17 & 0.01 &      0.0575 &    0.9995\\
42 & 41 & 0 & 0.0064    & 0.9977\\
42 & 41 & 0.001  &   0.0126 &     0.9992\\
52 &  51 & 0  &          0.0066  & 0.9972\\
52 & 51 & 0.0005 &    0.0101 &    0.9988\\
78 &  77 & 0 &            0.0060 &   0.9973\\
78 & 77 & 0.0001 &   0.0068 &    0.9984\\
\hline
\end{tabular}
\end{center}
\caption{Accuracy of the WFS representation}
\label{table1}
\end{table}

\subsection{Automatic optimal degree selection in IRF}
The IRF algorithm hierarchically builds the WFS from lower to higher
degree. In most previous SPHARM literature \cite{bulow.2004} \cite{gerig.2002}  \cite{gerig.2001}   \cite{gu.2004} \cite{shen.2006} \cite{shen.2004}, the issue of the
optimal degree has not been addressed. In SPHARM, the degree is simply selected based on a pre-specified error bound that depends on the size of anatomical structure. We present a statistical framework for automatically determining the optimal degree that does not depend on the size of anatomical structure. 

Although increasing the degree of WFS increases
the goodness-of-fit, it also increases the number of coefficients to be estimated quadratically. So it is necessary to
find the optimal degree where the goodness-of-fit and the number of
parameters balance out.

The Fourier coefficients $f_{lm}$ can be modeled to follow independent normal distribution
$N(\mu_{lm},\sigma_{l}^2)$. It is natural to assume an equal
variance within the same degree. This assumption is equivalent to the following $(k-1)$-th degree model: \bqn f(p_i) =
\sum_{l=0}^{k-1} \sum_{m=-l}^l
e^{-\lambda(\lambda+1)t}\mu_{lm}Y_{lm}(p_i) +
\bm{\epsilon}(p_i),\label{eq:model1}\eqn where $\bm{\epsilon}$ is a zero
mean isotropic Gaussian random field. Then at each iteration, we test if adding the $k$-th degree terms  in the
$(k-1)$-th degree model is statistically significant by testing the null hypothesis
$$H_0: \mu_{km} = 0 \mbox{ for } |m| \leq k.$$
The test statistic is constructed from the sum of squared errors
(SSE). Let the $k$-th degree {\em sum of squared errors} (SSE) be
$$\mbox{SSE}_k = \sum_{i=1}^n \bm{r}_k^2(p_i).$$
The plot of the  {\em root mean squared errors}  (RMSE), $\sqrt{\mbox{SSE}_k/n}$ for varying degree $5 \leq k \leq 85$ is shown in Figure \ref{fig:SSE}. As the degree $k$ increases, the root mean squared errors keep decreasing until it flattens out. So it is reasonable to stop the iteration when the decrease in error is no longer significant. Under $H_0$, the test statistic is
$$F =\frac{(\mbox{SSE}_{k-1} -\mbox{SSE}_{k})/(2k+1)}{\mbox{SSE}_{k-1}/(n-(k+1)^2)}
\sim F_{2k+1,n-(k+1)^2},$$ the $F$-distribution with $ 2k+1$ and
$n-(k+1)^2 $ degrees of freedom. We compute the $F$ statistic at
each degree and stop the IRF procedure if the corresponding P-value
first becomes bigger than the pre-specified significance $\alpha$
(traditionally at 0.01 or 0.05 and $\alpha=0.01$ in this study). For bandwidth $ t=0.0001$, the
optimal degree is determined to be $k=78$ (Figure \ref{fig:IRF}).






\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{thick.smooth.eps}
\end{center}
\caption{Cortical thickness measures projected onto an average
surface for better visualization. As the bandwidth increases from
$t=0.0001$ to $t=0.01$, the amount of smoothing also increases. The
first image shows the cortical thickness obtained from the
traditional deformable surface algorithm \cite{chung.ni.2005}
\cite{chung.ni.2003} \cite{macdonald.2000}. \label{fig:thick}}
\end{figure*}


\subsection{Computational cost and numerical accuracy}

The computation is done with MATLAB 6.1 in a Pentium 1.4 Ghz Intel M-processor laptop with 512MB memory. The complete MATLAB codes are freely available at http://www.stat.wisc.edu/$\sim$mchung/softwares/weighted-SPHARM/weighted-SPHARM.html. The IRF algorithm requires computing and storing the SPHARM bases in a hard drive. Sampling the SPHARM bases up to 85 degree at approximately 40,000 points on a unit sphere requires 2.4 GB of space and 16 minutes of computation. The SPHARM bases at each degree are stored as a single file and loaded into the computer at each degree in the IRF algorithm. For up to $k=78$ degree, there are total $3(k+1)^2 = 18,723$ unknown Fourier coefficents corresponding to the three Cartesian coordinates of a cortical surface. The average computational times for estimating all the Fourier coefficients are 354, 345 and 313 seconds respectively for 1, 24 and 48 cortical surfaces. The average computational times decreases as the number of cortical surfaces increases since the SPHARM bases are recycled for all surfaces at a given degree.

The numerical implementation of WFS is validated against the known analytical value. From the orthonormality of SPHARM, for $l \leq k$, $$e^{l(l+1)t}\mathcal{F}_t^k[Y_{lm}](p) =  Y_{lm}(p).$$
WFS of the spherical harmonic $Y_{lm}$ multiplied by $e^{l(l+1)t}$ should theoretically results in the original spherical harmonics. This exact analytical relation gives a criteria for numerical accuracy. The difference between the analytical value and the WFS representation gives the measure of numerical accuracy. Table \ref{table1} shows the result for selective $l$ and $m$ with various bandwidth $t$. The fourth column shows the mean absolute error computed over all mesh vertices as 
$$\frac{1}{n}\sum_{i=1}^n \Big|Y_{lm}(p_i) - e^{l(l+1)t}\mathcal{F}_t^k[Y_{lm}](p_i)\Big|.$$
As expected, the mean absolute error decreases as the degree increases. For the 78th degree with $t=0.0001$, the error is smaller than 2 decimal places. We have also estimated the Fourier coefficient $f_{lm}=\langle Y_{lm}, Y_{lm} \rangle$ to check if it is 1 in the fifth column. The estimation is accurate up to 2 decimal places for all degrees. A similar validation technique was introduced in \cite{chung.isbi.2006} for determining the accuracy of a numerical integration.


We have also checked the ringing artifacts (Gibbs phenomenon) \cite{gelb} in our numerical implementation. In both the traditional SPHARM and WFS, continuous measurements will have rapidly decaying Fourier coefficients so they converge faster. For continuous measurements such as the Cartesian coordinates and cortical thickness, we do not expect significant ringing artifacts for fairly large degree of harmonics. Note that we are using sufficiently high degree up to 78th harmonics, which corresponds to 6,241 independent basis functions. 

On the other hand, discontinuous measurements will have slowly decaying Fourier coefficients and thus the representations converge slowly. The ringing artifacts would be expected for the SPHARM representation as demonstrated in Figure \ref{gibbs}. In the figure, discontinuous measurements are constructed as a step function of value 1 in the circular band $\frac{1}{8} < \theta < \frac{1}{4}$ and 0 outside of the band. The SPAHRM  representation of the step function resulted in significant ringing artifacts even for fairly high degrees up to $78$. In comparison, the WFS representation of the step function does not exhibit any serious ringing artifacts as seen in Figure \ref{gibbs}. The superior performance of WFS can be easily explained in terms of convergence. The WFS representation additionally weights Fourier coefficients with exponentially decaying weights, which contributes more rapid convergence even for discontinuous measurements. This robustness of WFS is also related to the fact that it is a PDE-based data smoothing technique while the traditional SPHARM is more of interpolation or reconstruction technique.









\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{scale2.eps}
\end{center}\caption{Multi-scale representation of surface registration toward the average template. The top is the inner surface while the bottom is the outer surface. $\alpha=0$ is the
surface of one particular subject while $\alpha=1$ is the average
surface of 24 subject.  The scale space $\alpha \in [0,1]$ is searched for the maximal discrepancy between two groups for increased statistical sensitivity.\label{fig:multi-scale}}
\end{figure*}


\subsection{Computing FWHM}
Since WFS is a kernel smoothing method, it is useful to know the full width at the half maximum (FWHM) of the underlying kernel. The computed FWHM is later used in the random field based
multiple comparison corrections \cite{chung.ni.2005} \cite{chung.ni.2003} \cite{worsley.2005} \cite{worsley.HBM.1996}. Computing the FWHM of the heat kernel used in WFS is not trivial since there is no known close form expression for FWHM as a function of bandwidth $t$. Therefore, FWHM is computed numerically.

 For $p, q, r \in S^2$, let us define the Cartesian inner product $\cdot$ as
$p \cdot q = \cos (\theta),$ where $\theta$ is an angle between $p$ and $q$. The heat kernel (\ref{eq:heat.kernel.sphere}) is symmetric along the
geodesic circle. If $p \cdot q = p \cdot r$, we have $K_{t}^{k}(p,q) = K_{t}^{k}(p,r)$. This property
can be used to simplify the expansion (\ref{eq:heat.kernel.sphere}) using the harmonic addition theorem  \cite{groemer} \cite{wahba.1990}.

\begin{theorem} (Harmonic addition theorem) \label{theorem:HAT} \bqn \sum_{m=-l}^l Y_{lm}(p)Y_{lm}(q) =
\frac{2l+1}{4\pi}P_l^0(p \cdot q).\label{eq:addition} \eqn
\end{theorem}
\begin{proof} Fix the azimuthal angel
$\varphi=0$ and $p$ be the north pole, i.e. $p=(0,0,1)$. Now by
varying $q= (\sin \theta, 0, \cos \theta)$ for $0 \leq \theta \leq \pi$,
we have $Y_{lm}(\theta,\varphi)=0$ if $m \neq 0$. Then we have

$$\sum_{m=-l}^l Y_{lm}(p)Y_{lm}(q) = Y_{l0}(p)Y_{l0}(q)=
\frac{2l+1}{4\pi}P_l^0(1)P_l^0(\cos \theta).$$
Note that $P_l^0(p) = 1$. This implies that the sum of product of the spherical harmonics is
a function of the inner product between $p$ and $q$ only. From symmetry, if we rotate $p$ back to the original position from the north pole, the same result should hold.
\end{proof}

Using the harmonic addition theorem, we simplify the heat kernel in the following theorem.

\begin{theorem} \label{theorem:heat.kernel}  For any $p, q \in S^2$, \bqn K_{t}^k(p,q) = \sum_{l=0}^{k}
\frac{2l+1}{4\pi}e^{-l(l+1)t}P_l^0(p \cdot q). \eqn
\end{theorem}
Theorem \ref{theorem:heat.kernel} is used to plot the shape of the
heat kernel by fixing $p$ to be the north pole and by varying
$\theta = \cos^{-1}(p \cdot q)$ (Figure
\ref{fig:heat.kernel.shape}). Similar result is also given in \cite{bulow.2004}. The maximum of the kernel is obtained at $\theta = 0$.  Then the FWHM is solved numerically for $\theta$
in
$$  \frac{1}{2}\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-l(l+1)t} =
\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-l(l+1)t}P_l^0(\cos
\theta).$$ The FWHM is then $2\theta$. Figure \ref{fig:FWHM} shows
the nonlinear relationship between bandwidth $t$ and the
corresponding FWHM. When $t=0.0001$, the corresponding FWHM is 0.0597. This is the FWHM we have used in the subsequent cortical thickness analysis.





\section{Application to Cortical morphometry}
In this section, we present applications of WFS in localizing the regions of abnormal amount of gray matter in autistic subjects. Surface normalization issue is discussed in detail along with a stochastic model for WFS.


\subsection{Data Set}

$n_1=$ 12 high functioning autistic (HFA) and $n_2=$12 normal control
(NC) male subjects were screened to be right-handed. The autistic
subjects were diagnosed via the Autism Diagnostic Interview -
Revised (ADI-R) by a trained and certified psychologist at
the Waisman center at the University of Wisconsin-Madison
\cite{dalton.2005}. Age distributions for HFA and NC were compatible
at $15.93 \pm 4.71$ and $17.08 \pm 2.78$ respectively.

High resolution anatomical magnetic resonance images (MRI) were
obtained using a 3-Tesla GE SIGNA scanner with a quadrature head RF
coil. A three-dimensional, spoiled gradient-echo (SPGR) pulse
sequence was used to generate $T_1$-weighted images. Image intensity
nonuniformity was corrected using the nonparametric nonuniform
intensity normalization method \cite{sled} and then the image was
spatially normalized into the Montreal neurological institute (MNI)
stereotaxic space using a global affine transformation
\cite{collins.1994}. Afterwards, an automatic tissue-segmentation
algorithm based on a supervised artificial neural network classifier
was used to classify each voxel as cerebrospinal fluid (CSF), gray
matter, or white matter \cite{kollakian}. Subsequently a deformable surface algorithm \cite{macdonald.2000} was used to generate the outer and the inner cortical meshes.



\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{plot.gray.eps}
\end{center}\caption{Left: plot of image registration variability
$c(\alpha)$. Middle: plot of the maximum of $T$ statistic at each
scale. Right: plot of corrected P-value corresponding to the maximum
of $T$ statistic. At $\alpha=0.4$, the minimum P-value of less than
0.1 is obtained so we choose $\alpha=0.4$ to be the optimal scale
that separates the two groups. The lines are the best fitting
quadratic curve in the least squares sense.\label{fig:plot.gray}}\end{figure*}

\subsection{Stochastic Model}
Let $\mathcal{M}_o$ and $\mathcal{M}_i$ be the outer (pial) and
inner surfaces of the brain respectively. The unit sphere $S^2$ is
realized as a triangle mesh and deformed to match the outer and
inner surfaces in such a way that anatomical homology and the
topological connectivity of meshes are preserved
\cite{macdonald.2000}. The cortical surfaces can be assumed to be
smooth 2-dimensional Riemannian manifolds parameterized by two
parameters \cite{davatzikos.1995} \cite{joshi.1995}. Based on the
deformable algorithm \cite{macdonald.2000} that establishes the
homology between the $S^2$ mesh and the outer cortical surface, the
Cartesian coordinates of the mapping are discretely parameterized
by the spherical parametrization (\ref{eq:laplace2}) as $$v =
(v_1(p),v_2(p),v_3(p)).$$ The inner surface is parameterized
similarly as $$w=(w_1(p),w_2(p),w_3(p)).$$ These discrete coordinate
functions are further smoothed by the WFS: \bqn v_i(p)=\sum_{l=0}^k
\sum_{m=-l}^l e^{-l(l+1)t}f_{lm}^iY_{lm}(p). \label{eq:weighted}\eqn
We model $v_i$ stochastically as (\ref{eq:model1}) by assuming
$f_{lm}^i$ to
 follow independent normal distribution $N(\mu_{lm}^i,\sigma_{l}^2)$
 for coordinate $i$, degree $l$, and order $m$.   This assumption is
equivalent to modeling $v_i$ as the sum of signal plus
 noise:
 $$v_i(p) = \sum_{l=0}^k \sum_{m=-l}^l e^{-l(l+1)t} \mu_{lm}^i Y_{lm}(p) + \epsilon_i(p),$$
 where $\epsilon_i$ is a zero mean Guassian random field with a certain isotropic covariance
 function. A similar
stochastic modeling approach has been used in \cite{miller.1997}
where the canonical expansion of a Gaussian random field is used to
model the component of a deformation field.


 The mean and the variance functions of the surface are given by \bqn \mathbb{E} v_i(p) &=& \sum_{l=0}^k \sum_{m=-l}^l e^{-l(l+1)t}\mu_{lm}^iY_{lm}(p),\label{eq:mean.surface}\\
\mathbb{V} v_i(p) &=& \sum_{l=0}^k \sum_{m=-l}^l
e^{-2l(l+1)t}\sigma_l^2Y_{lm}^2(p)\nonumber \eqn The total
variability of the surface is then measured by
$$\int_{S^2} \mathbb{V} v_i \; d\mu(p) =  \sum_{l=0}^k \sum_{m=-l}^l e^{-2l(l+1)t}\sigma_l^2$$
indicating the increase of smoothing bandwidth decreases the total
variability. If
\bqn v_{ij}(\theta,\varphi) = \sum_{l=0}^k \sum_{m=-l}^l e^{-l(l+1)t}
 f_{lm}^{ij}Y_{lm}(\theta,\varphi)\label{eq:WFS.subject}\eqn is the WFS for
 the $j$-th subject $(1 \leq j \leq s)$, the unknown parameters $\mu_{lm}^i$ and $\sigma_l^2$
 are estimated as the sample mean and the sample variance:
\bqn\widehat{\mu_{lm}^i} &=& \frac{1}{s}\sum_{j=1}^n f_{lm}^{ij}, \label{eq:sample.mean}\\
\widehat{\sigma_{l}^2} &=&      \frac{1}{(2l+1)(s-1)}\sum_{m=-l}^l
\sum_{j=1}^s (f_{lm}^{ij} - \mu_{lm}^i)^2.\nonumber\eqn The inner
surface is stochastically modeled similarly as \bqn w_i(p) =
\sum_{l=0}^k \sum_{m=-l}^l
e^{-l(l+1)t}g_{lm}^iY_{lm}(p).\label{eq:WFS.inner}\eqn









\subsection{Surface normalization}

Previously cortical surface normalization was performed 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}
\cite{robbins.2003} \cite{thompson.1996}. 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. 


Given two WFS surfaces $v_{i1}$ and
$v_{i2}$, the displacement field $d_i$ that minimizes the integral
of the squared errors between deformed surface $d_i(v_{i1})$ and $v_{i2}$ is simply given by the following theorem.\\
\begin{theorem} \label{theorem:LSE.normal}$$ v_{i2} -v_{i1} = \arg \min_{d_i \in \mathcal{H}_k} \int_{S^2}
[v_{i1} + d_i(v_{i1})-v_{i2}]^2\; d \mu(p). $$
\end{theorem}
\begin{proof}
If $d_i \in \mathcal{H}_k$, the deformation is of the form $$v_{i1} + d_i(v_{i1})=\sum_{l=0}^k \sum_{m=-l}^l [e^{-l(l+1)t}f_{lm}^{i1}+\beta_{lm}^i]Y_{lm}$$ for unknown $\beta_{lm}^i$, which need to be determined. Substituting the term in the integral, we obtain the expression
$$I = \int_{S^2}  \Big[\sum_{l=0}^k \sum_{m=-l}^l  \alpha_{lm}^iY_{lm}\Big]^2 \;d \mu(p),$$
where
$\alpha_{lm}^i = e^{-l(l+1)t}f_{lm}^{i1}+\beta_{lm}^i - e^{-l(l+1)t}f_{lm}^{i2}$. From the orghonormality of spherical harmonics, the expression is simplified as
$$I =   \sum_{l=0}^k \sum_{m=-l}^l  \alpha_{lm}^i^2.$$
This is minimum if all $\alpha_{lm}^i$ vanishes and we obtain the result. \end{proof}



Theorem \ref{theorem:LSE.normal} shows that the optimal displacement
in the least squares sense is obtained by simply taking the
difference between two WFS representations. Unlike other surface
registration methods used in warping surfaces between
subjects \cite{chung.ni.2005} \cite{robbins.2003}
\cite{thompson.1996}, it is not necessary to consider an additional
cost function that guarantees the smoothness of the displacement
field since the displacement field $v_{i2}-v_{i1}$ is already a
linear combination of smooth basis functions. Based on this idea, we
normalize WFS surfaces. We will refer to the surface correspondence obtained by Theorem \ref{theorem:LSE.normal} as the {\em WFS-correspondence}.

Let $\bar{v_i}$ be the mean surface obtained by replacing
$\mu_{lm}^i$ in (\ref{eq:mean.surface}) with the sample mean
(\ref{eq:sample.mean}). Figure \ref{fig:multi-scale} shows the mean
surface for 24 subjects used in the study. The mean surface serves
as a template for statistical analysis later. For subject $j$, the
displacement from surface $v_{ij}$ to the template is $$ \Delta
v_{ij} = \bar{v_i} - v_{ij}.$$ Consider  the trajectory of
the deformation from $v_{ij}$ to the template $\bar v_i$ parameterized by $\alpha  \in [0,1]$:
\bqn {\bf v}_{ij}(\alpha)=v_{ij} + \alpha \Delta v_{ij} = (1-\alpha)v_{ij} +
\alpha \bar{v_i}\label{eq:scale1}.\eqn When $\alpha =0$, ${\bf
v}_{ij}(\alpha)$ is the $j$-th subject surface, while when $\alpha=1$,
it is the template. The parameter $\alpha$ controls the
amount of registration from the coarse to fine scale toward the
template. Figure \ref{fig:multi-scale} shows ${\bf v}_i(\alpha)$ at
11 different scales between 0 and 1 with 0.1 increment for a single
subject. The larger the value of $\alpha$, the smaller the image
registration variability across the subjects with respect to the
template. This is shown from the total variability computed at each
scale $\alpha$:
$$\int_{S^2} \mathbb{V}({\bf v}_{ij}) \;d\mu(p)=
c(\alpha) \sum_{l=0}^k \sum_{m=-l}^l e^{-2l(l+1)t} \sigma_{l}^2,$$
where $$c(\alpha) = \frac{n-1}{n^2}\alpha^2 +
\big(1-\frac{n-1}{n}\alpha\big)^2$$ is decreasing for $0 \leq
\alpha \leq 1$ (Figure \ref{fig:plot.gray}).

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{contour.eps}
\end{center}\caption{Left: contour plot of the average distance map
in mm at scale $\alpha=0$. Right: Gaussian kernel smoothing of the
gray matter density map with 10mm
FWHM.\label{fig:contour}}\end{figure}

\subsection{Gray Matter Density}

Gray matter density is a 3D measure defined as the probability of a particular voxel belonging to gray matter. It has been used in various anatomical studies:  normal development   \cite{good.2001}  \cite{paus.1999}, autism \cite{chung.ni.2004}, depression \cite{pizzagalli.2004}, epilepsy \cite{mcmillan.2004} and Alzheimer's disease \cite{johnson.2004} \cite{thompson.2003}. There are many different techniques for obtaining the density depending on how it is defined. In VBM, it is modeled as a Gaussian mixture on tissue intensity values \cite{ashburner.2000} \cite{good.2001}. In modulated-VBM  \cite{good.2001}, the density obtained from the standard VBM is rescaled by the Jacobian determinant of image registration to preserve the total amount of gray matter. This is related to the RAVENS (regional analysis of volumes examined in normalized space) approach \cite{davatzikos.2001}. In this study, we have avoided modulating the density  with the Jacobian determinant. Since only the surface deformation is known, we do not have any information about the Jacobian determinant inside the gray matter regions. 

Paus et al. modeled the density as a Bernoulli random variable taking value 1 inside the gray matter segmentation and 0 outside the segmented regions \cite{paus.1999}. In a slightly different formulation, Thompson computed the density as the  fraction of gray matter within a ball of radius 15mm along a cortical surface \cite{thompson.2003}. This approach is equivalent to convoluting  the binary mask of the gray matter with a uniform probability distribution of radius 15mm and interpolating voxel values to the cortical surface mesh. This equivalence relation is the basis of how we project a 3D density map to a 2D cortical surface and  compare them with cortical thickness in our study.

\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.8\linewidth]{tstat.gray.eps}
\end{center}\caption{T statistic map for gray matter density and cortical thickness thresholded at between -4 and 4 for the comparision purpose. In the gray matter density results, signals are mainly detected in either the inner surface or the outer surface but not in the middle surface. This indicates that density-based morphometries will mainly detect signal near tissue boundaries. The gray matter density and the thickness results do not overlap reconfirming that the thickness and gray matter density are not positively correlated.  
\label{fig:tstat.gray}}\end{figure*}

In this study, the gray matter density was constructed using the 3D Euclidian
distance map of the surfaces. For the outer surface
$\mathcal{M}_o$, the distance map at each
voxel $x$ is defined as
$${\tt dist}_o(x) = \min_{y \in \mathcal{M}_o}\|x-y\|,$$ where $\| \cdot \|$
is the Euclidian norm. The minimum is found using the nearest
neighbor search algorithm on an optimized k-D tree
\cite{friedman.1997}. Similarly we denote the distance map for the
inner surface $\mathcal{M}_i$ as ${\tt dist}_i(x)$. Then the
average distance map is defined as $${\tt dist}(x) = \frac{{\tt dist}_o(x) +
{\tt dist}_i(x)}{2}.$$ The average distance map for a subject is shown in
Figure \ref{fig:contour}. The minimum of the average distance is
always obtained in the middle of the outer and the inner surfaces,
where the probability of a voxel belong to the gray matter class
should be the highest. Then we define the gray matter density as
\bqn {\tt density}(x) = \exp \Big[-\frac{{\tt dist}^2(x)}{2\rho^2}\Big],
\label{eq:gray.density}\eqn where parameter $\rho^2$ controls the
spread of density ($\rho^2=3$ in our study). The gray
matter density is always between 0 and 1 and it obtains its maximum
in the interior of the gray matter region, where the average
distance map obtains the minimum. The density map is further
convoluted with the 3D Gaussian kernel $K$ with 10mm FWHM to
increase the smoothness and normality of data \cite{ashburner.2000}
\cite{chung.ni.2004} (Figure \ref{fig:contour}). The smoothed
density map $K*{\tt density}(x)$ is stochastically modeled as a
Gaussian random field. 


\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{gray.to.eps}
\end{center}\caption{Top: gray matter density projected onto
inner, mid and outer surfaces. On the inner surface, the deep sulcal
regions show the low density while the gyral ridges show high
density. On the outer surface, this is opposite. The deep sulcal
regions show high density while the gyral ridges show lower density.
The middle surface shows high density. Bottom: Scatter plot of gray
matter density over thickness. They show negative correlations.
\label{fig:gray.to.surface}}
\end{figure*}


Afterwards, the two sample t-test statistic $T$ with the equal
variance assumption is computed on the convoluted gray matter
density maps at each voxel  \cite{chung.ni.2004} . The resulting $T$
random field $T(x)$ is distributed as a student $t$ distribution
with $\nu=n_1+n_2-2$ degrees of freedom at each voxel $x$. Based on
the random field theory \cite{worsley.2005} \cite{worsley.HBM.1996}
\cite{worsley.1994}, the test statistic, which accounts for the multiple
comparison correction, is the maxima of $T$ field over the gray
matter $\mathcal{M}_{g}$. The corresponding {\em corrected P-value}
is computed using the following formula: \bq && P\Big(\sup_{x \in
\mathcal{M}_{g}} T(x) > h\Big)\\
& \approx &
\frac{\mbox{Vol}(\mathcal{M}_{g})}{\mbox{FWHM}^3}\frac{(4\ln
2)^{3/2}}{(2\pi)^2}\Big(\frac{\nu-1}{\nu}h^2
-1\Big)\Big(1+\frac{h^2}{\nu}\Big)^{-\frac{\nu-1}{2}},\eq where
$\mbox{Vol}(\mathcal{M}_{g}) = 2.13 \times 10^5 \mbox{mm}^3$ is the
volume of the gray matter of the template. The
gray matter volume is estimated by computing the volume bounded by
the outer and inner meshes \cite{chung.ni.2003}.  
Restricting the search region from the whole brain volume to the
gray matter boosts the signal detection power. 

We performed the above procedure for each image registration scale from 0 to 1 at 0.1 increment.  Increasing surface registration toward the template reduces the registration variability while increasing the gray density variability due to the misalignment of gray matter. So the maxima of $T$ field will not be obtained at the two extremes of registration scale ($\alpha=0,1$) as demonstrated in Figure \ref{fig:plot.gray}.  To increase statistical sensitivity and specificity, it is necessary to find an optimal registration scale that provides the maximal discrepancy between the groups as measured by the maxima of $T$ field. Hence, the minimum corrected
P-value can be chosen as a criteria for determining the optimal
scale. Our approach is similar to the scale-space search method \cite{worsley.1995}, where the maxima of $T$ field is searched over smoothing scale. 

The maximum $T$ statistic value and its corresponding
corrected P-value at each scale are plotted in Figure \ref{fig:plot.gray} showing that
the optimal scale is obtained when $\alpha=0.4$. At this scale, the maximum T-stat. was 5.43
while the minimum T-stat. was -5.04. The random field theory based thresholding of
$h= \pm 5.35$ gives the corrected P-value of 0.1. Figure \ref{fig:tstat.gray} shows the optimally constructed T-stat. map thresholded at $\pm 4.0$ and interpolated into the nearest point in the cortical surface showing increased gray matter density in the localized areas of the autistic subjects.






\subsection{Cortical Thickness}

The previously available approaches for computing the cortical
thickness in discrete triangle meshes produce noisy thickness
measures \cite{chung.ni.2005} \cite{fischl.2000}
\cite{macdonald.2000}. So it is necessary to smooth the thickness
measurements along the cortex via PDE based smoothing techniques \cite{andrade.2001} \cite{cachia.2003} \cite{chung.ni.2003}. On the other hand, the WFS provides smooth
functional representation of the outer and inner surfaces so that
the distance measures between the surfaces should be already smooth.Hence, the WFS avoids this additional step of thickness smoothing done in most of thickness analysis literature \cite{chung.ni.2003} \cite{chung.ni.2005} . It is not necessary to perform data smoothing in the WFS formulation.  


Using Theorem
\ref{theorem:LSE.normal} and the resulting WFS-correspondence, we establish the homology between the
outer and the inner surfaces in the least squares fashion. For the outer
surface (\ref{eq:weighted}) and the inner surface
(\ref{eq:WFS.inner}), the cortical thickness is defined to be the
Euclidean distance between the WFS-correspondence:
$${\tt thick}(p) = \Big[\sum_{l=0}^k \sum_{m=-l}^l e^{-2l(l+1)
t}(g_{lm}^i-f_{lm}^i)^2\Big]^{1/2}.$$
A similar approach has been proposed for measuring the closeness  between
two surfaces \cite{gerig.2001}; however, this is the first study using  harmonics in defining the cortical thickness. Figure \ref{fig:thick} shows the comparison of
cortical thickness computed from the traditional deformable surface
algorithm \cite{macdonald.2000} and the WFS-correspondence. The cortical thickness obtained from the traditional approach introduces a lot of triangle mesh noise into its estimation while the WFS-correspondance approach does not. The spatial smoothness of the thickness is controlled by the bandwidth $t$. 

For the group comparision between the autistic and the normal control groups, a two sample t-test was performed. 
The corresponding corrected P-value was computed using the
following formula: \bq && P\Big(\sup_{p \in S^2} T(p) > h\Big)\\
& \approx & \frac{1}{\mbox{FWHM}^2}\frac{(4\ln
2)}{(2\pi)^{1/2}}\frac{\Gamma(\frac{\nu+1}{2})}{(\frac{\nu}{2})^{1/2}
\Gamma(\frac{\nu}{2})}h\Big(1+\frac{h^2}{\nu}\Big)^{-\frac{\nu-1}{2}},\eq
where FWHM is from the heat kernel used in the WFS. We used FWHM $=0.0597$ corresponding to bandwidth $t=0.0001$. However, we did not detect any statistically significant region at $\alpha = 0.05$ level. The minimum T-stat. was -4.73 while the maximum T-stat.
was 4.83. Figure \ref{fig:tstat.gray} shows the T-stat. map thresholded at $\pm 4$. 





\subsection{Comparing Cortical Thickness and Gray Matter Density}
Most morphometric studies \cite{ashburner.2000} \cite{bookstein.2001} \cite{chung.ni.2005} 
\cite{chung.ni.2003} \cite{fischl.2000} \cite{good.2001} have never compared cortical thickness and gray matter density together so it is not clear if the two anatomical indices are positively correlated. Comparing the SPM of density and thickness in Figure \ref{fig:tstat.gray}, no statistically significant regions overlap. Since both metrics have been assumed to be the indicators of the amount of gray matter, the result is paradoxical. We have correlated these two metrics within a subject (Figure \ref{fig:gray.to.surface}) and across subjects (Figure \ref{fig:corr.gray.thick}). Surprisingly the scatter plot in Figure \ref{fig:gray.to.surface} shows negative correlation. By assigning the density value of a voxel that contains a vertex of a cortical mesh to the vertex, we can project the gray matter density onto inner, middle and outer surfaces. The middle surface is obtained by averaging the inner and the outer surfaces in the WFS-correspondance. Figure \ref{fig:corr.gray.thick} shows negatively correlated scatter plots. Figure \ref{fig:corr.gray.thick} shows the complex pattern of nonuniformity of density. On the outer surface, deep sulci have higher density compared to gyri while on the inner surface, the pattern is opposite. The middle surface shows higher density compared to the outer and the inner surfaces as expected. These complex patterns of the nonuniformity of density is due to the folding pattern of the cortex. Since the sulci on the outer surface and the gyri on the inner surface are highly folded,  these regions should have more  gray matter within the sphere of fixed radius as illustrated in Figure \ref{fig:schematic}. On the other hand, thin cortical regions will fold more than thick cortical regions.  This inverse geometric relation is causing the negative correlation between density and thickness and, in turn, the resulting SPM differ in the regions of statistically significant difference. We further computed the correlation between these two measures across 24 subjects. Figure \ref{fig:gray.to.surface} shows a similar result showing negative correlation across subjects in most regions of cortex.





\section{Conclusions and Discussions}
In this paper, we have presented a unified theoretical framework for WFS and the 
detailed numerical implementation issues. WFS is used as a smooth global parametrization of 
cortical surfaces. It is a very flexible functional estimation technique for
scalar and vector data projected onto a unit sphere. WFS is shown to be a solution of a Cauchy problem in PDE, and for a specific weights, it becomes diffusion smoothing \cite{chung.ni.2003}. As a special case of WFS when the bandwidth vanishes, the traditional SPHARM can be incorporated into this more general framework. However, WFS was shown to perform better than SPHARM when data are more noisy and discontinuous by not having the significant ringing  artifacts. 

As an application of this novel approach, 
we used WFS as a tool for comparing the gray matter and the cortical thickness in a single mathematical framework. Using the WFS representation as the ground truth, cortical thickness and gray matter density are constructed and compared. In the cortical thickness analysis, the thickness is defined using the WFS-correspondence. Afterwards, the SPM of thickness and gray matter density are compared to show the statistically significant regions do not overlap.  This surprising result is caused by the negative correlation between density and thickness. Increased folding increases the gray matter density while decreasing thickness. This should serve as a spring board for more thorough investigation on comparing cortical thickness and density based morphometric techniques such as VBM.  

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{corr.gray.thick.eps}
\end{center}\caption{Correlation of thickness and gray matter
density for 24 subjects mapped on both the inner and the outer
surfaces. Most of both inner and outer surfaces show negative
correlation. Thicker cortical regions are less convoluted so the
gray matter density tend do be lower.
\label{fig:corr.gray.thick}}\end{figure}


\section*{Acknowledgment}

The authors wish to thank Steven Robbins of the Montreal Neurological Institute, John Ashburner of the Wellcome Department of Imaging Neuroscience,  Martin A. Syner of the Department of Computer Science, University of North
Carolina at Chapel Hill, and Shubing Wang of the Department of
Statistics, University of Wisconsin-Madison for valuable discussions
and comments.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.9\linewidth]{grid.eps}
\end{center}\caption{Simple 2D schematic showing the negative correlation between thickness and
gray matter density. Gray colored pixels are gray matter. The black circle is the contour of heat kernel.  There are more gray matter pixels in region (a) than region (c) although the thickness in region (c) is thicker than that of region (a). The gray matter density in the middle of the gray matter (b) is 1 for almost all subject indicating very small between-subject and between-group variability. Because of the small between-group variability, VBM does not usually detect signal in the middle of the gray matter. Most of significant signal detected in VBM is near the tissue boundary where the between-group variability is high. 
\label{fig:schematic}}\end{figure}

\bibliography{reference.MAY.2006}
\bibliographystyle{plain}

\end{document}


\subsection{Jacobian Determinant}
In the whole brain volume tensor-based morphometry (TBM), the
Jacobian determinant of the deformation field has been used to
quantify the amount of gray matter \cite{}. Note that this 3D
Jacobian is different from 2D Jacobian that measures the surface
area expansion and shrinking along the cortex \cite{}.

We will restrict the TBM to the cortical surfaces. We will compute
the Jacobian determinant at each scale $\alpha$ in
(\ref{eq:scale1}).


 The displacement
vector from surface ${\bf v}_{ij}(\alpha)$ to the template can be
written as
$$\bar v_{i} - {\bf v}_{ij}(\alpha) = (1-\alpha)(\bar v_i -
v_{ij}) = .$$ Then the Jacobian determinant $J_j$ of warping the
$j$-th subject surface ${\bf v}_{ij}(\alpha)$ to the template is
approximated as the trace of the gradient of the displacement tensor
\cite{}:
$$J(\alpha) \approx 1 + \sum_{i=1}^3 \partial_i ( \bar
v_{i} - {\bf v}_{ij}(\alpha) ) $$


\bq \partial_1 &=& \cos \theta \cos \varphi d\theta - \sin \theta
\sin \varphi d\varphi\\
\partial_2 &=& \cos \theta \cos \varphi d\theta - \sin \theta
\cos \varphi d\varphi \\
\partial_3 &=& -\sin \theta d\theta \eq

 So we have a following general linear model (GLM):
$$K*{\tt density}(p) = \beta_0 + \beta_1 \cdot {\tt thick}(p) + \beta_2 \cdot{\tt group} +\epsilon(p)$$
where the categorical variable ${\tt group}$ is 1 for an autistic
subject and 0 for a control subject. Then we test for the density
difference $H_0: \beta_2=0$ for all $p \in \mathcal{M}_o$. This has
an effect of testing the density difference while controlling for
the thickness difference between the groups. The test statistic
$F(p)$ at each point $p$ is based on the ratio of the residual sum
of squares between the null and the alternate models
\cite{chung.ni.2004}. $F(p)$ is a $F$ random field distributed as a
$F$-distribution with 1 and $n_1+n_2-3$ degrees of freedom
\cite{worsley.HBM.1996} \cite{worsley.1994}
