% This is mymiccaipaper.tex the demonstration file of
% MICCAI 2006 for Latex2e
% adapted from LLNCS.DEM from
% the LaTeX macro package from Springer-Verlag
% for Lecture Notes in Computer Science,
% version 2.2 for LaTeX2e
%
\documentclass{llncs}
%
%\usepackage{makeidx}  % allows for indexgeneration
%
%\usepackage{graphicx}
%\graphicspath{{pics/}{figs/}}

\usepackage{epsfig,epsf,psfrag}
\usepackage{amsmath,amsfonts}
% list here all the paths to your figure folders

%\usepackage{times}
%\usepackage{epsfig}
%\usepackage{graphicx}
%\usepackage{amsmath}
%\usepackage{amssymb}

\newcommand{\ba}{\left( \begin{array}}
\newcommand{\ea}{\end{array} \right)}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}
\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}

\begin{document}

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

\vspace{2cm}
\begin{center}
\Large
{TECHNICAL REPORT NO. 1118}\\
\large March 14, 2006

\vspace{0.5cm} {\Large \bf Weighted Spherical Harmonic Representation\\
and Its Application to Cortical Analysis}
% Make title bold and 14 pt font (Latex default is non-bold, 16pt)

\vspace{1cm}

Moo K. Chung$^{1}$, Li Shen$^2$, Kim M. Dalton$^1$, Daniel J. Kelley$^1$, Steven M. Robbins$^3$, Alan C. Evans$^3$, Richard J. Davidson$^{1}$\\

\vspace{0.5cm}
%$^1$Department of Statistics, Biostatistics, and Medical Informatics\\
$^1$Waisman Laboratory for Brain Imaging and Behavior\\
%$^3$Departmnet of Psychology and Psychiatry\\
University of Wisconsin-Madison\\
Madison, WI 53706. USA\\
%Computer and Information Science Department\\
$^2$University of Massachusetts-Dartmouth, MA 02747 USA\\
$^3$Montreal Neurological Institute, McGill University, Canada\\
\vspace{0.3cm} \tt{mchung@stat.wisc.edu}}
\end{center}
%\thispagestyle{empty}
\newpage
%
\frontmatter          % for the preliminaries
%
%\pagestyle{headings}  % switches on printing of running heads
%\pagestyle{empty}  % switches off printing of running heads
%
\mainmatter              % start of your contributions
%
\title{ Weighted Spherical Harmonic Representation and Its Application to Cortical
Analysis}
%
%\titlerunning{Short Title}  % abbreviated title (for running head)
%

\author{Moo K. Chung$^{1}$, Li Shen$^2$, Kim M. Dalton$^1$, Daniel J. Kelley$^1$, Steven M. Robbins$^3$, Alan C. Evans$^3$, Richard J. Davidson$^{1}$\\
%$^1$Department of Statistics, Biostatistics, and Medical Informatics\\
$^1$Waisman Laboratory for Brain Imaging and Behavior\\
%$^3$Departmnet of Psychology and Psychiatry\\
University of Wisconsin-Madison\\
Madison, WI 53706. USA\\
%Computer and Information Science Department\\
$^2$University of Massachusetts-Dartmouth, MA 02747 USA\\
$^3$Montreal Neurological Institute, McGill University, Canada\\
\tt{mchung@stat.wisc.edu}} \institute{ }
%\author{Moo K. Chung\inst{1,2,3} \and Li Shen\inst{4}  \and Daniel Kelley\inst{3} }
%\authorrunning{R. Larsen et al.}   % abbreviated author list (for running head)
%%%% modified list of authors for the TOC (add the affiliations)
%\tocauthor{Rasmus Larsen (Technical University of Denmark), author2
%(Affiliation of Author2), author3 (Affiliation of Author3),}
%
%\institute{Informatics and Mathematical Modelling, Technical
%University of Denmark, Denmark, \email{rl@imm.dtu.dk}
%  \and
%  Department of Radiology, University Hospital, Denmark,
%  \and
% Innovative Solutions Corp., USA}

% The following lines are used to remove authors and affiliations from
% the front page and running titles
%\author{blinded}
%\authorrunning{-}
%\tocauthor{-}
%\institute{-}

\maketitle              % typeset the title of the contribution



\begin{abstract} We present a novel weighted spherical harmonic (SPHARM) representation of cortical
surfaces and its application to a cortical thickness analysis in
autism. The weighted-SPHARM is a hierarchical smoothing technique
given as the solution to a parabolic partial differential equation.
The weighted-SPHARM generalizes the classical-SPHARM with an
additional parameter that modulates the high frequency content of
data. We introduce a new algorithm called the iterative residual
fitting (IRF) and address the problem of determining the optimal
degree of the weighted-SPHARM. As an illustration, our unified
framework has been applied in detecting the regions of abnormal
behavior-structure correlation in a group of autistic subjects.
\end{abstract}

\section{Introduction}
\label{sec:intro} \vspace{-0.3cm}
 The {\em cortical thickness} has been widely used as an anatomical index for
quantifying the amount of gray matter in the brain. Magnetic
resonance images (MRI) are segmented into three tissue types:
cerebrospinal fluid (CSF), grey matter, and white matter. The
CSF/grey matter interface is called the outer cortical surface while
the grey/white matter interface is called the inner cortical
surface. The distance between the outer and the inner surfaces is
the cortical thickness. The cortical surfaces are represented 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
estimated by computing the shortest distance between vertices of the
two triangle meshes  \cite{fischl.2000} \cite{macdonald.2000}. The
mesh construction, discrete thickness computation procedures
introduce substantial noise in the thickness measure (Figure
\ref{fig:thickness}). So it is necessary to increase the
signal-to-noise ratio (SNR) and smoothness of data for the
subsequent statistical analysis. For smoothing cortical data,
diffusion equation based methods have been used \cite{andrade.2001}
\cite{cachia.2003b} \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 image processing step on top of the cortical thickness
estimation. In this paper, we will present a more direct novel
approach of representing the cortical surfaces using the weighted
spherical harmonics (SPHARM). Since the cortical surfaces are
represented as a weighted linear combination of smooth basis
functions, the resulting cortical thickness measurements are smooth
bypassing the data smoothing problem. The weighed-SPHARM differs
from the traditional SPHARM
\cite{gerig.2001}\cite{gu.2004}\cite{shen04} as a more general
representation that is formulated as the solution to a parabolic
partial differential equation (PDE). We will refer the traditional
SPHARM as the classical-SPHARM to distinguish it from our new
approach.

The theoretical construction and the numerical implementation issues
are explained in the sections \ref{sec:spharm} and
\ref{sec:numerical} respectively. Then the technique is applied in
detecting the regions of abnormal structure and behavior correlation
in a group of autistic subjects in the section \ref{sec:autism}.

\vspace{-0.1cm}
\section{Weighted spherical harmonic representation}
 \label{sec:spharm} \vspace{-0.3cm}
 For the parametrization $ (u_1,u_2,u_3)=(\sin \theta \cos \varphi, \sin \theta \sin \varphi,
\cos \theta)\label{eq:laplace}$ of the unit sphere $S^2$ with
$p=(\theta, \varphi) \in \mathcal{N} =[0,\pi]\otimes [0,2\pi)$, the
corresponding spherical Laplacian 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}.$ Let $\Delta^{2s+1}$ be the $(2s+1)$-th iterated spherical Laplacian.
 The even orders are not considered
for the convergence of the subsequent weighted-SPHARM.  The
eigenvalues and eigenfunctions of the iterated spherical Laplacian
satisfies $\Delta^{2s+1} Y_{lm} + \lambda_l^{2s+1} Y_{lm} =0.$ There
are $2l+1$ eigenfunctions $Y_{lm} ( -l \leq m \leq l)$,
corresponding to the same eigenvalue $\lambda_l^{2s+1}
=[l(l+1)]^{2s+1}$ \cite{wahba.1990}. $Y_{lm}$ is called the {\em
spherical harmonic} of degree $l$ and order $m$  and given
explicitly in
\cite{gerig.2001}\cite{gu.2004}\cite{shen04}\cite{wahba.1990}.

For $f,h \in L^2(S^2)$, the space of square integrable functions in
$S^2$, the inner product is defined as $\langle f, h \rangle =
\int_{0}^{2\pi} \int_0^{\pi} f(\theta,\varphi) h(\theta,\varphi)
\;d\mu(\theta,\varphi),$ where $d\mu(\theta,\varphi)=\sin \theta
\;d\theta \;d\varphi$. The norm is defined as $\|f\|=\sqrt {\langle
f, f \rangle}$. With respect to the inner product, the spherical
harmonics form a complete orthonormal basis in $L^2(S^2)$. So for
any $f \in L^2(S^2)$, we have the Fourier series expansion of $f$.
This is stated as Theorem 1. Let $\mathcal{H}_k = \{ \sum_{l=0}^k
\sum_{m=-l}^l\beta_{lm} Y_{lm}: \beta_{lm} \in \mathbb{R}\}$ be the
subspace spanned by up to $k$-th degree spherical harmonics.
\vspace{0.2cm}
\\{\bf Theorem 1.} $\sum_{l=0}^{k} \sum_{m=-l}^l f_{lm}Y_{lm} =
\arg\min_{h \in \mathcal{H}_k} \|f-h\|^2,$ where
$f_{lm} =\langle f, Y_{lm} \rangle$ are the Fourier coefficients.\\
This is the basis of the classical-SPHARM representation for
anatomical boundaries that are topologically equivalent to a unit
sphere \cite{gerig.2001} \cite{gu.2004} \cite{shen04}. An anatomical
boundary is mapped to a unit sphere and its Cartesian coordinates
$v(\theta,\varphi) = (v_1,v_2,v_2)$ are represented by SPHARM in
component wise fashion as $v_i =\sum_{l=0}^{k} \sum_{m=-l}^l
f_{lm}^iY_{lm}$. As an alternative to this classical-SPHARM, we
present the weighted version of the SPHARM.
\begin{figure}[t]
  \begin{center}
\includegraphics[width=1\linewidth] {MSE3.eps} &
\caption{ \label{fig:MSE} Plots of the square-root of the mean SSE
of the weighted-SPHARM for varying $t \;(0.01,0.001, 0.0001, 0)$.
$t=0$ is the classical-SPHARM. The cortical surfaces correspond to
the weighted-SPHARM at $85$-th order for different $t$. As $t \to
0$, the weighed-SPHARM converges to the
classical-SPHARM.}\end{center}
\end{figure}
The weighted-SPHARM can be formulated in a PDE-based image smoothing
framework. As a way to smooth data in $S^2$, we are interested in
the following parabolic-PDE: \bqn
\partial_t g + \Delta^{2s+1} g =0,\;
g(p,t=0)=f(p),\label{eq:cauchy}\eqn where $f$ is the given
functional data to be smoothed. Increasing the order of the iterated
Laplacian has the effect of attenuating the high frequency content
of the data faster. When $s=0$, we have the usual isotropic heat
diffusion. The solution to this Cauchy problem is solved by the
eigenfunction expansion technique and presented as the following
theorem:
%For $g \in L^2(S^2)$, the solution can be written as
%$$g(p,t) = \sum_{l=0}^{\infty} \sum_{m=-l}^l c_{lm}(t)Y_{lm}(p).$$
%Substituting it into (\ref{eq:cauchy}), we have
%$$\frac{d c_{lm}(t)}{dt} + \lambda_l^{2s+1}c_{lm}(t)=0.$$
%Solving the ordinary differential equation, we obtain $c_{lm}(t) =
%b_{lm}e^{-\lambda_l^{2s+1} t}$. By letting $t=0$ and solving for the
%initial condition, we obtain $b_{lm} = f_{lm} \rangle$.
\vspace{0.2cm}
\\{\bf Theorem 2.} The unique solution to the Cauchy problem (\ref{eq:cauchy})
 is given by
\bqn g(p,t) = \sum_{l=0}^{\infty} \sum_{m=-l}^l e^{-\lambda_l^{2s+1}
t} f_{lm}Y_{lm}(p). \label{eq:weighed}\eqn We will refer the finite
expansion $\sum_{l=0}^{k} \sum_{m=-l}^l e^{-\lambda_l^{2s+1} t}
f_{lm}Y_{lm}(p)$ as the weighted-SPHARM. When $t=0$, the
weighted-SPHARM collapses to the classical-SPHARM. The parameter $t$
is a smoothing parameter that determines the amount of smoothing.
This can be seen by reformulating the weighted-SPHARM as kernel
smoothing by rearranging the summation and the integral as $g(p,t)
%&=& \int_{S^2} f(q) \sum_{l=0}^{\infty} \sum_{m=-l}^l
%e^{-\lambda_l^{2s+1} t} Y_{lm}(q)Y_{lm}(p) \;d\mu(q).
%\label{eq:kernel}\\
= \int_{S^2} f(q) K_t(p,q)\;d\mu(q),$ where the positive symmetric
kernel is rewritten as \bq K_t(p,q) = \sum_{l=0}^{\infty}
\sum_{m=-l}^{l} e^{-\lambda_l^{2s+1}t}Y_{lm}(p)Y_{lm}(q). =
\frac{1}{4\pi}\sum_{l=0}^{\infty} (2l+1)
e^{-\lambda_l^{2s+1}t}P_l^0(\cos p \cdot q).\eq The simplification
from the double summation to the single summation is due to the
harmonic addition theorem \cite{wahba.1990}. The kernel $K_t$
obtains its maximum when $p \cdot q = 0$. The {\em full width at the
half maximum} (FWHM) of the kernel, needed in the statistical
analysis later, is obtained numerically by solving for $p \cdot q$
in
$$  \frac{1}{2}\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-\lambda_l^{2l+1}t} =
\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-\lambda_l^{2s+1}t}P_l^0(p \cdot
q).$$ Then FWHM is $2(p\cdot q)$, which is a nonlinear function of
$t$. For instance, when $s=0$ and $t=0.001$, FWHM is 0.2262.




%This can be stated as the following theorem.\\
%At each fixed point $p$, smooth representation $h$ of functional
%data $f$ is searched in the function space $L^2(\mathcal{M})$ that
%minimizes the integral of the weighted squared distance between $f$
%and $h$. This is formulated as the following regularization problem:
%\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

%The minimizer of (\ref{eq:min2}) is given by the following theorem.\\
%{\bf Theorem 2.} \bq \frac{K_t*f(p)}{\int_{\mathcal{M}}
%K_t(p,q)\;d\mu(q)} =  \arg\min_{h \in L^2(\mathcal{M})}
%\int_{\mathcal{M}} K_t(p,q)\big[f(q) - h(p)\big]^2 \; d\mu(q).\eq
%{\em Proof.} Since the integral is quadratic in $h$, the minimum
%exists and obtained by $$\frac{\partial}{\partial h}
%\int_{\mathcal{M}} K_t(p,q)[f(q) - h(p)]^2 \;d\mu(q) = -2
%\int_{\mathcal{M}} K_t(p,q)[f(q) - h(p)] \;d\mu(q)=0.$$ Solving
%this, we obtain the result.
\begin{figure}[t]
  \begin{center}
\includegraphics[width=0.9\linewidth] {model.eps} &
\caption{\label{fig:optimal} Plots showing the P-value of the
optimal degree selection procedure for different $t$. The iterative
procedure stops when the P-value $\geq 0.01$. The outer cortical
surfaces are the corresponding optimal weighted-SPHARM. The
optimal-SPHARM is obtained at $k=18 \;(t=0.01), k=42 \;(t=0.001),
k=52\; (t=0.0005), k=78\; (t=0.0001)$.}\end{center}
\end{figure}
\vspace{-0.1cm}
\section{Numerical Implementation}
\label{sec:numerical} \vspace{-0.3cm} The previous numerical
approach for constructing SPHARM has been mainly based on solving a
system of linear equations called the {\em normal equations}
\cite{gerig.2001} \cite{shen04}. Theorem 1 shows that the SPHARM
coefficients can be numerically estimated via the least squares
estimation technique as follows. At each point $p_i$ in $S^2$, we
have a normal equation $f(p_i) = \sum_{l=0}^{k} \sum_{m=-l }^{l}
f_{lm} Y_{lm}(p_i).$ This is given as the following matrix form:
$$\underbrace{\left(%
\begin{array}{c}
  f(p_1) \\
%  f(p_2)\\
  \vdots \\
f(p_n) \\
\end{array}%
\right)}_{\bf F}= \underbrace{\left(%
\begin{array}{ccccccc}
 Y_{0,0}(p_1) & Y_{1,-1}(p_1) & Y_{1,0}(p_1) &\cdots& Y_{k,-k}(p_1)&\cdots & Y_{k,k}(p_1)\\
%Y_{00}(p_2) & Y_{1,-1}(p_2) & Y_{1,0}(p_2) & \cdots & Y_{kk}(p_2)\\
  \vdots  & \vdots  & \vdots & \ddots & \vdots &\ddots & \vdots  \\
Y_{0,0}(p_n) & Y_{1,-1}(p_n) & Y_{1,0}(p_n) &\cdots&Y_{k,-k}(p_n) &\cdots & Y_{k,k}(p_n)\\
\end{array} \right)}_{{\bf Y}_0 \hspace{2cm}  {\bf Y}_1 \hspace{2cm}  \cdots \hspace{1cm} {\bf Y}_k}\underbrace{\left(%
\begin{array}{c}
  f_{00} \\
%  \beta_{1-1} \\
  \vdots \\
  f_{kk} \\
\end{array}\right)}_{\bf b}.$$
Then the least squares estimation of the Fourier coefficients
$f_{lm}$ is given by $\widehat{\bf b}= ({\bf Y}'{\bf Y})^{-1}{\bf
Y}'{\bf F},$ where ${\bf Y} = [{\bf Y}_0 {\bf Y}_1 \cdots {\bf
Y}_k]$ and ${\bf Y}_j$ is the submatrix consisting of the $k$-th
degree harmonics only. The problem with this widely used formulation
is that the size of the matrix ${\bf Y}$ is $n \times (k+1)^2$ which
can possibly reach the RAM memory limit of the most desktop
computers for large $n$ and $k$. This is mainly true for many
cortical surface extraction algorithms  that produces no less than
$n> 100,000$ nodes \cite{fischl.2000}. The problem with the
extremely large matrix can be overcome by decomposing the subspace
$\mathcal{H}_k$ into smaller subspaces. We will refer our new
algorithm as the {\em iterative residual fitting} (IRF). Since the
classical-SPHARM is a special case of the weighted version, we will
only describe the weighed version.



\vspace{-0.2cm}
\subsection{Iterative residual fitting algorithm}
\vspace{-0.1cm} 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 $\mathcal{I}_j =
\{ \sum_{m=-j}^j \beta_{jm} Y_{jm}(p): \beta_{lm} \in \mathbb{R} \}$
is the subspace spanned by the $l$-th degree spherical harmonics
only. Then at given degree $j$, we estimate the $2j+1$ Fourier
coefficients $f_{j,-j},\cdots,f_{j,j}$ simultaneously within the
smaller subspace $\mathcal{I}_j$. Instead of fitting the original
data $f$, we fit the the $(j-1)$-th degree residual defined as ${\bf
e}_{j-1} = f -\sum_{l=-(j-1)}^{j-1} \sum_{m=-l}^l
e^{-\lambda_l^{2s+1}t}f_{lm} Y_{lm}$ in the normal equation. This
can be formulated as \vspace{0.2cm}
\\{\bf Theorem 3.}
$\sum_{m=-j}^{j} f_{jm} Y_{jm} = \min_{h \in \mathcal{I}_{j}} \|{\bf
e}_{j-1}-h\|^2.$\\ Based on Theorem 3, we have the following
iterative algorithm that hierarchically constructs the
weighted-SPHARM from lower to higher degrees.


\vspace{0.3cm} \hspace{-0.5cm}{\bf Algorithm}: Iterative Residual
Fitting (IRF)\\
{\tt  Let $j=0$ and $f_{00} \leftarrow \frac{\sum_{i=1}^n
f(p_i)Y_{00}(p_i)}{\sum_{i=1}^n
 Y_{00}^2(p_i)}.$\\
While $j \leq k$ do}\\
\vspace{-0.4cm} {\tt \begin{enumerate}
$j \leftarrow j+1$.\\
${\bf e}_{j-1} \leftarrow f -
\sum_{l=-(j-1)}^{j-1} \sum_{m=-l}^l e^{\lambda_l^{2s+1}t}f_{lm} Y_{lm}.$\\
\hspace{0.7cm} $(f_{j,-j},\cdots,f_{j,j})' \leftarrow ({\bf Y}_j
{\bf Y}_j)^{-1}{\bf Y}_j({\bf e}_{j-1}(p_1),\cdots {\bf
e}_{j-1}(p_n))'$
\end{enumerate}}


Let the sum of squared error(SSE) corresponding to the $j$-th degree
be $\mbox{SSE}_j = \sum_{i=1}^n {\bf e}_j^2(p_i).$
Figure~\ref{fig:MSE} shows the plot of the square-root of
 the mean SSE given by $\sqrt{\mbox{SSE}_j/n}$ for $5 \leq j \leq
 85$ for various $t$. As $t$ decreases, the weighted-SPHARM converges to the
classical-SPHARM ($t=0$). At $t=0.0001$, the difference is
negligible.

\begin{figure}[t]
  \begin{center}
\includegraphics[width=0.9\linewidth] {thickness3.eps} &
\caption{\label{fig:thickness} Cortical thickness projected onto the
average outer cortex for various $t$ and corresponding optimal
degree: $k=18 (t=0.01), k=42 (t=0.001), k=52 (t=0.0005), k=78
(t=0.0001)$. The average cortex is constructed by averaging the
coefficients of the weighted-SPHARM. The highly noise first image
shows thickness measurements obtained by computing the distance
between two triangle meshes [blinded].}\end{center}
\end{figure}
\vspace{-0.2cm}
\subsection{Optimal weighted-SPHARM}
\vspace{-0.1cm} Although increasing the degree of SPHARM increases
the goodness-of-fit, it also increases the number of parameters 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. Consider the following $(k-1)$-th degree model
$$f(p_i) = \sum_{l=0}^{k-1} \sum_{m=-l}^l e^{-\lambda_l^{2s+1}t}f_{lm}Y_{lm}(p_i) +
\epsilon(p_i),\; i=1,\cdots, n$$ where $\epsilon$ is Gaussian random
fields. Testing if the $k$-th degree model is better than the
previous $(k-1)$-th degree model can be done by testing $H_0: f_{km}
= 0 \mbox{ for all } $-k \leq m \leq k.$ Then under the null
hypothesis, 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 $0.01$
(Figure \ref{fig:optimal}).





\begin{figure}[t]
  \begin{center}
\includegraphics[width=0.95\linewidth]{spmplot2.eps}
\caption{ \label{fig:spm} {\bf Top}: scatter plots of reaction time
(vertical) vs. thickness (horizontal). The plots corresponded to the
circled regions in the partial correlation map in the middle. {\bf
Middle}: Partial correlation maps for each group projected onto the
average outer cortex {\bf Bottom}: P-value map showing the regions
of abnormal behavior-structure correlation in autism.}\end{center}
\end{figure}

\vspace{-0.1cm}
\section{Autism structure and behavior correlation}
\label{sec:autism} \vspace{-0.3cm} We have applied the
weighted-SPHARM in detecting the regions of the abnormal abnormal
cortical thickness and behavior correlation in a group of
autistic subjects. We briefly describe how the analysis is done.\\
{\em Subjects.} 12 high functioning autistic (HFA) and 12 normal
control (NC) subjects used in this study were screened to be
right-handed males. Age distributions for HFA and NC are $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.\\
{\em Cortical thickness.} Both the outer and inner cortical surfaces
were extracted for each subject via a deformable surface algorithm
[blinded]. 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 pattern of gyral ridges are matched smoothly [blinded]. The
weighted-SPHARM with $s=0, t=0.001$ is used for representing both
surfaces with the optimal degree determined to be $k=42$. Then the
cortical thickness was computed for each subject. If
$v_i=\sum_{l=0}^k \sum_{m=-l}^l e^{-\lambda_l t}f_{lm}^iY_{lm}$ are
the coordinates for the outer surface and $w_i =\sum_{l=0}^k
\sum_{m=-l}^l e^{-\lambda_l t}g_{lm}^iY_{lm}$ are the coordinates
for the inner surface, the distance between two surface is measured
by $[\sum_{l=0}^k \sum_{m=-l}^l e^{-\lambda_l
t}(g_{lm}^i-f_{lm}^i)^2]^{1/2}$ \cite{gerig.2001} [blinded]. Figure
\ref{fig:thickness} shows the thickness projected
onto the average outer cortex.\\
{\em Behavioral measure.} The subjects were asked to decide whether
a picture of a human face was either emotional (happiness, fear or
anger) or neutral (showing no obvious emotion) by pressing one of
two buttons. The faces were black and white photographs taken from
the Karolinska Directed Emotional Faces set. The response time (ms)
for HFA and NC are $1329.8 \pm 206.7$ and $1110.9 \pm 182.3$ and
respectively. The response time has been previously used as a
behavioral index for characterizing autism [blinded].
\\ {\em Partial correlation
maps.} In correlating the behavioral measure with the cortical
thickness, we removed the confounding effect of age and the total
outer cortical surface area using the partial correlation mapping
technique. Let $Y=(Y_1,Y_2)$ be cortical thickness and response
time, and $X=(X_1,X_2)$ be age and total cortical area respectively.
The covariance matrix of $(Y,X)'$ is denoted as $\left(%
\begin{array}{cc}
  \Sigma_{YY} & \Sigma_{YX} \\
  \Sigma_{XY} & \Sigma_{XX} \\
\end{array}\right).$  The partial covariance of
$Y$ given $X$ is $\Sigma_{YY} -
\Sigma_{YX}\Sigma_{XX}^{-1}\Sigma_{XY} =(\sigma_{ij}).$ The {\em
partial correlation}
$\rho_{Y_i,Y_j|X}=\sigma_{ij}/\sqrt{\sigma_{ii}\sigma_{jj}}$
measures the correlation between variables $Y_i$ and $Y_j$ while
removing the effect of variables $X$. If we denote the partial
correlations as $\rho_1$ and $\rho_2$ for HFA and NC respectively,
we test if $H_0: \rho_1(p) = \rho_2(p) \mbox{ for all } p \in S^2$.
Based on the Fisher transform, the test statistic under $H_0$ is $
Z(p) =\frac{1}{2}\ln
\Big(\frac{1+r_1}{1-r_1}\cdot\frac{1-r_2}{1+r_2}
\Big)/\sqrt{\frac{1}{n_1-4}+\frac{1}{n_2-4}},$ where $r_i$ is the
sample partial correlation for group $i$. $Z(p)$ is an approximate
Gaussian random field. The type I error $\alpha$ for testing one
sided test is then given by: $\alpha = P\Big(\sup_{p \in
\partial \Omega} Z(p) > h \Big) \approx \sum_{d=0}^2 \phi_d(\partial
\Omega)\rho_d(h),$ where $\phi_d$ are the $d$-dimensional Minkowski
functionals of $\partial \Omega$ and $\rho_d$ are the
$d$-dimensional Euler characteristic (EC) density of correlation
field \cite{worsley.1994}. The Minkowski functionals are $\phi_0
=2,\phi_1 =0, \phi_2 = \mbox{area}(S^2)/2 = 2\pi.$ The EC densities
are: $\rho_0(h) = \int_{h}^{\infty}
\frac{1}{\sqrt{2\pi}}e^{-u^2/2}\;du$ and $\rho_2(h) = \frac{4\ln
2he^{-h^2/2}}{(2\pi)^{3/2}\mbox{FWHM}^2},$ where FWHM is computed
numerically to be 0.2262 for $s=0$ and $t=0.001$. The resulting
P-value map is given in Figure \ref{fig:spm} showing the regions of
statistically significant abnormal behavior-structure correlation in
autism.



\vspace{-0.3cm}
\section{Conclusion}
\label{sec:results} \vspace{-0.3cm} We have presented a new unified
cortical analysis framework using the weighted-SPHARM that extends
the traditional SPHARM. The weighed-SPHARM is formulated as the
solution to a parabolic PDE in such a way that the time parameter
controls the amount of smoothing in estimating the underlying
functional data. To address the computational issues, we developed
the IRF algorithm with the optimal degree selection technique. The
methodology has been successfully applied in localizing the cortical
regions of abnormal  behavior-structure correlation in autism.

\vspace{-0.3cm}
%\bibliographystyle{splncs}
%\bibliography{mymiccaipaper}


\begin{thebibliography}{10}



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

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

\bibitem{chung.ni.2003}
Chung, M., Worsley, K., Robbins, S., Paus, T., Taylor, J., Giedd,
J., Rapoport,
  J., Evans, A.
\newblock Deformation-based surface morphometry applied to gray matter
  deformation.
\newblock \textbf{18} (2003)  198--213

\bibitem{davatzikos.1995}
Davatzikos, C., Bryan, R.
\newblock Using a deformable surface model to obtain a shape representation of
  the cortex.
\newblock Proceedings of the IEEE International Conference on Computer Vision
  (1995)

\bibitem{fischl.2000}
Fischl, B., Dale, A.
\newblock Measuring the thickness of the human cerebral cortex from magnetic
  resonance images.
\newblock PNAS \textbf{97} (2000)  11050--11055



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

\bibitem{gu.2004}
Gu, X., Wang, Y., Chan, T., Thompson, T., S.T., Y.
\newblock Genus zero surface conformal mapping and its application to brain
  surface mapping.
\newblock IEEE Transactions on Medical Imaging \textbf{23} (2004)  1--10

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

\bibitem{shen04}
Shen, L., Ford, J., Makedon, F., Saykin, A.
\newblock surface-based approach for classification of 3d neuroanatomical
  structures.
\newblock Intelligent Data Analysis \textbf{8} (2004)

\bibitem{wahba.1990}
Wahba, G.
\newblock Spline models for observational data.
\newblock SIAM, New York (1990)

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

\end{thebibliography}


\end{document}
