\documentclass[10pt,twocolumn]{article}

%\usepackage{cvpr}
\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}}

\setlength{\textheight}{8.875in} \setlength{\textwidth}{6.875in}
\setlength{\columnsep}{0.4in} \setlength{\topmargin}{0in}
\setlength{\headheight}{0in} \setlength{\headsep}{0in}
\setlength{\parindent}{1pc}
\setlength{\oddsidemargin}{-.1875in}  % Centers text.
\setlength{\evensidemargin}{-.1875in}

% Add the period after section numbers.  Adjust spacing.
\newcommand{\Section}[1]{\vspace{-8pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
\newcommand{\SubSection}[1]{\vspace{-3pt}\subsection{\hskip -1em.~~#1}
        \vspace{-3pt}}


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

\vspace{2cm}
\begin{center}
\Large
{TECHNICAL REPORT NO. 1116}\\
\large January 28, 2005

\vspace{0.5cm} {\Large Model Building in Two-sphere via
Gauss-Weierstrass Kernel Smoothing and Its Application\\
 to Cortical Analysis, Part I.}
% Make title bold and 14 pt font (Latex default is non-bold, 16pt)

\vspace{1cm}

Moo K. Chung\vspace{0.5cm}

Department of Statistics\\
Department of Biostatistics and Medical Informatics\\
Keck Laboratory for Functional Brain Imaging and Behavior\\
%$^4$Department of Psychology, $^5$Department of Psychiatry\\
University of Wisconsin-Madison\\
Madison, WI 53706\\
%$^6$Montreal Neurological Institute, McGill University\\
%Montreal, Quebec, Canada
\end{center}

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

\thispagestyle{empty}
\newpage


\begin{document}

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


\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$Departmnet of Psychology and Psychiatry\\
University of Wisconsin-Madison\\
Madison, WI 53706. USA\\
$^4$Montreal Neurological Institute, McGill University, Canada\\
\tt{mchung@stat.wisc.edu}}

%\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\\
%{\small\url{http://www.author.org/~second}}}
%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
  via the concept of 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
cortical 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 the asymptotic limit 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 vertical
direction on a cortical surface.  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, in turn, 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
regularizing a cost function. Unlike the classical SPHARM, we weight
the measurement such that the closer measurements are weighed more.
The weighted-SPHARM is related to both the classical SPHARM and an
isotropic heat diffusion on a unit sphere as asymptotic limits.

Let us overview some of previous studies 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 \cite{ashburner.1998} do not directly
quantify the amount of tissue growth and atrophy. 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 in Riemannian manifolds. 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 surface in
reducing the data dimension \cite{shen04}. Then they performed the
linear discriminant analysis and support vector machines in surface
classification. 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 noise will be in the high degree spherical
harmonics \cite{gu.2004}. It will be shown that our weighted-SPHARM
more penalizes the high degree spherical harmonics compared to the
classical SPHARM.

Bulow (2004) used spherical harmonics in developing an isotropic
heat diffusion via Fourier transform on a unit sphere as form of
hierarchical surface representation \cite{bulow.2004}. We will show
that our weighted-SPHARM representation is related to Bulow' 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.


Many 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
the mathematical simplicity of manipulating them. In our study, we
will only use real-valued spherical harmonics with different
normalizing constants. This has the effect of decreasing the number
of coefficients to be estimated by half.

Once the differentiable smooth 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 a quadratic
polynomial \cite{Dale, joshi.1997}. The SPHARM-based global
parametrization tend to be computationally expensive compared to the
local quadratic polynomial based parametrization 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}
In this section, we introduce mathematical notations and basic
concepts in SPHARM that are needed in describing the
weighted-SPHARM.
\subsection{Surface 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. 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}).

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 unit sphere. The computation
of the Riemannian metric tensors and 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 \cite{MacDonald:00}. For the visualization purpose, the mean curvature was
   computed and thresholded 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, wahba.1990}. For the completeness of the
exposition, we write the explicit form of the $2l+1$ spherical
harmonics of degree $l$ 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$. 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 $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 and mainly refered as the genearlized Fourier series
expansion. This is the basis of the classical SPHARM representation
of closed 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$. One the other hand, 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 asymptotically related to the weighted-SPHARM.

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
regularization 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 $\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 equation (\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}^d\sum_{m=-l}^l
\beta_{lm}Y_{lm}(p)Y_{l'm'}(p).\label{eq:right}\eqn Integrate the
both sides of the above equation with respect to measure $\mu(p)$.
Then using the othonormal conditition \bq \int_{S^2}
Y_{ij}(p)Y_{lm}(p)\;d\mu(p) =
\left\{\begin{array}{cc} 1, &  i=l, j=m, \\
0,& \mbox{ otherwise},\\
\end{array}
\eq 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}

We briefly explain the numerical implementation issues. We only need
to numerically compute the Fourier coefficients $\langle f,
Y_{lm}\rangle$ in the weighted-SPHARM. The eigenvalues
$\lambda_{lm}$ are given analytically from an analytical 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 and the Riemann sum approximation should
converges to the integral as the mesh resolution increases.

Let $n$ be the total number of nodes in the mesh. Assume triangle
elements $T_{k_1}, \cdots, T_{k_m}$ are adjacent to each other at
node $q_k$ ($1\leq k \leq n$). Then the Fourier coefficients
$\langle f, Y_{lm}\rangle$ is approximated as a limit of the Riemann
sum over triangle elements:
$$\int_{S^2}f(q)Y_{lm}(q) \;d\mu(q) \doteq \frac{1}{3}\sum_{k=1}^n \sum_{i=1}^{k_m}f(q_k)Y_{lm}(q_k)|T_{k_i}|.$$
The weighted-SPHARM is constructed by iteratively adding each term
in Theorem 2.

\subsection{Gauss-Weistrass kernel and heat flow}

For the choice of eigenvalues \bqn \lambda_{lm} = e^{-l(l+1)\sigma},
\label{eq:eigen}\eqn 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. As $\sigma \to \infty$, $\lambda_{lm} \to 1$ and the
weighted-SPHARM representation converges to the classical SPHARM
representation. It is interesting to note that even though the
regularizing cost functions are different in Theorem 1 and Theorem
2, they are related asymptotically.

Another interesting property is observed by noting that
$K_{\sigma}*f$ is the unique solution to the following isotropic
heat diffusion \bqn \frac{\partial g}{\partial t} = \Delta g, \;
g(p,t=0) = f(p) \label{eq:diffusion}\eqn at time $t=\sigma^2/2$
\cite{chung.ni.2005, rosenberg}. This is easily seen from the fact
that $K_{\sigma}$ is the Green's function of equation
(\ref{eq:diffusion}). From this property combined with Theorem 3, we
conclude that the weighted-SPHARM is the finite approximation to the
isotropic heat flow in $S^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. In each case }
\end{table}

We have compared the numerical implementation of the weighted-SPHARM
result 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),$$ where $\delta_{l l'}$ is the Kroneker's delta. 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 heat 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 $\mu$ 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 eigenvalue
(\ref{eq:eigen}) corresponding to the Gauss-Weistrass kernel. The
estimation of the Riemannian metric tensors requires partial
derivatives of $\widehat{\nu}$. 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

%\partial_2 \widehat{ \nu} &=& \sum_{l=0}^{k} \sum_{m=-l}^l \lambda_{lm} \langle
%Z, Y_{lm} \rangle
%\partial_{\varphi} Y_{lm}(\theta,\varphi). \eq
The partial derivatives of spherical harmonics are computed
iteratively. The associated Legendre polynomials in the spherical
harmonic basis are defined as
$$P_l^m(\cos \theta) = \sin ^m \theta \frac{d^m}{dx^m} P_l(x)\Big|_{x=\cos \theta},$$ where $P_l(x)$ is 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} =
\sqrt{g_{11}\cdot g_{22} - g_{12}^2}.$$ 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 of 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 \cite{wahba.1990},
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 $\theta$ in
$$  \frac{1}{2}\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma} =
\sum_{l=0}^{k} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}P_l^0(\cos
\vartheta).$$ The FWHM is then $2\vartheta$. Table 2 in Figure
\ref{fig:fwhm} show the nonlinear relationship between $\sigma$ and
the corresponding FWHM for $k=20$.



\section{Application to autism study} Three Tesla $T_1$-weighted MR scans were acquired for
16 autistic and 12 control right handed 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 group
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. %[blinded]. %\cCollins {\em et al.}, 1994).
Afterwards, an automatic tissue-segmentation algorithm based on a
supervised artificial neural network classifier was used to classify
each voxel as cerebrospinal fluid (CSF), gray matter, or white
matter \cite{Kollakian}. 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. The surface-to-surface registration is completely
automated and has been validated in \cite{Robbins.2003}. 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.

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

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
\cite{chung.ni.2005}. 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.

The weighted-SPHARM is a very flexible functional estimation
technique for scalar and vector data defined 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 the choice of the
Gauss-Weistrass kernel. We can extend this argument further. By
choosing the Green's function of a particular partial differential
equation (PDE) as the kernel in the weighted-SPHARM, we can
construct the finite least squares estimation of the solution of the
PDE without numerically solving PDE using the finite element
technique. This should serve as a spring board for investigating the
wide variety of PDE based data smoothing technique in the local
polynomial 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$ corresponing to
   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.ni.2005}
M.K. Chung, S.~Robbins, Davidson R.J. Alexander~A.L. Dalton, K.M.,
and A.C. Evans.
\newblock Cortical thickness analysis in autism with heat kernel smoothing.
\newblock {\em NeuroImage}, 25:1256-1265, 2005.

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

\bibitem{MacDonald: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} \}.$$
