
%\documentclass[journal]{IEEEtran}

\documentclass[11pt]{article}
%\usepackage{cite}      % Written by Donald Arseneau
                        % V1.6 and later of IEEEtran pre-defines the format
                        % of the cite.sty package \cite{} output to follow
                        % that of IEEE. Loading the cite package will
                        % result in citation numbers being automatically
                        % sorted and properly "ranged". i.e.,
                        % [1], [9], [2], [7], [5], [6]
                        % (without using cite.sty)
                        % will become:
                        % [1], [2], [5]--[7], [9] (using cite.sty)
                        % cite.sty's \cite will automatically add leading
                        % space, if needed. Use cite.sty's noadjust option
                        % (cite.sty V3.8 and later) if you want to turn this
                        % off. cite.sty is already installed on most LaTeX
                        % systems. The latest version can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/cite/



\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{psfrag}    % Written by Craig Barratt, Michael C. Grant,
                        % and David Carlisle
                        % This package allows you to substitute LaTeX
                        % commands for text in imported EPS graphic files.
                        % In this way, LaTeX symbols can be placed into
                        % graphics that have been generated by other
                        % applications. You must use latex->dvips->ps2pdf
                        % workflow (not direct pdf output from pdflatex) if
                        % you wish to use this capability because it works
                        % via some PostScript tricks. Alternatively, the
                        % graphics could be processed as separate files via
                        % psfrag and dvips, then converted to PDF for
                        % inclusion in the main file which uses pdflatex.
                        % Docs are in "The PSfrag System" by Michael C. Grant
                        % and David Carlisle. There is also some information
                        % about using psfrag in "Using Imported Graphics in
                        % LaTeX2e" by Keith Reckdahl which documents the
                        % graphicx package (see above). The psfrag package
                        % and documentation can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/psfrag/

%\usepackage{subfigure} % Written by Steven Douglas Cochran
                        % This package makes it easy to put subfigures
                        % in your figures. i.e., "figure 1a and 1b"
                        % Docs are in "Using Imported Graphics in LaTeX2e"
                        % by Keith Reckdahl which also documents the graphicx
                        % package (see above). subfigure.sty is already
                        % installed on most LaTeX systems. The latest version
                        % and documentation can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/subfigure/

%\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}
%\newtheorem{proof}{\bf Proof}  %for TMI
\newenvironment{proof}[1][Proof]{\begin{trivlist}
\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}
\begin{document}
%


\hspace{-0.6cm}DEPARTMENT OF STATISTICS\\
University of Wisconsin\\
1300 University Ave.\\
Madison, WI 53706

\vspace{2cm}
\begin{center}
\Large
{TECHNICAL REPORT NO. 1122}\\
\vspace{0.2cm} \large June 20, 2006

\vspace{0.5cm} {\Large Unified Cortical Surface Morphometry and Its
Application to Quantifying Amount of Gray Matter}
% Make title bold and 14 pt font (Latex default is non-bold, 16pt)

\vspace{1cm}

Moo K. Chung\vspace{0.5cm}

{\tt mchung@stat.wisc.edu}\\ {\tt http://www.stat.wisc.edu/$\sim$mchung}\\
\vspace{0.2cm}
Department of Statistics\\
Department of Biostatistics and Medical Informatics\\
Keck Laboratory for Functional Brain Imaging and Behavior\\
%$^4$Department of Psychology, $^5$Department of Psychiatry\\
University of Wisconsin-Madison\\
Madison, WI 53706\\
%$^6$Montreal Neurological Institute, McGill University\\
%Montreal, Quebec, Canada
\end{center}

%Version 2.0 Kernel definition changed with the addition of term
%$u_0=\det g^{1/2}$.

\thispagestyle{empty}
\newpage


% paper title
\title{Unified Cortical Surface Morphometry and Its Application to Quantifying 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,%\IEEEmembership{Member,~IEEE,},
        ~Kim~M.~Dalton,%~\IEEEmembership{Fellow,~OSA,}
        ~Li~Shen,\\
        ~Alan~C.~Evans,
        ~Richard~J.~Davidson}


       %\thanks{M.K. Chung, K. Dalton, R. Davidson is
       %with the University of Wisconsin, Madison. Steve Robinson and A.
       %Evans are with the Montreal Neurological Institute.}%~\IEEEmembership{Life~Fellow,~IEEE}% <-this % stops a space

% The paper headers
%\markboth{IEEE Transactions on Medical Imaging}{Shell \MakeLowercase{\textit{et al.}}:
%Comparative Analysis on Voxel-based and Tensor-based Morphometries}
%\markboth{Technical Report , Department of Statistics, University of
%Wisconsin-Madison}


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

% make the title area
\maketitle


\begin{abstract}

In quantifying the amount of gray matter in a population,
voxel-based morphometry (VBM) and cortical thickness analysis are the
two most widely used techniques. There are still many unanswered
methodological questions regarding these methods. It is not clear
what is the optimal amount of registration needed in VBM. It is assumed that the gray matter density obtained in
the VBM will be positively correlated to
cortical thickness but there is no systematic study that compared
these two disparate measurements. The gray matter density is a 3D
measurement defined in voxels in 3D whole brain volume while the
thickenss is a 1D measurement defined along the 2D cortical surface.
Due to this disparate dimensionality, it is not clear how to compare
them in a systematic fashion.

In this paper, we present a novel weighted Fourier Series (WFS)
representation for the cortical surface that enables us to address
these questions in a unified mathematical framework. The WFS
representation is a data smoothing technique formulated as a
solution of the Cauchy problem. This representation provides the
explicit smooth functional estimation of true unknown cortical
boundary as a linear combination of basis functions. Based on this
new representation, cortical thickness and gray matter density can
be defined and compared naturally. The basic theory underlying the
weighted Fourier Series representation and its numerical
implementation issues are presented in detail.

As an
illustration, this unified approach is applied in the cortical thickness analysis and the VBM in a group of autistic subjects. The thickness analysis and VBM results are compared in parallel.
\end{abstract}

%\begin{keywords}
%Heat kernel, Laplace-Beltrami operator, surface based morphometry
%\end{keywords}
% Note that keywords are not normally used for peerreview papers.
%\IEEEpeerreviewmaketitle



\section{Introduction}

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. Among many magnetic resonance
imaging (MRI) morphometric techniques proposed in quantifying the
amount of gray matter locally, voxel-based morphometry (VBM)
\cite{ashburner.2000} \cite{davatzikos.2001} \cite{wright.1995} and cortical thickness analysis
\cite{chung.ni.2005} \cite{chung.ni.2003} \cite{fischl.2000}
\cite{miller.2000} have been the two most widely used techniques so far.


Voxel-based morphometry (VBM) is a fully automated image analysis
technique allowing identification of regional differences in gray
and white matters between groups of subjects without a prior region
of interest (ROI). Although there are many variations in VBM, the
underpinning procedures for different VBM are identical. The 3D whole
brain MRI is normalized into a template then each voxel is assigned
a probability of that voxel belonging to a particular tissue class.
There are three main tissue classes: cerebrospinal fluid (CSF),
gray matter, and white matter. Afterwards the tissue probability is used in a general linear model
(GLM) \cite{friston.2002} to characterizing the amount of tissue concentration at each
voxel.

On the other hand, cortical thickness analysis requires an
additional step of segmenting the cortical surface. The CSF/gray matter interface is
called the {\em outer surface} (pial surface) while the gray/white
matter interface is called the {\em inner surface} \cite{macdonald.2000}. Then the
distance between the outer and inner surfaces is defined as the {\em
cortical thickness}. The cortical thickness 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}.

Although VBM and cortical thickness analysis are popular techniques,
in almost all studies, only a single morphometric technique is used
to quantify gray matter. Consequently, each study provides only a
small vignette of the overall picture of a population. Further,
studies using different morphometric techniques often result in
contradictory findings. Therefore it is necessary to be able to compare and the VBM to the
cortical thickness analysis directly in a unified methodological framework. Studies that effectively compare different morphometric techniques and anatomical measures will provide
more insight into understanding the population.

As a basic mathematical tool for achieving these goals, we present a
new representation technique called the weighed Fourier series (WFS)
representation. This is an explicit data smoothing technique
formulated as a solution to a partial differential
equation (PDE) and related to the spherical harmonic (SPHARM)
representation \cite{gerig.2001} \cite{shen.2004} indirectly. The basic
theoretical properties of WFS  and its numerical implementation
issues will be presented in great detail so that readers should be able to implement and possibly modify
the WFS for their own use. Based on WFS, VBM and cortical
thickness analysis will be performed in parallel and the results
will be compared and combined together.

The two major contributions of the paper is the theoretical and computational
development of WFS, and showing various applications toward VBM and cortical
thickness analysis in a unified framework.





\section{Literature Review}

We will mainly review literatures that are directly related to our
methodology and address what is our new contributions in the context of the previous literatures. Although there are many different morphometric techniques such as deformation-based morphometry (DBM) \cite{chung.ni.2003} \cite{davatzikos.1996} and tensor-based morphometry (TBM) \cite{ashburner.2000b} \cite{chung.mmbia.2006} \cite{chung.ni.2003} \cite{thompson.2003}  \cite{thompson.1996}, we will not review or address these methods.

\subsection{Voxel-based Morphometry} VBM as implemented in the statistical parametric mapping
(SPM) software (http://www.fil.ion.ucl.ac.uk/spm) starts with
normalizing each structural MRI to the standard SPM template and
segmenting it into white and gray matters, and CSF based on a Gaussian mixture model \cite{ashburner.2000} \cite{davatzikos.2001}  \cite{good.2001} \cite{wright.1995}. In a slightly different formulation, the
tissue density is generated by convoluting the binary mask of the
tissue with a 3D Gaussian kernel \cite{paus.1999}. The resulting
density maps are warped into a normalized space and compared across subjects. A modified version of VBM has been also
performed along the cortex, where a fraction of gray matter within a
ball of radius 15mm is taken as gray matter density
\cite{thompson.2003}. This 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 will project the 3D density maps to the 2D cortical surface and  compare them with the cortical thickness in our study.


VBM has been applied to 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}. Most previous studies in VBM did not compare the tissue density to cortical thickness or any other morphometric measures so our study is the first to present the comparative analysis with some counterintuitive results.

The objective of VBM is to compare regional difference in relative
tissue concentration. It is not necessary for image registration
used in VBM to match every cortical features exactly, but merely
corrects for global brain shape differences. If the image
registration was exact, all the segmented images would appear
identical and no statistically significant differences would be
detected \cite{ashburner.2000}. The amount of image registration
needed in VBM has been a contentious issue that has yet to be
addressed quantitatively \cite{ashburner.2001}  \cite{bookstein.2001} \cite{crum.2003} \cite{mehta.2003}. The WFS representation will enable us to address this issue in a unified fashion.






\subsection{Cortical Thickness Analysis}

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}. Unlike 3D whole brain volume based VBM, 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 time 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 the WFS, the cortical surfaces are estimated as a weighted linear combination of smooth basis functions so that most algebraic operations on the WFS will also be smooth. This is the first study comparing the result of VBM to a cortical thickness analysis. Taking the WFS as some sort of ground truth, gray matter density and cortical thickness are correlated both within a subject and between subjects.

\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 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{chung.mmbia.2006} \cite{gu.2004} \cite{shen.2006}. Gu et al. presented the SPHARM representation 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 functions are used in  constructing the Fourier series expansion of the mapping from the cortex to a unit sphere. So the SPHARM representation is more of an interpolation technique than a smoothing technique. On the other hand, WFS is a kernel smoothing technique given as a solution to a particular PDE. The solution to the PDE is expanded in basis functions similar to the Fourier series expansion with weights. WFS offers many advantages over previous PDE-based smoothing techniques. The PDE-based smoothing methods tend to suffer a numerical convergence problem \cite{andrade.2001} \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 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}. 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}.


The SPHARM representation will be shown to be the special case of WFS.  In SPHARM, 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 \cite{chung.mmbia.2006}. In most SPHARM literatures, the degree of the Fourier series expansion has been arbitrary 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.



\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 the 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 noises much faster than low
frequency noises. 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 method 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}.

The 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 the
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.

\subsection{Isotropic Diffusion on Unit Sphere}
Let us apply the WFS theory to a unit sphere.  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 with different normalizing
constants are used throughout the paper for its convenience for 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$. This form of WFS is called the
{\em weighted-SPHARM} and it has been used as a global differential parameterization of the cortex for tensor-based morphometry \cite{chung.mmbia.2006}. 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 satifies
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}


\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 degree $k=85$ 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=1\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 statistically 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{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 integrate 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, due to its
brute force nature of the technique, the computation is extremely
slow. The second method is based on the fast Fourier
transform (FFT)  \cite{bulow.2004} \cite{gu.2004}. The drawback of the 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 SPHARM literatures. 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 most of computer memories. 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 for WFS will be given in the following subsection. Although IRF was first introduced in
\cite{shen.2006}, the correctness of the algorithm was not given in \cite{shen.2006}. In this paper, we present Theorem \ref{theorem:IRF} that proves the correctness of the IRF for the first time. The IRF can be also used in estimating SPHARM coefficients by letting the bandwidth $t=0$ in the algorithm.

\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 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 IRF procedure. 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{figure}[t]
\begin{center}
\includegraphics[width=0.9\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.
 \label{fig:heat.kernel.shape}}
\end{figure}

\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{Automatic optimal degree selection in IRF}
The IRF algorithm hieratically build the WFS from lower to higher
degree. In most previous SPHARM literatures \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. We present a statistical framework for automatically determining optimal degree in the Algorithm \ref{algorithm:IRF}. 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 the equal
variance within the same degree. This assumption is equivalent to
modeling WFS as the sum of signal plus noise: \bqn f(p_i) =
\sum_{l=0}^{k} \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$
($\alpha=0.01$ in this study). For bandwidth $ t=0.0001$, the
optimal degree is determined to be $k=78$ (Figure \ref{fig:IRF}).



\subsection{Computing FWHM}
Since the WFS representation 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 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 $\thata = 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)\sigma} =
\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}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.2262. This is the FWHM we have used in the subsequent cortical thickness analysis.



\begin{figure*}[t]
\begin{center}
\includegraphics[width=1\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*}

\section{Unified Surface-based Morphometry via WFS}
In this section, we show how the WFS representation is used to quantify
cortical surface shape variations in a group of autistic subjects.
The WFS provides a unified framework for comparing VBM and the cortical thickness analysis in
parallel.

\subsection{Data Set}

$n_1=$ 12 high functioning autistic (HFA) and $n_2=$12 normal control
(NC) subjects were screened to be right-handed males. The autistic
subjects were diagnosed via The Autism Diagnostic Interview -
Revised (ADI-R) used 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 are 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} is used to generate the outer and the inner cortical meshes.


\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_2(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 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



\begin{figure*}[t]
\begin{center}
\includegraphics[width=1\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 amount of registration serves as a
multi-scale representation where the optimal scale for
differentiating two groups should be searched
for.\label{fig:multi-scale}}
\end{figure*}

\begin{figure*}[t]
\begin{center}
\includegraphics[width=1\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{Surface normalization}

Previously cortical surface normalization is 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.

Given surfaces $v_{i1}$ and
$v_{i2}$ as in (\ref{eq:WFS.subject}), the displacement field $d_i$ that minimizes the integral
of the squared errors of warping $v_{i1}$ to $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}
[d_i(v_{i1})-v_{i2}]^2\; d \mu(p). $$
\end{theorem}
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.

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 a 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 surface \bqn {\bf
v}_{ij}(\alpha)=v_{ij} + \alpha \Delta v_i = (1-\alpha)v_{ij} +
\alpha \bar{v_i},\label{eq:scale1}\eqn which is the trajectory of
the deformation from $v_{ij}$ to the template $\bar v_i$
parameterized by $\alpha \in [0,1]$. When $\alpha =0$, ${\bf
v}_{ij}(\alpha)$ is the $j$-th subject surface while when $\alpha=1$
it is the template surface. 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 over $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{Optimal Voxel-based Morphometry}
Based on the multi-scale representation of surface registration, the
optimal amount of registration needed in VBM is determined. The gray matter density is constructed using the 3D Euclidian
distance map of the surfaces at each scale. For the outer surface
$\mathcal{M}_o(\alpha)$ at scale $\alpha$, the distance map at each
voxel $x$ is defined as
$${\tt dist}_o(x) = \min_{y \in \mathcal{M}_o(\alpha)}\|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(\alpha)$ 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. In this paper, we used $\rho^2=3$. 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.

At each scale, we construct the smoothed density map for each
subject. Then 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. The optimal image
registration scale is determined to be the one that provides the
maximal discrepancy between the groups. Hence, the minimum corrected
P-value can be chosen as a criteria for determining the optimal
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. is 5.43
while the minimum T-stat. is -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.

\begin{figure*}[t]
\begin{center}
\includegraphics[width=0.95\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 comparison 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 validates the conjecture in VBM that the signal will
be mainly detected around the tissue boundary. The gray matter
density results and the thickness results do not overlap
reconfirming that the thickness and gray matter density are not
positive correlated measurements.  Our study directly demonstrates
that the VBM and cortical thickness analysis SPMs will not overlap.
\label{fig:tstat.gray}}\end{figure*}

\subsection{Cortical Thickness Analysis}

\begin{figure*}[t]
\begin{center}
\includegraphics[width=1\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*}

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.9\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 to be lower.
\label{fig:corr.gray.thick}}\end{figure}



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 literatures \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}, we establish the homology between the
outer and the inner surfaces in the least squares fashion. We will
call this homology as the {\em WFS-correspondence}. 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
the spherical 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 dose 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, two sample t-test is performed.
The corresponding corrected P-value is 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.2262$ corresponding to bandwidth $t=0.0001$. Then for the threshold $h=\pm 4.5$, we obtain the corrected
P-value of 0.1. The minimum T-stat. is -4.73 while the maximum T-stat.
is 4.83. Figure \ref{fig:tstat.gray} shows the T-stat. map thresholded at $\pm 4$.





\begin{figure*}
\begin{center}[t]
\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 the gray matter. The
black circles are the contour of heat kernel.  There are less gray
matter pixels in region (c) than region (a) 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 close to 1 for 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 the significant signal detected in VBM is near the tissue
boundary where the between-group variability is high.
\label{fig:schematic}}\end{figure*}

\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} perform VBM and cortical
thickness analysis separately and it is not clear if two anatomical
indices measure the same anatomical characteristic. Although both
the gray matter density and the cortical thickness are presumed to
measure the amount of gray matter, it is unclear if these two
measures 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. So we have correlated these two metrics within a
subject (Figure \ref{fig:gray.to.surface}) and across subjects
(Figure \ref{fig:corr.gray.thick}) to further investigate the
relationship. Surprisingly the scatter plot in Figure
\ref{fig:gray.to.surface} shows negative correlation within a
subject. 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. Surfaces show different patterns of a
negatively correlated scatter plots. Figure
\ref{fig:corr.gray.thick} shows 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 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. Based on this
primary result, it is clear that gray matter density and cortical
thickness should be analyzed together in a multivariate fashion
rather than analyzing separately.





\section{Conclusions and Discussions}
In this paper, we presented a unified theoretical framework for the
WFS and detailed numerical implementation issues. The WFS is used as
a smooth global parametrization of the cortex. It is a very flexible
functional estimation technique for scalar and vector data projected
onto a unit sphere. The WFS can be also viewed as a cortical data
smoothing technique. The WFS is shown to be a solution of the Cauchy
problem in PDE and for a specific weights, it becomes diffusion
smoothing \cite{chung.ni.2003}. As an application of this novel
approach, we used the WFS as a tool for comparing the VBM and the
cortical thickness analysis. Using the WFS representation as the
ground truth, cortical thickness and gray matter density are
constructed, and morphometric analysis on these indices are
performed.


In the VBM, the problem of determining optimal
amount of image registration has been addressed. It was shown
that the opimal scale is obtained somewhere between global affine
registration and the WFS-correspondance based nonlinear warping. This multi-scale VBM incorporates the convoluted nature of the gray
matter using the WFS more accurately than the previous 3D-based VBM.
The explicit mathematical representation of
the WFS-based surface-to-surface registration enabled us to
construct the trajectory of the deformation field. This trajectory
is used as a parameter for controlling the amount of image
registration in a multi-scale fashion. Then the optimal VBM is
chosen that gives the maximal discrimination between the two
clinical groups. In the cortical thickness analysis, the thickness is established using the same WFS- correspondence. Afterwards, the SPMs of the VBM and the thickness 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 investigating further on comparing the VBM and the cortical thickness analysis.


\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.


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

\begin{thebibliography}{10}

\bibitem{andrade.2001}
A.~Andrade, F~Kherif, J.~Mangin, K.J. Worsley, A.~Paradis, O.~Simon,
  S.~Dehaene, D.~Le~Bihan, and J-B. Poline.
\newblock Detection of fmri activation using cortical surface mapping.
\newblock {\em Human Brain Mapping}, 12:79--93, 2001.

\bibitem{ashburner.2000}
J.~Ashburner and K.~Friston.
\newblock Voxel-based morphometry - the methods.
\newblock {\em NeuroImage}, 11:805--821, 2000.

\bibitem{ashburner.2001}
J.~Ashburner and K.~Friston.
\newblock Why voxel-based morphometry should be used.
\newblock {\em NeuroImage}, 14:1238--1243, 2001.

\bibitem{ashburner.2000b}
J~Ashburner, C.~Good, and K.J. Friston.
\newblock Tensor based morphometry.
\newblock {\em NeuroImage}, 11S:465, 2000.

\bibitem{bookstein.2001}
F.L. Bookstein.
\newblock Voxel-based morphometry should not be used with imperfectly
  registered images.
\newblock {\em NeuroImage}, 14:1454--1462, 2001.

\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{bulow.2004}
T.~Bulow.
\newblock Spherical diffusion for 3d surface smoothing.
\newblock {\em IEEE Transactions on Pattern Analysis and Machine Intelligence},
  26:1650--1654, 2004.

\bibitem{cachia.2003}
A.~Cachia, J.-F. Mangin, Rivi{\'e}re D., D.~Papadopoulos-Orfanos,
F.~Kherif,
  I.~Bloch, and J.~R{\'e}gis.
\newblock A generic framework for parcellation of the cortical surface into
  gyri using geodesic vorono{\"\i} diagrams.
\newblock {\em Image Analysis}, 7:403--416, 2003.

\bibitem{chung.isbi.2006}
M.K. Chung.
\newblock Heat kernel smoothing on unit sphere.
\newblock In {\em Proceedings of IEEE International Symposium on Biomedical
  Imaging (ISBI)}, 2006.

\bibitem{chung.ni.2004}
M.K. Chung, K.M. Dalton, A.L. Alexander, and R.J. Davidson.
\newblock Less white matter concentration in autism: 2D voxel-based
  morphometry.
\newblock {\em NeuroImage}, 23:242--251, 2004.

\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{chung.mmbia.2006}
M.K. Chung, S.~Robbins, K.M. Dalton, S.~Wang, A.C. Evans, and R.J.
Davidson.
\newblock Tensor-based cortical morphometry via weighted spherical harmonic
  representation.
\newblock In {\em IEEE Computer Society Workshop on Mathematical Methods in
  Biomedical Image Analysis (MMBIA)}, 2006.

\bibitem{chung.ni.2003}
M.K. Chung, K.J. Worsley, S.~Robbins, T.~Paus, Taylor, J.N. J.,
Giedd, J.L.
  Rapoport, and A.C. Evans.
\newblock Deformation-based surface morphometry applied to gray matter
  deformation.
\newblock {\em NeuroImage}, 18:198--213, 2003.

\bibitem{collins.1994}
D.L. Collins, P.~Neelin, T.M. Peters, and A.C. Evans.
\newblock Automatic 3d intersubject registration of mr volumetric data in
  standardized talairach space.
\newblock {\em J. Comput. Assisted Tomogr.}, 18:192--205, 1994.

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

\bibitem{crum.2003}
W.R. Crum, L.D. Griffin, D.L.G. Hill, and D.J. Hawkes.
\newblock Zen and the art of medical image registration: correspondence,
  homology, and quality.
\newblock {\em NeuroImage}, 20:1425--1437, 2003.

\bibitem{dale.1999}
A.M. Dale and B.~Fischl.
\newblock Cortical surface-based analysis i. segmentation and surface
  reconstruction.
\newblock {\em NeuroImage}, 9:179--194, 1999.

\bibitem{dalton.2005}
K.M. Dalton, B.M. Nacewicz, T.~Johnstone, H.S. Schaefer, M.A.
Gernsbacher, H.H.
  Goldsmith, A.L. Alexander, and R.J. Davidson.
\newblock Gaze fixation and the neural circuitry of face processing in autism.
\newblock {\em Nature Neuroscience}, 8:519--526, 2005.

\bibitem{davatzikos.1995}
C.~Davatzikos and R.N. Bryan.
\newblock Using a deformable surface model to obtain a shape representation of
  the cortex.
\newblock {\em Proceedings of the IEEE International Conference on Computer
  Vision}, 9:2122--2127, 1995.

\bibitem{davatzikos.2001}
C.~Davatzikos, A.~Genc, D.~Xu, and S.M. Resnick.
\newblock Voxel-based morphometry using the ravens maps: Methods and validation
  using simulated longitudinal atrophy.
\newblock {\em NeuroImage}, 14:1361--1369, 2001.

\bibitem{davatzikos.1996}
C.~Davatzikos, M.~Vaillant, S.M. Resnick, J.L. Prince, S.~Letovsky,
and
  N~Bryan.
\newblock A computerized approach for morphological analysis of the corpus
  callosum.
\newblock {\em Journal of Computer Assisted Tomography}, 20:88--97, 1996.

\bibitem{fan}
J.~Fan and I.~Gijbels.
\newblock {\em Local Polynomial Modelling and Its Applications}.
\newblock Chapman $\&$ Hall/CRC, 1996.

\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{fischl.1999}
B.~Fischl, M.I. Sereno, R.~Tootell, and A.M. Dale.
\newblock High-resolution intersubject averaging and a coordinate system for
  the cortical surface.
\newblock {\em Hum. Brain Mapping}, 8:272--284, 1999.

\bibitem{friedman.1997}
J.H. Friedman, J.L. Bentley, and R.A. Finkel.
\newblock An algorithm for finding best matches in logarithmic expected time.
\newblock {\em ACM transactions on mathematics software}, 3:209--226, 1997.

\bibitem{friston.2002}
K.J. Friston.
\newblock A short history of statistical parametric mapping in functional
  neuroimaging.
\newblock Technical Report Technical report, Wellcome Department of Imaging
  Neuroscience, ION, UCL., London, UK., 2002.

\bibitem{gerig.2001}
G.~Gerig, M.~Styner, D.~Jones, D.~Weinberger, and J.~Lieberman.
\newblock Shape analysis of brain ventricles using spharm.
\newblock In {\em MMBIA}, pages 171--178, 2001.

\bibitem{gerig.2002}
G.~Gerig, M.~Styner, and G.~Szekely.
\newblock Statistical shape models for segmentation and structural analysis.
\newblock I:467--473, 2004.

\bibitem{good.2001}
C.D. Good, I.S. Johnsrude, J.~Ashburner, R.N.A. Henson, K.J.
Friston, and
  R.S.J. Frackowiak.
\newblock A voxel-based morphometric study of ageing in 465 normal adult human
  brains.
\newblock {\em NeuroImage}, 14:21--36, 2001.

\bibitem{groemer}
H.~Groemer.
\newblock {\em Geometric Applications of Fourier Series and Spherical
  Harmonics,}.
\newblock Cambridge University Press, 1996.

\bibitem{gu.2004}
X.~Gu, Y.L. Wang, T.F. Chan, T.M. Thompson, and Yau. S.T.
\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{johnson.2004}
S.C. Johnson, L.C. Baxter, L.~Susskind-Wilder, D.J. Connor, M.N.
Sabbagh, and
  R.J. Caselli.
\newblock Hippocampal adaptation to face repetition in healthy elderly and mild
  cognitive impairment.
\newblock {\em Neuropsychologia}, 42:980--989, 2004.

\bibitem{joshi.1995}
S.C. Joshi, J.~Wang, M.I. Miller, D.C. Van~Essen, and U.~Grenander.
\newblock On the differential geometry of the cortical surface.
\newblock {\em Vision Geometry IV}, pages 304--311, 1995.

\bibitem{kelemen}
A.~Kelemen, G.~Szekely, and G.~Gerig.
\newblock Elastic model-based segmentation of 3-d neuroradiological data sets.
\newblock {\em IEEE Transactions on Medical Imaging}, 18:828--839, 1999.

\bibitem{kollakian}
K.~Kollakian.
\newblock Performance analysis of automatic techniques for tissue
  classification in magnetic resonance images of the human brain.
\newblock Technical Report Master's thesis, Concordia University, Montreal,
  Quebec, Canada, 1996.

\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{mcmillan.2004}
A.B. McMillan, B.P. Hermann, S.C. Johnson, R.R. Hansen,
M.~Seidenberg, and M.E.
  Meyerand.
\newblock Voxel-based morphometry of unilateral temporal lobe epilepsy reveals
  abnormalities in cerebral white matter.
\newblock {\em NeuroImage}, 23:167--174, 2004.

\bibitem{mehta.2003}
S.~Mehta, T.J. Grabowski, Y.~Trivedi, and H.~Damasio.
\newblock Evaluation of voxel-based morphometry for focal lesion detection in
  individuals.
\newblock {\em NeuroImage}, 20:1438--1454., 2004.

\bibitem{miller.1997}
M.I. Miller, A.~Banerjee, G.E. Christensen, S.C. Joshi, N.~Khaneja,
  U.~Grenander, and L.~Matejic.
\newblock Statistical methods in computational anatomy.
\newblock {\em Statistical Methods in Medical Research}, 6:267--299, 1997.

\bibitem{miller.2000}
M.I. Miller, A.B. Massie, J.T. Ratnanather, K.N. Botteron, and J.G.
Csernansky.
\newblock Bayesian construction of geometrically based cortical thickness
  metrics.
\newblock {\em NeuroImage}, 12:676--687, 2000.

\bibitem{muller.2005}
H.-G. Muller.
\newblock Functional modeling and classification of longitudinal data.
\newblock {\em Scandinavian Journal of Statistics}, 32:223--240, 2005.

\bibitem{paus.1999}
T.~Paus, A.~Zijdenbos, K.J. Worsley, D.L. Collins, J.~Blumenthal,
J.N. Giedd,
  J.L. Rapoport, and A.C. Evans.
\newblock Structural maturation of neural pathways in children and adolescents:
  In vivo study.
\newblock {\em Science}, 283:1908--1911, 1999.

\bibitem{pizzagalli.2004}
D.A. Pizzagalli, T.R. Oakes, A.S. Fox, M.K. Chung, C.L. Larson, H.C.
  Abercrombie, S.M. Schaefer, R.M. Benca, and R.J. Davidson.
\newblock Functional but not structural subgenual prefrontal cortex
  abnormalities in melancholia.
\newblock {\em Molecular Psychiatry}, 9:393--405, 2004.

\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}, 2006. in press.

\bibitem{quicken}
M.~Quicken, C~Brechbuhler, J.~Hug, H.~Blattmann, and G.~Szekely.
\newblock Parameterization of closed surfaces for parametric surface
  description.
\newblock {\em IEEE Computer Society Conference on Computer Vision and Pattern
  Recognition (CVPR)}, pages 354--360, 2000.

\bibitem{ramsay.1997}
J.O. Ramsay and B.W. Silverman.
\newblock {\em Functional Data Analysis}.
\newblock Springer, 1997.

\bibitem{robbins.2003}
S.M. Robbins.
\newblock Anatomical standardization of the human brain in euclidean 3-space
  and on the cortical 2-manifold.
\newblock Technical Report PhD thesis, School of Computer Science, McGill
  University, Montreal, Quebec, Canada, 2003.

\bibitem{shen.2006}
L.~Shen and M.K. Chung.
\newblock Large-scale modeling of parametric surfaces using spherical
  harmonics.
\newblock In {\em Third International Symposium on 3D Data Processing,
  Visualization and Transmission (3DPVT)}.

\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, 2004.

\bibitem{sled}
J.G. Sled, A.P. Zijdenbos, and A.C. Evans.
\newblock A nonparametric method for automatic correction of intensity
  nonuniformity in mri data.
\newblock {\em IEEE Transactions on Medical Imaging}, 17:87--97, 1988.

\bibitem{thompson.2003}
P.M. Thompson, Hayashi K.M., de~Zubicaray~G., Janke A.L., Rose S.E.,
Semple J.,
  Herman D., Hong M.S., Dittmer S.S., Doddrell D.M., and Toga A.W.
\newblock Dynamics of gray matter loss in alzheimer's disease.
\newblock {\em J. Neurosci.}, 23:994--1005, 2003.

\bibitem{thompson.1996}
P.M. Thompson and A.W. Toga.
\newblock A surface-based technique for warping 3-dimensional images of the
  brain.
\newblock {\em IEEE Transactions on Medical Imaging}, 15, 1996.

\bibitem{wahba.1990}
G.~Wahba.
\newblock {\em Spline models for observational data}.
\newblock SIAM, 1990.

\bibitem{worsley.1994}
K.J. Worsley.
\newblock Local maxima and the expected euler characteristic of excursion sets
  of $\chi2$, $f$ and $t$ fields.
\newblock {\em Advances in Applied Probability.}, 26:13--42, 1994.

\bibitem{worsley.HBM.1996}
K.J. Worsley, S.~Marrett, P.~Neelin, A.C. Vandal, K.J. Friston, and
A.C. Evans.
\newblock A unified statistical approach for determining significant signals in
  images of cerebral activation.
\newblock {\em Human Brain Mapping}, 4:58--73, 1996.

\bibitem{worsley.2005}
K.J. Worsley, J.E. Taylor, F.~Tomaiuolo, and J.~Lerch.
\newblock Unified univariate and multivariate random field theory.
\newblock {\em NeuroImage}, 12:in press, 2005.

\bibitem{wright.1995}
I.~C. Wright, P.~K. McGuire, J.-B. Poline, J.~M. Travere, R.~M.
Murray, C.~D.
  Frith, R.~S.~J. Frackowiak, and K.~J. Friston.
\newblock A voxel-based method for the statistical analysis of gray and white
  matter density applied to schizophrenia.
\newblock {\em NeuroImage}, 2:244--252, 1995.

\end{thebibliography}

\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}
