\documentclass[10pt,twocolumn,letterpaper]{article}


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

% Include other packages here, before hyperref.

% If you comment hyperref and then uncomment it, you should delete
% egpaper.aux before re-running latex.  (Or just hit 'q' on the first latex
% run, let it finish, and you should be clear).
\usepackage[pagebackref=true,breaklinks=true,letterpaper=true,colorlinks,bookmarks=false]{hyperref}


\cvprfinalcopy % *** Uncomment this line for the final submission

%\def\cvprPaperID{****} % *** Enter the CVPR Paper ID here
\def\httilde{\mbox{\tt\raisebox{-.5ex}{\symbol{126}}}}

% Pages are numbered in submission mode, and unnumbered in camera-ready
\ifcvprfinal\pagestyle{empty}\fi

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


\begin{document}

%%%%%%%%% TITLE

\title{Tensor-based Cortical Morphometry via\\
 Weighted Spherical Harmonic Representation}

%\author{First Author\\
%Institution1\\
%Institution1 address\\
%{\tt\small firstauthor@i1.org}
% For a paper whose authors are all at the same institution,
% omit the following lines up until the closing ``}''.
% Additional authors and addresses can be added with ``\and'',
% just like the second author.
% To save space, use either the email address or home page, not both
%\and
%Second Author\\
%Institution2\\
%First line of institution2 address\\
%



\author{Moo K. Chung$^{1,2}$, Steve Robbins$^4$, Kim M. Dalton$^2$\\ Shubing Wang$^1$,
Alan C. Evans$^4$, Richard J. Davidson$^{2,3}$\\
$^1$Department of Statistics, Biostatistics, and Medical Informatics\\
$^2$Waisman Laboratory for Brain Imaging and Behavior\\
$^3$Department of Psychology and Psychiatry\\
University of Wisconsin-Madison\\
%Madison, WI 53706. USA\\
$^4$Montreal Neurological Institute, McGill University, Canada\\
\tt{mchung@stat.wisc.edu}}

\maketitle \thispagestyle{empty}


%Moo K. Chung, Steve Robbins$^2$ and Alan C. Evans
%Department of Statistics, Department of Biostatistics and Medical Informatics\\
%Keck Laboratory for Functional Brain Imaging and Behavior\\
%University of Wisconsin-Madison\\
%Madison, WI 53706. USA\\
%$^2$Montreal Neurological Institute, McGill University, Canada\\
%\tt{mchung@stat.wisc.edu}}
 \maketitle
% \thispagestyle{empty}

%%%%%%%%% ABSTRACT
\begin{abstract}
  We present a new tensor-based morphometric framework that quantifies cortical shape variations
  using the local area element. The local area element is
obtained from the Riemannian metric tensors, which are, in turn,
obtained from the smooth functional parametrization of a triangle
mesh. For the smooth parametrization, we have developed a novel
weighted spherical harmonic (SPHARM) representation. The
weighted-SPHARM differs from the classical SPHARM in a regularizing
cost function. The classical SPHARM is a special case of the
weighted-SPHARM. Further, for a specific choice of weights, the
weighted-SPHARM is shown to be the finite least squares
approximation to the solution of an isotropic heat diffusion on a
unit sphere. The main aims of this paper are to present a
theoretical framework for the weighted-SPHARM, and to show how it
can be used in the tensor-based morphometry. As an illustration, the
methodology has been applied in the problem of detecting abnormal
cortical regions in a clinical population.
\end{abstract}

%%%%%%%%% BODY TEXT
\section{Introduction}
In many previous cortical morphometric studies, cortical thickness
have been mainly used to quantify cortical shape variations in a
population \cite{fischl.1999, MacDonald:00, miller.2000}. The
cortical thickness measures the amount of gray matter in the
vertical direction on the cortex.  We present a new tensor-based
morphometry (TBM) that quantifies the amount of gray matter along
the tangential direction of the cortex by computing the local area
element. The local area element is obtained from the Riemannian
metric tensors, which are computed from the smooth functional
parametrization of a cortical mesh. For this purpose, we present a
novel weighted spherical harmonic (SPHARM) representation that
differs from the classical SPHARM \cite{gerig.2001, shen04} in a
regularizing cost function. Unlike the classical SPHARM, we weigh
measurements in such a way that the closer measurements are weighed
more. The weighted-SPHARM is mathematically related to both the
classical SPHARM and an isotropic heat diffusion on a unit spher.


Let us overview previous literatures that are related to our study.
Gerig et al. (2001) used the mean squared distance (MSD) of SPHARM
coefficients in quantifying ventricle surface shape in a twin study
\cite{gerig.2001}. The distance based metrics widely used in
deformation-based morphometry do not directly
quantify the amount of tissue growth and atrophy \cite{Chung:01a}.
For directly measuring the amount of tissue volume, the Jacobian
determinant of the deformation field is a better metric
\cite{Chung:01a}. Our local area element is the differential
geometric generalization of the Jacobian determinant. So the area
element will be able to quantify the cortical tissue growth/atrophy
directly.


Shen et al. (2004) used the principal component analysis technique
on the SPHARM coefficients of schizophrenic hippocampal surfaces in
reducing the data dimension \cite{shen04}. Then they classified the
hippocampal surfaces using the linear discriminant analysis and a
support vector machine. In a related work, Gu et al. (2004)
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}. It will be shown that the
weighted-SPHARM penalizes high degree spherical harmonics more than
the classical SPHARM does.

Bulow (2004) used the spherical harmonics in developing an isotropic
heat diffusion via the Fourier transform on a unit sphere as a form
of hierarchical surface representation \cite{bulow.2004}. We will
show that the weighted-SPHARM representation is related to the heat
diffusion asymptotically.

%\cite{joshi.1997} used spherical harmonics in expanding a Gaussian
%random field in a unit sphere and used it in the probabilistic
%representation of the hippocampus surface.


Most SPHARM literatures \cite{bulow.2004, gerig.2001, gu.2004,
shen04} use the both real- and imaginary-valued spherical harmonics.
However, the coefficients of imaginary-valued spherical harmonic
basis do not serve any purpose in SPHARM representation other than
providing mathematical simplicity. In this paper, we will use
real-valued spherical harmonics with different normalizing constants
than \cite{bulow.2004, gerig.2001, gu.2004, shen04}.

Once the differentiable parametrization of the cortex is established
by the weighted-SPHARM, we can compute the Riemmanian metric tensors
and local area element. Many previous differential geometric
cortical modeling is based on locally fitting quadratic polynomials
\cite{Dale, joshi.1997}. The SPHARM-based global parametrization
tend to be computationally expensive compared to the local quadratic
polynomial fitting while providing more accuracy and flexibility for
hierarchical representation.

% Surface parameterization of the cortical
%surface has been done previously by \cite{Joshi}. From the surface
%parameterization, Gaussian and mean curvatures of the brain surface
%can be computed and used to characterize its shape \cite{Dale,
%Griffin, Joshi}. In particular, \cite{Joshi} used the quadratic
%surface in estimating the Gaussian and mean curvature of the
%cortical surfaces.

\section{Preliminary}
\subsection{Parametrization} Let $\mathcal{M}$ and $S^2$ be a cortical
surface and a  unit sphere respectively. $\mathcal{M}$ and $S^2$ are
realized as meshes with more than 80,000 triangle elements
\cite{MacDonald:00}. It is natural to assume the cortical surface to
be a smooth 2-dimensional Riemannian manifold parameterized by two
parameters \cite{Davatzikos}. This parametrization is constructed in
the following way. A point $p=(x,y,z) \in \mathcal{M}$ is mapped
onto $u=(u_1,u_2,u_3) \in S^2$ via a deformable surface algorithm
that preserves anatomical homology and the topological connectivity
of meshes (Figure \ref{fig:sphere}) \cite{MacDonald:00}. Let $U$ be
the inverse mapping from $S^2$ to $\mathcal{M}$. Point
$u=(u_1,u_2,u_3) \in S^2$ is parameterized by the spherical
coordinates:
$$(u_1,u_2,u_3)=(\sin \theta \cos \varphi, \sin \theta \sin \varphi,
\cos \theta)$$ with $(\theta, \varphi) \in \mathcal{N}
=[0,\pi]\otimes [0,2\pi)$. This mapping will be denoted as $X$, i.e.
$X: \mathcal{N} \to S^2$. Then we have composite mapping $Z=U \circ
X : \mathcal{N} \to \mathcal{M}$. $Z$ is a 3D vector and it will be
stochastically modeled as \bqn Z(\theta,\varphi) = \nu
(\theta,\varphi)+ \epsilon(\theta,\varphi),\label{eq:linear}\eqn
where $\nu$ is unknown true differentiable parametrization and
$\epsilon$ is a random vector field on $S^2$. The computation of the
Riemannian metric tensors and the local area element require
estimating differentiable function $\nu$.


\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{sphere.eps}
\end{center}
   \caption{Cortical manifold $\mathcal{M}$ (left) is mapped onto unit sphere $S^2$ (right)
   via a deformable surface algorithm that preserves anatomical
   homology and surface topology \cite{MacDonald:00}. For the visualization purpose, the mean curvature was
   computed and segmented to better represent sulci and
   gyri.} \label{fig:sphere}
\end{figure}



\subsection{Spherical harmonic representation} The basis functions on the unit sphere are given as
the eigenfunctions satisfying $\Delta f + \lambda f =0,$ where
$\Delta$ is the spherical Laplacian:
$$\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, denoted as  $Y_{lm} ( -l \leq m
\leq l)$, corresponding to the same eigenvalue $\lambda = l(l+1)$.
$Y_{lm}$ is called the {\em spherical harmonic} of degree $l$ and
order $m$ \cite{courant}. The explicit form of the
$2l+1$ spherical harmonics of degree $l$ are given 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 SPHARM literatures \cite{bulow.2004, gerig.2001,
gu.2004, shen04} that used the complex-valued spherical harmonics,
we use only real-valued spherical harmonics with different
normalizing constants since they are more convenient for a
real-valued stochastic model (\ref{eq:linear}).

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 Lebesgue measure
$d\mu(\theta,\varphi)=\sin \theta d\theta d\varphi$. With respect to this 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. 

Consider
subspace $$\mathcal{H}_k = \{ \sum_{l=0}^k \sum_{m=-l}^l\beta_i
Y_{lm}: \beta_i \in \mathbb{R}\} \subset L^2(S^2),$ which is spanned
by up to the $k$-th degree spherical harmonics. Then the least
squares estimation, denoted as $\widehat f$, of $f \in L^2(S^2)$ in
the subspace $\mathcal{H}_k$ is given by
$$\hat f(p) = \sum_{l=0}^{k} \sum_{m=-l }^{l}\langle f, Y_{lm}
\rangle Y_{lm}(p).$$ This can be stated as the following theorem.\\
{\bf Theorem 1.}
$$\sum_{l=0}^{k} \sum_{m=-l}^l \langle f, Y_{lm} \rangle Y_{lm} = \arg\min_{h \in \mathcal{H}_k}
\int_{S^2}\big[f(q)- h(q)\big]^2\;d\mu(q).$$ This theorem is a well
known result in Fourier analysis and mainly refered as the
genearlized Fourier series expansion. This is the basis of the
classical-SPHARM representation for anatomical boundaries
\cite{gerig.2001, gu.2004, shen04}.


\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{schmatic5.eps}
\end{center}
   \caption{The schematic comparison of the classical-SPHARM and
   weighted-SPHARM. The classical approach estimates functional data $f$ by minimizing the
   integrated squared distance between $f$ and smooth function
   $h$. This distance is indicated by an arrow. The weighted-SPHARM estimates $f$ locally at each fixed $p$ by
   minimizing the integrated weighted squared distance
   between $f$ and $h(p)$. The weighted-SPHARM can be
   viewed as a local kernel regression \cite{fan}.} \label{fig:schmatic}
\end{figure}

\section{Weighted-SPHARM}
\subsection{Basic theory.} 
The classical-SPHARM is only one possible
representation of functional data measured on the unit sphere. We
will present a more general representation technique in the
framework of a local kernel regression \cite{fan}. We will call this
technique as the {\em weighted-SPHARM} since the coefficients of
SPHARM are additionally weighted by the eigenvalues of a kernel. It
will be shown that the classical-SPHARM is a special case of the
more general weighted-SPHARM representation.

First, we start with the spectral representation of positive
definite kernel in $S^2$. Any positive definite kernel $K(p,q)$ in
$S^2$ can be represented as \bqn K(p,q) = \sum_{l=0}^{\infty}
\sum_{m=-l}^l \lambda_{lm} Y_{lm}(p)Y_{lm}(q),\label{eq:kernel}\eqn
where eigenvalues $\lambda_{00} \geq \lambda_{1m_1} \geq
\lambda_{2m_2} \geq \cdots \geq 0$ satisfy
$$\int_{S^2} K(p,q)Y_{lm}(q) \;d\mu(q) = \lambda_{lm} Y_{lm}(p).$$
This is the special case of the Mercer's theorem \cite{courant}.
Without loss of generality, we assume the kernel is normalized in
such a way that \bqn \int_{S^2} K_(p,q) \;d\mu(q)=1.
\label{eq:dist}\eqn  At each fixed point $p$, smooth representation
$h$ of functional data $f$ is searched in the subspace
$\mathcal{H}_k$ that minimizes the integral of the weighted squared
distance between $f$ and $h$. This is formulated as the following minimization problem: \bqn \min_{h \in \mathcal{H}_k}\int_{S^2}
K(p,q)\big[f(q) - h(p)\big]^2 \; d\mu(q).\label{eq:min}\eqn See
Figure \ref{fig:schmatic} for the schematic comparison of the
classical SPHARM and the weighed-SPHARM. The minimizer of
(\ref{eq:min}) is given by the following
theorem.\\
{\bf Theorem 2.} \bq &&\sum_{l=0}^k \sum_{m=-l}^{l} \lambda_{lm}\langle f, Y_{lm}\rangle Y_{lm}\\
&=&\arg\min_{h \in \mathcal{H}_k}\int_{S^2} K(p,q)\big[f(q) -
h(p)\big]^2 \; d\mu(q).\label{theorem:finite.polynomial}\eq {\em
Proof.} Let $$h(p)= \sum_{l=0}^k\sum_{m=-l}^l \beta_{lm}Y_{lm}(p)
\in \mathcal{H}_k.$$ The integral can be written as
\bq & &I(\beta_{00}, \beta_{1-1},\beta_{10},\beta_{11}, \cdots, \beta_{kk})\\
&=&\int_{S^2} K(p,q)\Big[f(q) - \sum_{l=0}^k\sum_{m=-l}^l
\beta_{lm}Y_{lm}(p)\Big]^2\;d\mu(q).\eq Since the functional $I$ is
quadratic in coefficients $\beta_{lm}$, the minimum exists and it is obtained
when \bqn \frac{\partial I}{\partial \beta_{l'm'}}=0 \mbox{ for all
} l' \mbox{ and } m'.\label{eq:condition}\eqn Then solving equation
(\ref{eq:condition}) with condition (\ref{eq:dist}), we have
\bqn &&Y_{l'm'}(p)\int_{S^2} K(p,q)f(q)\;d\mu(q) \label{eq:left}\\
&=& \sum_{l=0}^k\sum_{m=-l}^l
\beta_{lm}Y_{lm}(p)Y_{l'm'}(p).\label{eq:right}\eqn Integrate the
both sides of the equation with respect to measure $\mu(p)$.
we obtain \bq \beta_{l'm'} &=& \int_{S^2}
f(q)\;d\mu(q)\int_{S^2} Y_{l'm'}(p) K(p,q) \;d\mu(p)\\
&=& \sum_{l=0}^{\infty} \sum_{m=-l}^l \lambda_{lm}\int_{S^2}
f(q)Y_{lm}(q)\;d\mu(q)\\
&& \times \int_{S^2}Y_{lm}(p)Y_{l'm'}(p)\;d\mu(p)\\
&=& \lambda_{lm}\int_{S^2} f(q)Y_{l'm'}(q)\;d\mu(q).\eq This proves
the statement.


Now we show what happens as the dimension of $\mathcal{H}_k$
increases. Define kernel smoothing as the integral convolution \bqn
K*f(p)=\int_{S^2} f(q)K(p,q)\;d\mu(q).\label{eq:conv}\eqn Then it
can be shown that the weighted-SPHARM converges to kernel smoothing
(\ref{eq:conv}) as the dimension of subspace $\mathcal{H}_k$
increases. This can be stated differently as\\
{\bf Theorem 3.} \bq K*f(p)=  \arg\min_{h \in L^2(S^2)} \int_{S^2}
K(p,q)\big[f(q) - h(p)\big]^2 \; d\mu(q).\eq {\em Proof.} The
weighted-SPHARM representation can be rearranged as \bq & &
\sum_{l=0}^k \sum_{m=-l}^{l}
\lambda_{lm}\langle f, Y_{lm}\rangle Y_{lm}(p)\\
& =& \int_{S^2} f(q) \sum_{l=0}^k  \sum_{m=-l}^l
\lambda_{lm}Y_{lm}(p)Y_{lm}(q)\;d\mu(q)\\
&\to&  \int_{S^2} f(q) K(p,q)\;d\mu(q) \;\mbox{ as } k \to \infty\eq
The last line is from equation (\ref{eq:kernel}). On the other hand,
from the completeness of Hilbert space $L^2(S^2)$, \bq &&\lim_{k \to
\infty} \arg\min_{h \in \mathcal{H}_k}\int_{S^2}
K(p,q)\big[f(q) - h(p)\big]^2 \; d\mu(q)\\
& =& \arg\min_{h \in L^2(S^2)}\int_{S^2} K(p,q)\big[f(q) -
h(p)\big]^2 \; d\mu(q).\eq This proves the statement. Theorem 3
connects the weighted-SPHARM to kernel smoothing as the asymptotic
limit.

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{xyz2.eps}
\end{center}
   \caption{Top: the original inverse mapping $U$ is displayed in
   $S^2$. It shows the coordinate functions projected onto $S^2$.
   Bottom: the weighed-SPHARM representation of the coordinate
   functions. The color scale for coordinates is
   thresholded at $\pm 45 \mbox{ mm}$ to better show the smoothing pattern of the
   weighted-SPHARM representation.} \label{fig:sulcal}
\end{figure}


For the choice of eigenvalues $ \lambda_{lm} = e^{-l(l+1)\sigma}$, the corresponding kernel is called the {\em
Gauss-Weistrass kernel} and it will be denoted as \bqn
K_{\sigma}(p,q) = \sum_{l=0}^{\infty} \sum_{m=-l}^l
e^{-l(l+1)\sigma}Y_{lm}(p)Y_{lm}(q).\label{eq:GWK}\eqn The subscript
$\sigma$ is introduced to indicate the dependence of the additional
parameter. When $\sigma =0$, $\lambda_{lm} = 1$ and the
weighted-SPHARM becomes the classical-SPHARM. It is interesting to note that even though the
regularizing cost functions are different in Theorem 1 and Theorem
2, they are related. Another interesting property is that
$K_{\sigma}*f$ is the unique solution to the following isotropic diffusion equation 
\bqn \frac{\partial g}{\partial \sigma} = \Delta g, \;
g(p,\sigma=0) = f(p) \label{eq:diffusion}\eqn
\cite{rosenberg}. From this property combined with Theorem 3, we
conclude that the weighted-SPHARM is the finite approximation of the
isotropic diffusion in $S^2$.

\subsection{Numerical Implementation}
We only need to numerically estimate the Fourier coefficients $\langle f,
Y_{lm}\rangle$ in the weighted-SPHARM. The eigenvalues
$\lambda_{lm}$ are given analytically from a given kernel. The
computation for the Fourier coefficients are based on the direct
numerical integration over high resolution triangle meshes with more
than 80,000 triangles and the average inter-vertex distance of
0.0189 mm. The accuracy of the weighed-SPAHRM is only restricted to
the mesh resolution.  Then the Fourier coefficients
$\langle f, Y_{lm}\rangle$ is approximated as the Riemann
sum over triangle elements. The Riemann sum based approximation should
converge to the integral as the mesh resolution increases. The weighted-SPHARM is constructed by iteratively adding each term in Theorem 2.

\begin{table}[b]
\begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
$l'$ & $m'$ & integral & difference\\
\hline
20 & 4 &  1.0001  &  $9.7029 \cdot 10^{-5}$  \\
20   & 10 & 0.9999 &  $1.6212 \cdot 10^{-4}$ \\
20 & 20 &  0.9999 & $-1.1174 \cdot 10^{-4}$ \\
\hline
\end{tabular}
\end{center}
\caption{Numerical accuracy of the weighted-SPHARM with
$\sigma=0.01$ for degree 20, and order 4, 10 and 20. The third
column checks if $\langle Y_{l'm'},Y_{l'm'} \rangle = 1$. The last
column shows the average difference between the weighted-SPHARM and
the expected heat diffusion.}
\end{table}

We have compared the numerical results of the weighted-SPHARM against the analytical solution of (\ref{eq:diffusion}). Let $f=e^{l'(l'+1)}Y_{l'm'}$ be an analytic test function. Then
$K_{\sigma}*f$ can be written as
$$e^{l'(l'+1)}\sum_{l=0}^{\infty} \sum_{m=-l}^{l} e^{-l(l+1)\sigma}
Y_{lm}(p)\int_{S^2} Y_{lm}(q) Y_{l'm'}(q) \; d\mu(q)$$
$$=e^{l'(l'+1)}\sum_{l=0}^{\infty} \sum_{m=-l}^{l}\\
e^{-l(l+1)\sigma} Y_{lm}(p) \delta_{l l'} \delta_{m m'}
=Y_{l'm'}(p).$$ The
table 1 shows the comparative result for $l'=20$ and selective $m'$
with $\sigma=0.01$ and degree $k=20$. The third column shows the
numerical computation of integral $\int_{S^2} Y_{l'm'}^2(p)\;dp = 1$
showing the accuracy up to 3 decimal places. This shows our Riemann
sum approximation provides sufficiently good accuracy, which depends
on the mesh resolution. The fourth column shows the average
difference between the weighed-SPHARM and the isotropic diffusion.

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{metric.eps}
\end{center}
   \caption{Riemannian metric tensor estimation. The metric tensors $g_{ij}$ are
   estimated by differentiating the weighted-SPHARM representation. Afterwards the local area element
   $\sqrt {\det g}$ is computed. The local area
   element measures the amount of are expansion and shrinking with respect to $S^2$.} \label{fig:sulcal}
\end{figure}

\subsection{Riemannian metric tensor estimation}

The weighted-SPHARM estimation $\widehat{\nu}$ of the unknown true
parametrization $\nu$ in equation (\ref{eq:linear}) is given by
$$\widehat{\nu} = \sum_{l=0}^k \sum_{m=-l}^l \lambda_{lm} \langle Z,
Y_{lm} \rangle Y_{lm}.$$ For this study, we used eigenvalues corresponding to the Gauss-Weistrass kernel. Denoting partial differential
operators as $\partial_1 =
\partial_{\theta}$ and $\partial_2 =
\partial_{\varphi}$, we have derivative estimations \bq
\partial_i \widehat{ \nu} &=& \sum_{l=0}^{k} \sum_{m=-l}^l \lambda_{lm}
\langle Z, Y_{lm} \rangle \partial_{i} Y_{lm}(\theta,\varphi).\eq
The partial derivatives of spherical harmonics are iteratively computed. The associated Legendre polynomials in the spherical
harmonic basis are given by
$$P_l^m(\cos \theta) = \sin ^m \theta \frac{d^m}{dx^m} P_l(x)\Big|_{x=\cos \theta},$$ where $P_l(x)$ are the
Legendre polynomials defined in $(-1,1)$ with $P_0(x)=1$ and
$P_1(x)=x$. Then for $0 \leq m \leq l-1$, \bq \partial_{\theta} d
P_l^m(\cos \theta) &=& m \sin ^{m-1} \theta \cos \theta
\frac{d^m}{dx^m} P_l(x)\Big|_{x=\cos \theta}\\ &&- \sin ^{m+1}
\theta \frac{d^{m+1}}{dx^{m+1}}
P_l(x)\Big|_{x=\cos \theta}\\
&=& m\cot \theta P_l^m(\cos \theta) - P_l^{m+1}(\cos \theta).\eq For
$m=l$, since $P_l$ is the $l$-th order polynomial, the second term
vanishes. A similar recursive relationship for an alternate
definition for the associated Legendre polynomial is given in
\cite{khan.2003}. Based on this iterative relation, we can compute
the partial derivatives \bq \partial_{\theta} Y_{lm} =
\left\{\begin{array}{cc} c_{lm}\partial_{\theta}
P^{|m|}_l(\cos\theta)\cos
(|m|\varphi), &-l \leq m\leq -1, \\
 \frac{c_{lm}}{\sqrt{2}}\partial_{\theta}P_l^0(\cos \thata),& m=0,\\
c_{lm} \partial_{\theta}P^{|m|
}_l(\cos\theta)\sin (|m|\varphi),& 1 \leq m\leq l\\
\end{array}\eq
and
 \bq \partial_{\varphi} Y_{lm} = \left\{\begin{array}{cc}
|m|c_{lm}P^{|m|}_l(\cos\theta)\cos
(|m|\varphi), &-l \leq m\leq -1, \\
 0,& m=0,\\
-|m|c_{lm} P^{|m|
}_l(\cos\theta)\sin (|m|\varphi),& 1 \leq m\leq l.\\
\end{array}\eq
Then the Riemannian metric tensors are estimated as $g=(g_{ij})
=\langle \partial_i \widehat{\nu}, \partial_j \widehat{\nu} \rangle$
and the area element $G(\theta,\varphi) = \sqrt{\det g} .$ The area element measures
the transformed area in $\mathcal{M}$ of the unit square of the
parameterized space $\mathcal{N}$ via mapping $\nu$.

\section{Statistical inference in $S^2$}

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{FWHM3.eps}
\end{center}
   \caption{Plot of $\sigma$ (horizontal) vs. FWHM (vertical) showing
   the nonlinear functional relationship. Table 2. shows FWHM for different
   $\sigma$.} \label{fig:fwhm}
\end{figure}

For the $i$-th subject ($1 \leq i \leq m$), we denote the cortical
manifold as $\mathcal{M}_i$ and its area element as
$G_i(\theta,\varphi)$. The area element is influenced by the global brain size. If we
enlarge the cortical coordinates by the factor of $r$, the area
element changes by the factor of $r^2$. So it is necessary to
normalize $G_i$ such that it is invariant under scaling.  The affine
scale invariant area element is given by
$$\widetilde{G_i} (\theta,\varphi)=\frac{4\pi
G(\theta,\varphi)}{\mu(\mathcal{M}_i)},$$ where $\mu(\mathcal{M}_i)$
is the total cortical area. If we enlarge the the cortical
coordinates by the factor of $r$, $\mu(\mathcal{M}_i)$ increases by
the factor of $r^2$ making $\widetilde{G_i}$ invariant under affine
scaling. The constant $4\pi$ is multiplied so that the normalization
is with respect to the total surface area of $S^2$. Then we have the
following general linear model (GLM):
$$\widetilde{G_i}(\theta,\varphi) = \alpha_0 + \alpha_1
\cdot{\tt age}_i + \alpha_2(\theta,\varphi) \cdot{\tt group}_i +
\epsilon(\theta,\varphi),$$ where $\epsilon$ is a mean zero Gaussian
random field. ${\tt age}_i$ and ${\tt group}_i$ are the age and a
categorical dummy variable ($0$ for autism and $1$ for control)
respectively for subject $i$. Then we test if there is any group
difference in the local area element measure by testing \bq && H_0:
\alpha_2(\theta,\varphi) =0 \mbox{ for all } \theta \mbox{ and }
\varphi.\\
&\mbox{ vs. }& H_1: \alpha_2(\theta,\varphi) \neq 0 \mbox{ for some
} \theta \mbox{ and } \varphi.\eq At each point $(\theta,\varphi)$,
a $F$-statistic with 1 and $n-3$ degrees of freedom, denoted as
$F(\theta,\varphi)$ is used as a test statistic. The $F$-statistic
is constructed as a ratio of the residual sum of error of model fit
of $H_0$ and $H_1$. Since we need to perform the test at every
$(\theta,\varphi)$, this becomes a multiple comparison problem. We
used the random field theory \cite{taylor.2006, worsley:96a} based
thresholding to determine the statistical significance.

The probability of obtaining false positives ($\alpha$-level) for
the one sided alternate hypothesis in given by
$$P( \sup F(\theta,\varphi) \geq F_{\alpha}) = \sum_{i=0}^2
\frac{\mathcal{L}_i(S^2)}{\mbox{FWHM}^i}\rho_i(y),$$
 where $\mathcal{L}_i$ is the $i$-th Lipschitz-Killing curvature or
 Minkowski functional \cite{taylor.2006}, and $\rho_i$ is the $i$-dimensional
 EC-density \cite{worsley:96a}. FWHM denotes the full width at the half
 maximum of smoothing kernel $K_{\sigma}$ used in the weighted sphearical harmonic
 representation. For the unit sphere, the Lipschitz-Killing curvatures
 are
$$\mathcal{L}_0(S^2)=2, \;
\mathcal{L}_1(S^2)=0, \; \mbox{ and } \mathcal{L}_2(S^2)= 2\pi.$$
The EC-densities are $$\rho_0(y) &=& \int_y^{\infty}
\frac{\Gamma(\frac{m}{2})}{((m-1)\pi)^{1/2}\Gamma(\frac{m-1}{2})}
\Big(1+\frac{y^2}{m-1}\Big)^{-m/2}\; dy$$
$$
\rho_2(y) &=&\frac{4 \ln 2}{(2\pi)^{3/2}}
\frac{\Gamma(\frac{m}{2})y\Big(1+\frac{y^2}{m-1}\Big)^{-(m-2)/2}}{(\frac{m-1}{2})^{1/2}
\Gamma(\frac{m-1}{2})}.$$


\subsection{Computing FWHM}
The computation for the FWHM of the Gauss-Weistrass kernel in $S^2$
is not trivial due to the fact there is no known close form
expression for the FWHM as a function of $\sigma$. So the FWHM is
computed numerically.


The Gauss-Weistrass kernel can be simplified from equation
(\ref{eq:GWK}), via the harmonic addition theorem, as \bqn K_{\sigma}(p,q) = \sum_{l=0}^{\infty}
\frac{2l+1}{4\pi}e^{-l(l+1)\sigma}P_l^0(\cos
\vartheta),\label{eq:theorem}\eqn where $\vartheta$ is the angle
between $p$ and $q$. Using the vector inner product $\cdot$, the
angle can be written as $\cos \vartheta = p \cdot q$. The maximum of
the Gauss-Weistrass kernel is obtained when $\vartheta = 0$ and it
is given by
$$\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}.$$
Now we fix $\varphi=0$ and let $p$ be the north pole, i.e.
$p=(0,0,1)$. By varying $q= (\sin \vartheta, 0, \cos \vartheta)$ for
$0 \leq \vartheta = \cos^{-1} (p \cdot q) \leq \pi$, we have $Y_{lm}
=0$ if $m \neq 0$. Note $P_l^0(1) = 1$ for all $l$. Then we solve
numerically for $\vartheta$ 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
\vartheta).$$ The FWHM is then $2\vartheta$. Table 2 in Figure
\ref{fig:fwhm} shows the nonlinear relationship between $\sigma$ and
the corresponding FWHM for $k=20$.
 \begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{sulcal.eps}
\end{center}
   \caption{Demonstration of cortical surface normalization in $S^2$ showing the nonlinear alignment
   of central and superior temporal sulci for 149 subjects. Left: before normalization. Right: after normalization. The probability of
   matching increases after normalization.} \label{fig:sulcal}
\end{figure}
 
\section{Application to autism study} Three Tesla $T_1$-weighted MR scans were acquired for
16 autistic and 12 control  males. 16 autistic subjects were diagnosed with high functioning autism. The average ages are
$17.1 \pm 2.8$ and $16.1 \pm 4.5$ for control and autistic groups
respectively. Image intensity nonuniformity was corrected using
nonparametric nonuniform intensity normalization method and then the
image was spatially normalized into the Montreal neurological
institute (MNI) stereotaxic space using a global affine transformation.  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. Triangular meshes for outer cortical surfaces were
obtained by a deformable surface algorithm \cite{MacDonald:00}. The mesh
starts as an ellipsoid located outside the brain and is shrunk to
match the cortical boundary. By performing an affine transform on
this ellipsoid, we obtain $S^2$ mesh, which is used in the
weighted-SPHARM.

The segmented cortical meshes are normalized via a nonlinear
surface-to-surface registration  \cite{Robbins.2003}.  Cortical curvatures of
two surfaces are mapped onto the sphere and they are aligned by
solving a regularization problem that tries to minimize the
discrepancy between two curvatures while maximizing the smoothness
of the alignment in such a way that the pattern of gyral ridges are
matched smoothly. 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:sulcal} demonstrates the effectiveness of this surface
registration algorithm by showing the increased matching probability
of superior temporal and central sulci for 149 subjects.

Afterwards, the weighted-SPHARM representation is used as an
estimate for differentiable smooth parametrization of the meshes. We
used parameters $k=20$ and $\sigma=0.001$. The corresponding FWHM is
0.2262 mm. The Riemannian metric tensors and the area element are
computed simultaneously. Based on the general linear model
framework, the statistical parametric map (F-statistic) is computed
and projected on both $S^2$ and the average cortical surface (Figure
\ref{fig:spm2}). The average cortical surface is constructed by
averaging the anatomically corresponding vertices in the meshes. 
It serves as an anatomical landmark for showing where
signals are detected. Then we performed the random field theory
based multiple comparison correction on the F-statistic at $\alpha
=0.05$ level but we did not detect any statistically significant
cortical regions of local area difference although we observed
maximum F-value of 9.3 (uncorrected p-value of 0.0054) at the left
temporal lobe.

\section{Conclusions and discussions}
In this paper, we presented a theoretical framework for the
weighted-SPHARM and its application in TBM. The weighted-SPHARM is
used as a differentiable parametrization of the cortex. This enable
us to compute the Riemannian metric tensors and local area element.
The local area element is used in determining the statistical
significance of the abnormal cortical tissue expansion/shrinking for
a clinical population. Unfortunately, we did not detect any statistically significant signal possibly due to low degree ($k=20$) used. So it is necessary to increase the computational efficiency using a fast Fourier coefficient estimation technique such as the iterated residual fitting (IRF) \cite{chung.2006.1118}.

The weighted-SPHARM is a very flexible functional estimation
technique for scalar and vector data obtained in $S^2$. It was shown
that the solution to the isotropic heat diffusion in $S^2$ is the
asymptotic limit of the weighted-SPHARM for particular weights. We can extend this argument further.  
The solution of  any linear self-adjoint partial differential equation (PDE) can be expressed as using the weighted-SPAHRM representation so we can avoid using numerical technique such as the finite element method in solving PDEs in $S^2$. This should serve as a spring board for investigating the
wide variety of PDE based data smoothing technique in the kernel regression framework \cite{fan}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\linewidth]{spm2.eps}
\end{center}
   \caption{The final statistical analysis result. The $F$ statistic
   is computed at every vertices of the triangle mesh in $S^2$ (left) and mapped
   onto the average cortex (right). The average cortex serves as an anatomical
   reference. The random field theory based multiple comparison at
   $\alpha=0.05$ did not show any significant result although the
   maximum signal is detected at the left temporal lobe ($\sup F=9.3$ corresponds to
   the uncorrected $p$-value of 0.0054).} \label{fig:spm2}

\end{figure}




 \small
\bibliographystyle{ieee}
%\bibliography{egbib}
\begin{thebibliography}{10}

%\bibitem{ashburner.1998} J. Ashburner, C. Hutton, R.S.J. Frackowiak,
%I. Johnsrude, C. Price, and K.J. Friston. Identifying global
%anatomical differences: Deformation-based morphometry. {\em Human
%Brain Mapping}. 6:348-357. 1998.

\bibitem{bulow.2004}
T.~Bulow,
\newblock  Spherical diffusion for 3d surface smoothing,
\newblock {\em IEEE Transactions on Pattern Analysis and Machine Intelligence},
  vol. 26, pp. 1650--1654, 2004.

\bibitem{Chung:01a}
 M.K. Chung, K.J. Worsley, T. Paus, C. Cherif, J.N.  Giedd, J.L.
Rapoport and A.C. Evans.
\newblock A Unified Statistical Approach to Deformation-Based Morphometry,
\newblock {\em NeuroImage} 14:595-606, 2001.

\bibitem{chung.2006.1118}
M.K. Chung, L. Shen, K.M. Dalton, D.J. Kelley, S.M. Robbins, A.C. Evans, R.J. Davidson,
\newblock Weighted spherical harmonic representation and its application
to cortical analysis. \newblock TR 1118. Dept. of Stat., Univ. of Wisconsin-Madison, 2006.\\
%{\tt http://www.stat.wisc.edu/$\sim$mchung/papers/TR1118.2006.pdf}, 2006

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


\bibitem{Dale}
A.M. Dale and B. Fischl.
\newblock Cortical Surface-based Analysis I. Segmentation and Surface
  Reconstruction.
\newblock {\em NeuroImage}, 9:179--194. 1999

\bibitem{Davatzikos}
C. Davatzikos and R.N. Bryan.
\newblock Using a Deformable Surface Model to Obtain a Shape Representation of
  the Cortex.
\newblock In {\em Proceedings of the IEEE International Conference on Computer
  Vision}, pages 2122--2127. 1995



%\bibitem{geric.2002}
%M.S. G.~Gerig and G.~Szekely,
%\newblock ``Statistical shape models for segmentation and structural
%  analysis,''
%\newblock in {\em Proceedings of IEEE International Symposium on Biomedical
%  Imaging (ISBI)}, 2004.




\bibitem{joshi.1997}
S. Joshi, U. Grenander, and M.I. Miller, On the geometry and shape
of brain sub-manifolds, {\em International Journal of Pattern
Recognition and Artificial
    Intelligence: Special Issue on Processing of MR Images of the
    Human}. 11:1317--1343. 1997.


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

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

\bibitem{gu.2004}
X.~Gu, Y.L.~Wang, T.F.~Chan, P.M.~Thompson,
S.T.~Yau. Genus Zero Surface Conformal Mapping and its Application
to Brain Surface Mapping, {\em IEEE Transactions on Medical
Imaging}, 23:1--10. 2004

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

\bibitem{fischl.1999}
B.~Fischl, and M.I.~Sereno, and R.~Tootell, and A.M.~Dale,
High-resolution intersubject averaging and a coordinate system for
the cortical surface {\em Hum. Brain Mapping}, vol. 8, pp. 272--284,
1999.

\bibitem{khan.2003}
T.~Khan, and A.~Thomas. On derivation of the radiative transfer
equation and its spherical harmonics approximation for scattering
media with spatially varying refractive indices. Technical Report.
2003

\bibitem{MacDonald:00}
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{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{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{rosenberg}
S.~Rosenberg,
\newblock {\em The Laplacian on a Riemannian Manifold},
\newblock Cambridge University Press, 1997.

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

%\bibitem{styner.2004}
%M. Styner, G. Gerig: Correction scheme for multiple correlated
%statistical tests in local shape analysis. Proc. SPIE Vol. 5370,
%Medical Imaging 2004: Image Processing; J. Michael Fitzpatrick,
%Milan Sonka; Eds, p. 233-240,2004.

\bibitem{taylor.2006}
J. E. Taylor. A Gaussian Kinematic Formula. {\em Annals of
Probability.} 2006. in press.

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

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

\bibitem{worsley:96a}
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.

%\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage
%\rule{0pt}{1pt}\newpage

\end{thebibliography}
\end{document}


The relationship between the SPRM and our new weighted SPRM is that
our approach reduces the least squares estimation error. Note that
\bq && \int_{S^2} K_{\sigma}(p,q)\big[f(q) - h(p)\big]^2 \; d\mu(q)\\
&\leq &\int_{S^2} K_{\sigma}(p,q)\; d\mu(q)\int_{S^2} \big[ f(q) -
h(p)\big]^2 \;d\mu(q)\\
&=& \int_{S^2} \big[f(q) - h(p)\big]^2 \;d\mu(q).\eq
 Hence, \bq &&\min_{h
\in \mathcal{H}_k}\int_{S^2} K_{\sigma}(p,q)\big[f(q)
- h(p)\big]^2 \; d\mu(q) \\
&\leq&  \min_{h \in \mathcal{H}_k}\int_{S^2} \big[f(q) - h(p)\big]^2
\;d\mu(q).\\
&\leq& \min_{h \in \mathcal{H}_k}\int_{S^2} \big[f(q) - h(q)\big]^2
\;d\mu(q). \eq This result shows that the residual error of weighted
SPRM is always smaller than that of SPRM. An additional properties
of the weighted SPRM can be obtained if we specify the analytical
form of kernel $K_{\sigma}$.


Probability of Type I error increases with the number of tests on
the unit sphere. For $n$ tests on the unit sphere, there are $n$
corresponding p-values $p_1,p_2,\cdots,p_n$. The Benjamini-Hochberg
procedure (1995) orders the p-values in in increasing order
$p_{(1)},p_{(2)},\cdots,p_{(n)}$. Then the Benjamini-Hochberg
threshold is given by
$$\max_{1  \leq i \leq n}\{ P_{(i)}: P_{(i)} \leq \alpha \frac{i}{n} \}.$$
