\documentclass{llncs}
\usepackage{llncsdoc}
\usepackage{epsfig,epsf,psfrag}
\usepackage{amsmath,amsfonts}
\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}

\begin{document}

\title{Unified Statistical Approach\\ to Cortical Thickness Analysis}
%\titlerunning{Cortical Thickness Analysis}
\author{Moo K. Chung\inst{1} \and Steve Robbins\inst{2} \and Alan C. Evans\inst{2}}
\institute{Department of Statistics,
%Waisman Laboratory for Brain Imaging and Behavior\\
University of Wisconsin-Madison
%,Madison, WI 53706\\
\and Montreal Neurological Institute, McGill University, Canada\\
\email{mchung@stat.wisc.edu, stever, alan@bic.mni.mcgill.ca}}
%\authorrunning{Chung et al.}
\maketitle

\begin{abstract}
This paper presents a unified image processing and analysis
framework for cortical thickness in characterizing a clinical
population. The emphasis is placed on the development of data
smoothing and analysis framework. The human brain cortex is a highly
convoluted surface. Due to the convoluted non-Euclidean surface
geometry, data smoothing and analysis on the cortex are inherently
difficult. When measurements lie on a curved surface, it is natural
to assign kernel smoothing weights based on the geodesic distance
along the surface rather than the Euclidean distance. We present a
new data smoothing framework that address this problem implicitly
without actually computing the geodesic distance and present its
statistical properties. Afterwards, the statistical inference is
based on the random field theory based multiple comparison
correction. As an illustration, we have applied the method in
detecting the regions of abnormal cortical thickness in 16 high
functioning autistic children.
\end{abstract}

%\keyword{Cortical thickness, Heat kernel, Diffusion, Image analysis,
%Cortex}
\section{Introduction}
The human cerebral cortex has the topology of a 2D highly convoluted
grey matter shell of average thickness of 3mm. The thickness of the
grey matter shell is usually referred as the {\em cortical
thickness} and can be obtained from magnetic resonance images (MRI).
The cortical thickness can be used as an anatomical index for
quantifying cortical shape variations. The thickness measures are
obtained after a sequence of image processing steps which are
described briefly here. The first step is to classify each voxel
into three different tissue types: cerebrospinal fluid (CSF), grey
matter, and white matter. The CSF/grey matter interface is called
the {\em outer cortical surface} while the grey/white matter
interface is called the {\em inner cortical surface}. These two
surfaces bound the gray matter. The mainstream approach in
representing the cortical surface has been to use a fine triangular
mesh that is constructed from deformable surface algorithms
 \cite{davatzikos.1995} \cite{MacDonald.2000}. Cortical thickness is estimated by computing
the distance between the two triangular meshes
\cite{Jones.2000} \cite{MacDonald.2000}. %\cite{miller.2000}. \cite{fischl.2000}  \cite{dale.1999}
In our study, we have used the method presented in
\cite{MacDonald.2000}. In order to compare cortical thickness
measures across subjects, it is necessary to align the cortical
surfaces via surface registration algorithms
  \cite{Robbins.2003} \cite{thompson.toga.1996}. %\cite{dale.1999}\cite{davatzikos.1997}\cite{fischl.1999}
For cross-comparison between subjects, surfaces are registered into
the {\em template surface} which serves as reference coordinates.

The image segmentation, thickness computation and surface
registration procedures are expected to introduce noise in the
thickness measure. In order to increase the signal-to-noise ratio
(SNR) and smoothness of data for the random field theory, some type
of data smoothing is necessary. %\cite{kiebel.1996}.
For 3D whole brain MRIs, Gaussian kernel smoothing is widely used to
smooth data, in part, due to its simplicity in numerical
implementation. The Gaussian kernel weights an observation according
to its Euclidean distance. However, data residing on the convoluted
brain surface fails to be isotropic in the Euclidean sense. On the
curved surface, a straight line between two points is not the
shortest distance so one may incorrectly assign less weights to
closer observations. So when the observations lie on the cortical
surface, it is more natural to assign the weights based on the
geodesic distance along the surface. Previously {\em diffusion
smoothing} has been developed for smoothing data along the cortex
before the random field based multiple comparison correction
\cite{andrade.2001} \cite{chung.taylor.2004} \cite{chung.cvpr.2003}.
By solving a diffusion equation on a manifold, Gaussian kernel
smoothing can be indirectly generalized. Although diffusion
smoothing has been used widely in image analysis starting with
\cite{perona}, most of previous work is about surface
fairing \cite{sochen}. %\cite{malladi} \cite{tang} \cite{taubin}.
There is a very few publications that smooth out observations
defined on surface for data analysis \cite{andrade.2001}
\cite{cachia.2003b}  \cite{chung.taylor.2004}
\cite{chung.cvpr.2003}.
%\cite{memoli.2004}.
The drawback of the previous diffusion smoothing approach is the
need for setting of up a finite element method (FEM) to solve the
diffusion equation numerically and making the algorithm converges
\cite{chung.taylor.2004}. To address this problem, we have developed
a simpler and more efficient method based on the heat kernel
convolution on a manifold.

As an illustration, the method was applied to groups of autistic and
normal subjects, and we were able to detect the regions of
statistically significant cortical thickness difference between the
groups.



\begin{figure}[h]
\centering
\includegraphics[width=122mm]{simulation.eps} \caption{Top: Heat kernel
 smoothing of cortical thickness with $\sigma = 1$ and $k=20, 100, 200$ iterations.  Bottom: Heat kernel
smoothing on simulated data with $\sigma=1$ and $k=20, 200, 5000$
iterations. The mean thickness and the variance are estimated from
12 normal subject data and Gaussian white noise is added to the mean
function.} \label{fig:hk}
\end{figure}


\section{Heat kernel smoothing}
The cortical surface $\partial \Omega$ can be assumed to be a $C^2$
Riemannian manifold \cite{joshi.1995}. Let $p=X(u^1,u^2) \in
\partial \Omega$ be the parametric representation of $\partial \Omega$.
We assume the following model on thickness measure $Y$:
$$Y(p) = \theta(p) + \epsilon(p),$$
where $\theta(p)$ is a mean thickness function and $\epsilon(p)$ is
a zero-mean random field, possibly a Gaussian white noise process,
with covariance function $R_{\epsilon}(p,q)$. The Laplace-Beltrami
operator $\Delta$ corresponding to the surface parameterization
$p=X(u^1,u^2) \in
\partial \Omega$ is given by $$ \Delta = \frac{1}{\det g^{1/2}}\sum_{i,j=1}^2 \frac{\partial
}{\partial u^i}\Big(\det g^{1/2}g^{ij}\frac{\partial}{\partial
 u^j}\Big),$$
where $g=(g_{ij})$ is the Riemannian metric tensor. Solving equation
$\Delta \psi = \lambda \psi$, we order eigenvalues $0 = \lambda_0
\leq \lambda_1 \leq \lambda_2 \leq \cdots$ and corresponding
eigenfunctions $\psi_0,\psi_1,\cdots$. The eigenfunctions $\psi_j$
form orthonormal basis of $L^2(\partial \Omega)$, the $L^2$ space of
functions defined on $\partial \Omega$. On the unit sphere, the
eigenvalues are $m(m+n-1)$ and the corresponding eigenfunctions are
spherical harmonics $Y_{lm}$ ($|m| \leq l, 0\leq l$)
\cite{wahba.1990}. On an arbitrary surface, the explicit
representation of eigenvalues and eigenfunction are only obtained
through numerical methods. Based on orthonormal basis, the {\em heat
kernel} $K_{\sigma}(p,q)$ is analytically given as \bqn
K_{\sigma}(p,q) =\sum_{j=0}^{\infty}e^{-\lambda_j\sigma}
\psi_j(p)\psi_j(q)\label{eq:heatkernel},\eqn where $\sigma$ is the
bandwidth of the kernel \cite{berline} \cite{rosenberg} . When
$g_{ij} = \delta_{ij}$, the heat kernel becomes the Gaussian kernel,
which is the probability density of $N(0,\sigma^2)$.
%$$K_{\sigma}(p,q) = \frac{1}{(2\pi \sigma)^{1/2}} \exp \Big[
%-\frac{\|p-q\|^2}{2\sigma^2}\Big].$$ Hence the heat kernel is a
natural extension of the Gaussian kernel. This can be interpreted as
the transition probability density for an isotropic diffusion
process with respect to the surface area element \cite{wang.1997}.
The kernel is symmetric, i.e. $K_{\sigma}(p,q) = K_{\sigma}(q,p)$
and isotropic with respect to the geodesic distance $d(p,q)$.

%The
%property of a kernel being isotropic needs some explanation. Let
%$d(p,q)$ be the geodesic distance between $p$ and $q$. Then function
%$f$ is {\em isotropic} on surface $\partial \Omega$ if
%$f(p)=constant$ for all point $p$ on geodesic circle
%$d(0,p)=\mbox{constant}$. Since $K_{\sigma}(p,q)$ has two arguments
%while symmetric, the isotropic property holds for either one of the
%arguments.
\begin{definition}
Heat kernel smoothing of cortical thickness $Y$ is the convolution:
\bqn K_{\sigma} *Y (p)=\int_{\partial \Omega}
K_{\sigma}(p,q)Y(q)\;dq.\label{eq:conv}\eqn \end{definition} It can
be written in terms of basis function expansion: $$K_{\sigma}*Y(p) =
\sum_{j=0}^{\infty} \alpha_j \phi_j(p),$$ where $\alpha_j =
e^{-\lambda_j \sigma} \int_{\partial \Omega} \phi_j(q)Y(q)\;dq$. We
also define the {\em heat kernel estimator} of unknown signal
$\theta(p)$ to be $\hat \theta_{\sigma} (p) = K_{\sigma} * Y(p)$. As
$\sigma \to 0$, $K_{\sigma}(p,q)$ becomes the Dirac delta function
$\delta(p,q)$ so the heat kernel estimator becomes unbiased as
$\sigma \to 0$, i.e. $\lim_{\sigma \to 0} \mathbb{E} \hat
\theta_{\sigma}(p) = \theta(p)$. As $\sigma$ gets larger, the bias
increases. However the total bias over all cortex is always zero,
i.e. $\int_{\partial \Omega} [\theta(p) - \mathbb{E} \hat
\theta_{\sigma}(p)] \; dp = 0$.
%, which can be seen from
%$$\int_{\partial\Omega} \mathbb{E}
%\theta_{\sigma}(p) \; dp = \int_{\partial \Omega} \int_{\partial
%\Omega} K_{\sigma}(p,q) \theta(p) \; dq \;dp = \int_{\partial
%\Omega} \theta(p) \;dp.$$
Let us list important nontrivial properties of heat kernel
smoothing.

\begin{theorem} \label{theorem:diff} $K_{\sigma}*Y$ is the unique
solution of the following isotropic diffusion equation at time
$t=\sigma^2/2$: \bqn \frac{\partial f}{\partial t} = \Delta f,\;
f(p,0)=Y(p)\label{eq:diff}, p \in \partial \Omega\eqn
\end{theorem}
This is a well known result \cite{rosenberg}. This theorem implies
that the heat kernel smoothing isotropically assigns weights on
$\partial \Omega$.
\begin{theorem} \label{theorem:poly}
$$K_{\sigma}*Y(p) = \arg \min_{\theta(p)\in L^2(\partial \Omega)}
\int_{\partial \Omega} K_{\sigma}(p,q)\big[Y(q) -\theta(p)]^2
\;dq.$$
\end{theorem}
The proof can be found in \cite{chung.ni.2005}. This shows that the
heat kernel smoothing can be formulated as a regression on a
manifold.

\begin{figure}[t]
\centering
\includegraphics[width=120mm]{validation.eps}
\caption{Left: Within-subject variance plotted over the number of
iterations of heat kernel smoothing with $\sigma=1$. Decreasing
variance implies the convergence of the heat kernel smoothing to the
mean thickness (Theorem 4). Right: Between-subject variance plotted
over the number of iterations illustrating Theorem 5.
\label{fig:validation}}
\end{figure}



\begin{theorem}\label{theorem:iteration}
$$\underbrace{K_{\sigma}
*\cdots * K_{\sigma}}_{k  \mbox{ \small{times} }}*Y = K_{\sqrt k
\sigma} * Y.$$
\end{theorem}
This can be seen as a scale space property of diffusion. From
Theorem 1, $K_{\sigma}* (K_{\sigma}
* Y)$ can be taken as the diffusion of signal $K_{\sigma}*Y$ after
time $\sigma^2/2$ so that $K_{\sigma}* (K_{\sigma}
* Y)$ is the diffusion of signal $Y$ after time $\sigma^2$. Hence
$$K_{\sigma}*K_{\sigma}* Y = K_{\sqrt 2 \sigma} * Y.$$
Arguing inductively we see that the general statement holds. We will
denote the $k$-fold iterated kernel as $K_{\sigma}^{(k)} =
\underbrace{K_{\sigma} *\cdots * K_{\sigma}}_{k  \mbox{
\small{times} }}$. This is the basis of our iterated heat kernel
smoothing. Heat kernel with a large bandwidth will be performed by
iteratively applying heat kernel smoothing with a smaller bandwidth.
 For instance iterated heat kernel smoothing
with $\sigma=1$ and $k=200$ will generate heat kernel smoothing with
the effective bandwidth of $\sqrt{200}=14.14$mm. Figure \ref{fig:hk}
shows the process of iterated heat kernel smoothing.

\begin{theorem} $$\lim_{\sigma \to \infty} K_{\sigma}
*Y = \frac{\int_{\partial \Omega} Y(q) \; dq}{\mu(\partial
\Omega)}.$$
\end{theorem}
Here $\mu(\partial \Omega)$ is the total surface area of
$\partial\Omega$. This theorem shows that when we choose large
bandwidth, heat kernel smoothing converges to the sample mean of
data on $\partial \Omega$.
%This is proved in [blind]. %This is easily proved by noting
%$$\widetilde{K}_{\sigma}*Y(p) =
%\frac{\int_{B_p}
%\exp\big[-\frac{d^2(p,q)}{2\sigma^2}\big]Y(p)\;dq}{\int_{B_p}
%\exp\big[-\frac{d^2(p,q)}{2\sigma^2}\big]\; dq} \to
%\frac{\int_{B_p}Y(q)\;dq}{\mu(B_p)}$$ as $\sigma \to \infty$. Now by
%letting $B_p$ to cover the whole cortex $\partial \Omega$, we prove
%the property.
Figure \ref{fig:hk} (bottom) shows the convergence of heat kernel
smoothing to the within-subject mean cortex 4mm as the bandwidth
increases. Figure \ref{fig:validation} (left) shows the convergence
of the within-subject variance indirectly implying $K_{\sigma}*Y$
converges to a constant, which is the average thickness over the
cortex.


\begin{figure}[t]
\centering
\includegraphics[scale=0.15]{sulcal.pattern.eps}
\caption{Thickness maps are projected onto a unit square. Left:
original noisy thickness map. Right: Heat kernel smoothing with
 $\sigma=1$ and $k=200$ iterations.\label{fig:sulcal.pattern}}
\end{figure}


It is natural to assume the measurements $Y(p)$ and $Y(q)$ to have
less correlation when $p$ and $q$ are away so we assume the
covariance function to be $R_{\epsilon}(p,q)=\rho(d(p,q))$ for some
nondecreasing function $\rho$. Then we can show the variance
reduction property of heat kernel smoothing.
\begin{theorem}\label{theorem:shijie}
${\bf Var} [K_{\sigma}*Y(p)] \leq {\bf Var} Y(p) \mbox{ for each } p
\in \partial \Omega.$
\end{theorem}
%The proof can be found in [blind].
Figure \ref{fig:validation}
(right) shows the between-subject variance decreases as $\sigma$
increases.



\begin{figure}[t]
\centering
\includegraphics[scale=0.11]{sulcus.eps}
\caption{\label{fig:spam-cachia-labels}Automatically generated
traces of the central and superior temporal sulcal fundi
\cite{cachia.2003b}. The first column shows the traces generated for
the template surface. The second column shows the probability of
sulcal matching based on 149 normal subjects before any surface
normalization. The third column shows the probabilities after
surface normalization. The first row is the left hemisphere and the
second row is the right hemisphere. Note that the distribution is
much more spatially concentrated and the matching probabilities are
much greater after normalization.}
\end{figure}

\begin{figure}[h]
\centering
%\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.12]{matching.probability3.eps}
\caption{\label{fig:cs-spam}Probability of sulcal matching, after
normalization, for 39 manually identified central sulci defined as
the surface region surrounded by gyri, not just the fundus. The
views are illustrated on a slightly-opened version of the template
cortical surface in order to better view inside the sulcus. The
warping in 2D localizes the central sulcus nearly completely inside
the template central sulcus. Left (right) figure is the left (right)
central sulci.}
\end{figure}


The problem with the heat kernel smoothing on an arbitrary surface
is that the explicit analytic form of the heat kernel is unknown. To
address this problem we use the {\em parametrix expansion} of the
heat kernel \cite{rosenberg} \cite{wang.1997}: \bqn K_{\sigma}(p,q)
&=& \frac{1}{(2\pi\sigma)^{1/2}}\exp
\big[-\frac{d^2(p,q)}{2\sigma^2}\big][1 +
O(\sigma^2)\big]\label{eq:parametrix}\eqn for small $d(p,q)$. This
expansion spells out the exact form of the kernel for small
bandwidth. When the metric is flat, the heat kernel becomes a
Gaussian kernel, reconfirming that heat convolution is a
generalization of Gaussian kernel. The expansion is the basis of our
heat kernel smoothing formulation. Heat kernel smoothing with a
large bandwidth will be decomposed into iterated kernel smoothing.
We will truncate and normalize the heat kernel using the first order
term. For each $p \in \partial \Omega$, we define
 \bqn \widetilde{K}_{\sigma}(p,q) =
\frac{\exp\big[-\frac{d^2(p,q)}{2\sigma^2}\big]{\bf 1}_{B_p}(q)}
{\int_{B_p} \exp\big[-\frac{d^2(p,q)}{2\sigma^2}\big]\;dq}
\label{eq:kern},\eqn where ${\bf 1}_{B_p}$ is an indicator function
defined on a small compact domain containing $B$ such that ${\bf
1}_{B_p}(q) =1$ if $q \in B_p$ and ${\bf 1}_{B_p}(q) =0$ otherwise.
Note that for each fixed $p$, $\widetilde{K}_{\sigma}(p,q)$ defines
a probability distribution in $B_p$ and it converges to
$K_{\sigma}(p,q)$ as $\sigma \to 0$ in $B_p$ This implies
$$ \widetilde{K}_{\sigma}^{(k)}*Y(p) \to K_{\sigma}^{(k)}*Y(p) \mbox{ as } \sigma \to 0.$$
%We can proved this only when $g_{ij}=\sigma_{ij}$ since the exact
%analytical expression of the heat kernel is nonexistent and still an
%ongoing research problem \cite{rosenberg}  \cite{wang.1997}.
%
%\begin{theorem}\label{theorem:flat} When $g_{ij}=\delta_{ij}$,
%$$K_{\sqrt{k}\sigma}*Y(p) \leq
%\widetilde{K}_{\sigma}^{(k)}*Y(p) \leq \frac{1}{\alpha^k(p)}
%K_{\sqrt{k}\sigma}*Y(p),$$ where $\alpha(p) = \int_{B_p}
%K_{\sigma}(p,q)\;dq$.
%\end{theorem}
%The upper bound is still valid without the assumption of the flat
%metric. As $\sigma \to 0$, $\alpha(B_p) \to 1$ and the inequalities
%collapse proving the convergence. For sufficiently small $\sigma$
%and $B_p$, we can make $\alpha(B_p)$ as small as possible.  Hence,
%decreasing the size of bandwidth and increasing the number of
%iterations will perform better although it will be slow down the
%algorithm.

For a discrete triangular mesh, we can take $B_p$ to be a set of
points containing $p$ and its neighboring nodes $q_1,\cdots,q_m$,
and take a discrete measure on $B_p$, which still make
(\ref{eq:kern}) a probability distribution. This can be viewed as a
{\em Gaussian kernel Nadaraya-Watson} type smoothing extended to
manifolds \cite{chaudhuri.2000}. %\cite{fan}.
Figure
\ref{fig:sulcal.pattern} shows a flattened thickness map
illustrating how heat kernel smoothing can enhance the thickness
pattern by increasing the signal-to-noise ratio.

%Figure \ref{fig:hk}
%shows the heat kernel smoothing on both a single subject data and
%simulated data. The cortical thickness measures are projected onto a
%template that has less surface folding to show the progress of
%smoothing. The bottom figure shows the progress of heat kernel
%smoothing with $\sigma=1$ and upto $k=5000$ iterations applied to a
%simulated data.
%From 12 normal subjects, the mean thickness $\theta(p)$ and variance
%$R_{\epsilon}(p,p)$ functions are estimated and used to generate
%random fields with a Gaussian white noise. It begin to show the
%convergence to the sample mean thickness over all cortex.


%\begin{figure}
%\centering
%\includegraphics[width=130mm]{validation.eps} \caption{
%Demonstration of variance reduction properties of heat kernel
%smoothing. Left: sample variance of thickness for a single subject
%is plotted over the number of iterations. Right:
%
%} \label{fig:validation}
%\end{figure}


\section{Random field theory on cortical manifold}
Here we will describe how to perform multiple comparisons on
$\partial \Omega$ using the random field theory. The random field
theory based approach is widely used for correcting multiple
comparisons in 3D whole brain volume but rarely used on 2D cortical
manifolds \cite{andrade.2001} \cite{chung.cvpr.2003}
\cite{chung.ni.2003} \cite{worsley.2005}. First we combine both the
autistic and the control subjects in a single indexing $j$ and set
up a general linear model (GLM) on cortical thickness $Y_j$ for
subject $j$:

 \bqn K_{\sigma}*Y_j(p) = \lambda_1(p) +
\lambda_2(p) \cdot {\tt{age}}_j +\lambda_3(p) \cdot {\tt{volume}}_j
+ \beta(p) \cdot {\tt{group}}_j + \epsilon_j
\label{eq:ge-age-volume}\eqn is used. Here dummy variable
$\tt{group}$ is $1$ for the autistic subjects and $0$ for the normal
subjects. {\tt volume} is the total gray matter volume for subject
$j$. The total gray matter volume is estimated by computing the
volume bounded by the both outer and inner surfaces
\cite{chung.ni.2003}. The error is modeled as a smooth Gaussian
random field which is viewed as the heat kernel convolution with
Gaussian white noise, i.e. $\epsilon_j = K_{\sigma}*W$. Then we test
the group difference by performing a hypothesis testing:
$$H_0: \beta(p) =0 \mbox{ for all } p \in \partial \Omega$$
v.s.
$$H_1: \beta(p) \neq 0 \mbox{ for some } p
\in \partial \Omega.$$ The test statistic is the ratio of the sum of
the squared residual errors under the null and alternate models.
Under $H_0$, the test statistic is a $F$ random random field with 1
and $n=n_1+n_2-4$ degrees of freedom \cite{worsley.1994}. The null
hypothesis is the intersection of collection of hypothesis $H_0 =
\bigcap_{p \in \partial \Omega} H_0(p)$, where $H_0(p): \beta(p) =
0$ for each fixed $p$. The type I error for the multiple comparisons
is then given by
 \bq\alpha &=&
P\Big(\bigcup_{p \in
\partial \Omega} \{F(p) >
h\}\Big)=1-  P\Big(\bigcap_{p \in \partial \Omega} F(p) \leq h  \}\Big)\\
&=& 1 - P(\sup_{p \in \partial \Omega} F(p) \leq h) = P(\sup_{p \in
\partial \Omega} F(p) > h)\eq
for some $h$.  The resulting $p$-value is usually called the {\em
corrected p-value}. The distribution of $\sup_{p \in \partial
\Omega} F(p)$ is asymptotically given as \bqn P(\sup_{p \in \partial
\Omega} F(p) > h) \approx \sum_{d=0}^2 \phi_d(\partial
\Omega)\rho_d(h) \label{eq:sup} \eqn where $\phi_d$ are the
$d$-dimensional Minkowski functionals of $\partial \Omega$ and
$\rho_d$ are the $d$-dimensional Euler characteristic (EC) density
of $F$-field with $\alpha=1$ and $\beta=n$ degrees of freedom
\cite{worsley.1994}. The Minkowski functionals are $\phi_0 =2,\phi_1
=0, \phi_2 = \mbox{area}(\partial \Omega)/2 = 49,616
 \mbox{mm}^2$, the half area of the template cortex $\partial
\Omega$. The EC density is given by

\bq \rho_0(h) &=&\int_h^\infty {{\Gamma ({{\alpha + \beta}\over
2})}\over {\Gamma ({\alpha \over 2})\Gamma ({\beta \over 2})}}
{\alpha\over \beta}\left({{\alpha x}\over \beta }\right)^{(\alpha-2)
\over 2} \left(1+{{\alpha x}\over \beta
}\right)^{-{(\alpha+\beta)\over 2}} d x,\\
\rho_2(h)&=&\frac{\lambda}{2\pi}{{\Gamma ({{\alpha +\beta-2}\over
2})}\over {\Gamma ({\alpha \over 2})\Gamma ({\beta \over 2})}}
\left({{\alpha h}\over \beta }\right)^{(\alpha-2) \over 2}
\left(1+{{\alpha h}\over \beta }\right)^{-{(\alpha +\beta-2)\over 2}}\\
&\times&\left[(\beta -1){{\alpha h}\over \beta}-(\alpha-1)\right]
\eq where $\lambda$ measures the smoothness of fields $\epsilon$ and
given as $\lambda = 1/(2\sigma^2)$. The resulting corrected
$p$-values maps for $F$ field is shown in Figure
\ref{fig:autism-tstat}. The main use of the corrected $p$-value maps
are the localization and visualization of thickness difference.








\section{Application} $T_1$-weighted MR scans were acquired for 16 autistic and 12 control
subjects on a 3-Tesla GE SIGNA scanner. They are all right-handed
males. 16 autistic subjects were diagnosed with high functioning
autism (HFA). The average age is $17.1 \pm 2.8$ is for the control
subjects and $16.1 \pm 4.5$ for the autistic subjects. The complete
description of the data set, image acquisition parameters, the
subsequent image processing routines, and the interpretation of the
resulting statistical parametric maps is provided in
\cite{chung.ni.2005}. Each image underwent several image
preprocessing steps. Image intensity nonuniformity was corrected
using nonparametric nonuniform intensity normalization method
\cite{Sled}. Then using the automatic image processing pipeline,
%\cite (Zijdenbos {\em et al.}, 1998)
the image was spatially normalized into the Montreal neurological
institute (MNI) stereotaxic space using a global affine
transformation. %\cCollins {\em et al.}, 1994).
Subsequently, 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}. Brain substructures such as the brain stem
and the cerebellum were removed automatically. Triangular meshes for
inner and outer cortical surfaces were obtained by a deformable
surface algorithm \cite{MacDonald.2000}. Such a {\em deformable
surface} approach has the advantage that the surface topology can be
fixed to be spherical and the deformation process can maintain a
non-intersecting surface at all times, obviating the need for
topology correction \cite{dale.1999}. The mesh starts as an
ellipsoid located outside the brain and is shrunk to obtain the
inner cortical surface. Then the inner surface is expanded, with
constraints, to obtain the outer cortical surface. The triangular
meshes are not constrained to lie on voxel boundaries. Instead, the
triangular meshes can cut through a voxel, which serves to reduce
discretization error and partial volume effect. Thickness is
measured using the natural anatomical homology between vertices on
the inner and outer cortical surface meshes, since the outer surface
is obtained by deforming the inner surface.


\begin{figure}[t]
\centering
\includegraphics[scale=0.11]{autism-fstat3.eps}
\caption{Corrected $p$ value maps of $F$-test removing the effect of
age and relative gray matter volume difference projected onto the
average outer (top) and inner surfaces (bottom). It shows relatively
asymmetric thickness difference between two groups.
\label{fig:autism-tstat}}
\end{figure}
Afterwards, thickness measures are smoothed with heat kernel
smoothing with parameters $\sigma=1$ and $k=200$ giving the
effective smoothness of $\sqrt{200}=14.14$ mm. A surface-to-surface
registration to a template surface was performed to facilitate
vertex-by-vertex inter-subject thickness comparison. %\cite{robbins.2003}.
We have formulated it as a registration problem of two functional
data on a unit sphere \cite{thompson.toga.1996}.
%\cite{robbins.2003}.
First a mapping from a cortical surface onto
the sphere is established while recording the mapping. Then cortical
curvatures are mapped onto the sphere. The two curvature functions
on the sphere are aligned by solving a regularization problem that
tries to minimize the discrepancy between two functions while
maximizing the smoothness of the alignment in such a way that the
pattern of gyral
ridges are matched smoothly. %\cite{robbins.2003}.
This alignment is
projected back to the original surface using the recorded mapping.
This regularization mechanism produces a smooth deformation field,
with very little folding.  The deformation field is parameterized
using a triangulated mesh and the algorithm proceeds in a
coarse-to-fine manner, with four levels of mesh resolution. Figure
\ref{fig:spam-cachia-labels} and Figure \ref{fig:cs-spam} illustrate
the effectiveness of this surface registration algorithm by
computing the probability of matching superior temporal sulcal fundi
and central sulci.

After smoothing out thickness measurements, statistical analysis is
performed following the procedures described in the previous
section. The resulting corrected $p$-value map ($<0.1$) for the $F$
statistic is projected onto the template surface for visualization.
Figure \ref{fig:autism-tstat} shows statistically significant
regions of cortical thickness between two groups. After removing the
effect of age and total grey matter volume difference, the
statistically significant regions of thickness decreases are highly
localized at the right inferior orbital prefrontal cortex, the left
superior temporal sulcus and the left occipito-temporal gyrus in
autistic subjects.

\section{Conclusions} This paper has introduced heat kernel
smoothing and its statistical properties for data analysis on the
cortical manifolds. The technique can be used in smooth out data
that is necessary in the random field theory based multiple
comparison correction. We have applied the methodology in detecting
the regions of abnormal cortical thickness in a group of autistic
subjects; however, the approach is not limited to a particular
clinical population. The algorithm is implemented in {\tt MATLAB}
and freely available to download on the web {\tt
http://www.stat.wisc.edu/$\sim$mchung/softwares/hk}. A sample
cortical mesh for a subject and its thickness measures can be also
downloaded from the same website for other researchers.

\section*{Acknowledgement}
Authors wish to thank Kim Dalton and Richard J. Davidson of the
Waisman Laboratory for Brain Imaging and Behavior, University of
Wisconsin-Madison for providing the data illustrated in this study.

\bibliographystyle{plain}
\bibliography{hk.lncs}
\end{document}
