% Run following
%Prostscript from a DVI file is to use the dvips command with the
%-Ppdf flag. It is recommended that LaTex papers containing images
%should also use the -G0 flag with the dvips command to produce the
%best results.

\documentclass{article}
\usepackage{spconf,amsmath,epsfig}
\usepackage{amsfonts,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}}



\title{Heat kernel smoothing on unit sphere}

%\name{ Moo K. Chung\thanks{Thanks to XYZ agency for funding.}, Jonathan Taylor}
%\address{School\\
%        Department\\
%        Address}

\name{Moo K. Chung}
\address{Department of Statistics, Biostatistics, and Medical Informatics\\
Waisman laboratory for brain imaging and behavior\\
    University of Wisconsin\\
         Madison, WI 53706\\
         $\tt{mchung@stat.wisc.edu}$}


%)\thanks{Thanks to XYZ agency for funding.}}
%\address{Author Affiliation(s)}


%\twoauthors {Moo K. Chung}
%         {
 %{Jonathan Taylor}
 %        {Stanford University\\
 %        Department of Statistics\\
 %        Stanford, CA 94305\\
 %        $\tt{jonathan.taylor@stanford.edu}$}

\begin{document}

\maketitle
%

\begin{abstract}

In brain imaging, cortical data such as the cortical thickness,
cortical surface curvatures and surface coordinates have been mapped
to a unit sphere for the purpose of visualization, surface
registration and statistical analysis. Since the
 unit sphere provides a readily available parametrization and
basis functions, cortical data can be easily quantified with respect
to the spherical parametrization. For the cortical data on the unit
sphere, it is necessary to smooth them to increase the
signal-to-noise ratio and the smoothness for the subsequent
statistical analysis. We present a mathematical framework for
smoothing data on a unit sphere using the heat kernel. The heat
kernel is analytically constructed using the spherical harmonics and
$O(n)$ heat kernel smoothing algorithm is presented.

\end{abstract}

\section{Introduction}

Due to the highly convoluted nature of the human brain cortex, some
sort of normalization onto a parameterized standard surface is
necessary for quantifying cortical data. One such normalization
technique is the cortical flattening that maps the convoluted
cortical manifold into a unit sphere with the minimum geometric
distortion \cite{dale.1999, hurdal.1999}. Alternately, one may
flatten the cortex using the inverse process of a deformable surface
algorithm, which deforms a sphere to match the cortex
\cite{macdonald.2000}.

Once the cortical data is mapped onto the unit sphere, spherical
harmonics have been used in modeling the data. The major use of
spherical harmonics in medical imaging has been in representing the
coordinates of closed anatomical boundaries \cite{ bulow.2004,
geric.2002, shen04}. In this study, we use the spherical harmonics
for developing a kernel based smoothing technique for any type of
cortical data.

Previously cortical data have been smoothed by either solving
diffusion equations \cite{andrade.2001, cachia.2003b,
chung.cvpr.2003} or the iterative applications of the first order
heat kernel approximation \cite{chung.ipmi.2005}. However, the
diffusion smoothing approaches require setting up a finite element
scheme, which is computationally nontrivial, and making the
algorithm converges. The iterative kernel smoothing method is
simpler in comparison; however, since it is based on the repeated
applications of the first order approximation, the convergence is
very slow. To address the shortcomings of these methods, we propose
a new technique that construct the heat kernel analytically using
the spherical harmonics. The numerical implementation of the
smoothing is validated against the known analytic solutions.


\section{Spherical harmonics}
Consider the following parametrization of unit sphere $S^2$:
$$p =(x, y, z) = (\sin \theta \cos \varphi, \sin \theta \sin \varphi,
\cos \theta)$$ with $\theta \in [0,\pi], \varphi \in [0,2\pi)$. The
Laplacian corresponding to the parametrization is given by
$$\Delta_{S^2} = \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}.$$
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_{S^2} fh \;d\mu = \int_{0}^{2\pi} \int_0^{\pi} f(\theta,\varphi)h(\theta,\varphi) \sin \theta d\theta d\varphi,$$
where $d\mu = \sin \theta d\theta d\varphi$ is the area element. The
norm is defined as $\| f \|^2 = \langle f,f \rangle$. Consider the
following eigenequation $\Delta f + \lambda f =0.$ There are $2l+1$
eigenfunctions, denoted as  $Y_{lm} ( -l \leq m \leq l)$,
corresponding to eigenvalue $\lambda = l(l+1)$. $Y_{lm}$ is called
the {\em spherical harmonic} of degree $l$ and order $m$
\cite{courant}. The $2l+1$ spherical harmonics of degree $l$ are
 \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$ given by
$$P_l^m(x) = \frac{(1-x^2)^{m/2}}{2^l l!} \frac{d^{l+m} }{dx^{l+m}}
(x^2 -1)^l.$$ Unlike other previous studies \cite{bulow.2004,
geric.2002, shen04} that use the complex numbered spherical
harmonics, we are using a real function with different normalizing
constants.



Any $f \in L^2(S^2)$ can be expressed as \bq f=\sum_{l=0}^{\infty}
\sum_{m=-l}^{l} \hat f_{lm} Y_{lm}, \label{eq:fourier}\eq where the
Fourier coefficient $\hat f_{lm} = \langle f, Y_{lm} \rangle$. This
can be viewed as the general Fourier expansion of the Hilbert space
$L^2(S^2)$ on sphere. If we use up to the $d$-th degree spherical
harmonics for modeling the cortical data, we need total $(2\cdot 0 +
1) + \cdots + (2\cdot d+1) =(d+1)^2$ number of basis functions.

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{heat.kernel.shape.eps}
\caption{\label{fig:heat.kernel}Shape of the heat kernel with
different bandwidth $\sigma =0.01, 0.02, 0.05, 0.1, 0.5$ from top to
bottom. The horizontal axis is the $\theta$ value from the north
pole ($\theta=0$) to the south pole ($\theta=\pi$). As $\sigma$
becomes large, the heat kernel converges to constant value
$1/4\pi$.}
\end{figure}




\section{Heat kernel}

By extending Mercer's theorem \cite{courant} to $L^2(S^2)$, any
continuous symmetric non-negative definite kernel $K(p,q)$ can be
written as \bqn K(p,q) = \sum_{l=0}^{\infty} \sum_{m=-l}^l c_{lm}
Y_{lm}(p)Y_{lm}(q) \label{eq:mercer}\eqn for some constants
$c_{lm}$. The Green's function of the isotropic diffusion equation
\bqn \frac{\partial f}{\partial t} = \Delta f\label{eq:iso.diff}\eqn
in $S^2$ is called the {\em heat kernel} \cite{chung.ipmi.2005} and
it can be written in the Mercer's theorem with $c_{lm} =
e^{-l(l+1)\sigma}$. Then the heat kernel is given by \bqn
K_{\sigma}(p,q) = \sum_{l=0}^{\infty} \big[e^{-l(l+1)\sigma}
\sum_{m=-l}^{l} Y_{lm}(p)Y_{lm}(q)\big] \label{eq:hk}\eqn with
$\sigma = \sqrt{2t}$. It directly generalize the Gaussian kernel in
the Euclidean space to $S^2$ \cite{chung.ipmi.2005}. The kernel can
be also seen as the transition probability of the isotropic
diffusion process in $S^2$ and it should be symmetric along the
geodesic circle centered around $p$, defined by $\mathcal{C}_p =
\{q: p \cdot q = const.\}$, where we used the vector product $\cdot$
in the Cartesian coordinate system. This implies that the heat
kernel is isotropic with respect to $p$. Hence, if $p \cdot q = p
\cdot r$, we have $K_{\sigma}(p,q) = K_{\sigma}(p,r)$. This property
can be used to simplify the expansion (\ref{eq:hk}). Let us fix
$\varphi=0$ and $p$ be the north pole, i.e. $p=(0,0,1)$. Now by
varying $q= (\sin \theta, 0, \cos \theta)$ for $0 \leq \theta =
\cos^{-1} (p \cdot q) \leq \pi$, we have $Y_{lm} =0$ if $m \neq 0$.
Then the inner summation in (\ref{eq:hk}) can be written as\\
{\bf Theorem. 1.} \bqn \sum_{m=-l}^l Y_{lm}(p)Y_{lm}(q) =
\frac{2l+1}{4\pi}P_l^0(\cos \theta).\label{eq:addition} \eqn This
implies that the sum of products of the spherical harmonics is a
function of the angle between $p$ and $q$ only. For any $p$ other
than the north pole, the left hand side is also a function of the
angle between $p$ and $q$. The relationship (\ref{eq:addition}) is
usually called the harmonic addition theorem \cite{groemer}. Then we
have the following theorem:\\ {\bf Theorem. 2.} For any $p$ and $q$
in $S^2$, \bqn K_{\sigma}(p,q) = \sum_{l=0}^{\infty}
\frac{2l+1}{4\pi}e^{-l(l+1)\sigma}P_l^0(p \cdot
q).\label{eq:theorem} \eqn

Based on formulation (\ref{eq:theorem}), the shape of heat kernel is
plotted in Figure \ref{fig:heat.kernel} with varying bandwidths
$\sigma =$ 0.01, 0.02, 0.05, 0.10, 0.50. As $\sigma$ increases, the
kernel flattens out to a constant value, i.e. $\lim_{\sigma \to
\infty} K_{\sigma} (p,q) = \frac{1}{4\pi}$.

The full width at the half maximum (FWHM) of kernel has been widely
used as a unit for measuring the amount of Gaussian kernel
smoothing. For the usual 2D Gaussian kernel of the form
$\frac{1}{2\pi \sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}}$, it is
given by FWHM $=\sqrt{8\ln 2} \sigma$. For heat kernel defined in
$S^2$, since we can not find FWHM analytically in a close form, we
estimate it numerically. The maximum of the heat kernel is obtained
when $\theta = 0$ in (\ref{eq:theorem}). Note $P_l^0(1) = 1$ for all
$l$. Then we solve numerically for $\theta$ in
$$ \sum_{l=0}^{d} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}P_l^0(\cos \theta)
= \frac{1}{2}\sum_{l=0}^{d} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}.$$
The FWHM for the heat kernel is then $2\theta$. We used $d=20$
degree spherical harmonics in estimating the FWHM (Figure 2).


%This
%grantees the accuracy up to $10^{-3}$ for $\sigma > 0.02$.


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{FWHM.eps}
\caption{\label{fig:FWHM} Plot of FWHM (vertical) vs. bandwidth
$\sigma$ (horizontal). The blue line is for the heat kernel and red
line is for the isotropic Gaussian kernel in $\mathbb{R}^2$.}
\end{figure}

\section{Heat kernel smoothing}

The solution to the diffusion equation (\ref{eq:iso.diff}) with
initial functional measurement $h(p)$ at time $t=\sigma^2/2$ is
given by the convolution: \bqn K_{\sigma}*h(p) &=&
\int_{S^2} K_{\sigma}(p,q)h(q) \;d\mu(q)\\
&=&  \sum_{l=0}^{\infty} w_j
\mathcal{P}_l(p)\label{eq:convolution},\eqn where the weighting
factor $w_l = (2l+1)e^{-l(l+1)\sigma}$ and integral
$$ \mathcal{P}_l(q)=\frac{1}{4\pi}\int_{S^2} P_l^0(p\cdot
q)h(q)\;d\mu(q).$$ The heat kernel smoothing can be given a
geometric interpretation in the following way. Since the associated
Legendre polynomial $P_l^0$ for the first two terms are $P_0^0(p
\cdot q_k)=1$ and $P_1^0(p \cdot q_k)=p \cdot q_k$, we can see that
the first term
$$\mathcal{P}_0(p)
=\frac{1}{4\pi}\int_{S^2} h(q)\;d\mu(q)$$ measures the average
signal over the unit sphere while the second term $$\mathcal{P}_1
(p) =\frac{1}{4\pi}\int_{S^2} p \cdot q h(q)\;d\mu(q)$$ measures the
first moment of the signal with respect to $p$. The higher order
terms of $\mathcal{P}_l(p)$ can be computed using the recurrence
relation for the associated Legendre polynomial \bqn P_{l+1}^0(p
\cdot q) = \frac{2l+1}{l+1}p \cdot q P_l^0(p \cdot q) -
\frac{l}{l+1}P_{l-1}^0(p \cdot q) \label{eq:rec} \eqn so
$P_l^0(p\cdot q)$ is the $l$-th degree polynomial of $p \cdot q$.
Hence $\mathcal{P}_l(p)$ measures a signal quantity involving higher
order moments.

%Given functional measurement $h(p)$ in $S^2$, we model it
%stochastically as $h(p) = m(p) + \epsilon(p)$, where $m(p)$ is the
%unknown mean function and $\epsilon(p)$ is a zero-mean random field,
%possibly a Gaussian white noise process. Then we define the {\em
%heat kernel smoothing estimator} of $m(p)$ to be the above
%convolution, i.e. $\widehat m(p) = K_{\sigma} * h(p)$.


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.285]{thickness.smoothing.eps}
\caption{\label{fig:FWHM} Heat kernel smoothing on cortical
thickness. (a) original cortical thickness data mapped onto a unit
sphere, (b, c) smoothing with $d=40$, and $\sigma =0.0001, 0.001$
respectively. (d, e, f) smoothing with $d=20$, and $\sigma = 0.001,
0.01, 0.1$ respectively.}
\end{figure}

\section{Numerical Implementation}
The $S^2$ surface is realized as a mesh with 81,920 triangles. It is
constructed from the deformable surface algorithm that gives a
direct homological map from the human cortical surface to $S^2$
\cite{chung.cvpr.2003, macdonald.2000}. Integral of form $\int_{S^2}
f(q) \;d\mu(q)$ is approximated as a limit of the Riemann sum over
triangle elements following the framework of finite element method
(FEM) \cite{chung.cvpr.2003}. Let $n=40,962$ be the total number of
nodes in the mesh. Suppose that triangular elements $T_{k_1},
\cdots, T_{k_m}$ are adjacent to each other at given node $q_k$
($1\leq k \leq n$). Then the total area of the mesh is computed as
$$\int_{S^2} \;d \mu \doteq \frac{1}{3}\sum_{k=1}^n \sum_{i=1}^{k_m} |T_{k_i}| = 12.565$$
while the area of the unit sphere is $4\pi = 12.566$, the difference
of less than $0.0001 \%$. So our triangular mesh is sufficiently
fine enough to realize the $S^2$ surface accurately. The accuracy of
our heat kernel smoothing is only restricted to the mesh resolution
and the Riemann sum approximation should converges to the integral
as the mesh resolution decreases.

In our numerical implementation, we avoided using formulation
(\ref{eq:convolution}) due to the fact that the integral kernel
$P_l^0(p \cdot q)$ is not separable. In this case,  the running time
is $O(n^2)$ for smoothing on every nodes in the mesh. On the other
hand, if we use the form of the separable kernel (\ref{eq:hk}), the
smoothing can be done in $O(n)$. So our $O(n)$ implementation is to
use the form
$$K_{\sigma}*h(p) = \sum_{l=0}^{d} \sum_{m=-l}^{l} e^{-l(l+1)\sigma}
Y_{lm}(p)\int_{S^2} Y_{lm}(q) h(q) \; d\mu(q).$$ The integral is
approximated as
$$\int_{S^2}f(q) \;d\mu(q) \doteq \frac{1}{3}\sum_{k=1}^n \sum_{i=1}^{k_m}f(q_k)|T_{k_i}|.$$
For an illustration, we mapped the cortical thickness data
\cite{chung.ipmi.2005} obtained from MRI onto a unit sphere. Then we
performed the smoothing with various bandwidths $\sigma$ and $d$
(Figure 3).

The numerical implementation is validated against analytical
solutions. Let $h=Y_{l'm'}$ be an analytic test function. Then
$K_{\sigma}*h(p)$ can be written as
$$\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)$$
$$=\sum_{l=0}^{\infty} \sum_{m=-l}^{l}
e^{-l(l+1)\sigma} Y_{lm}(p) \delta_{l l'} \delta_{m m'} =
e^{-l'(l'+1)\sigma}Y_{l'm'}(p),$$ where $\delta_{l l'}$ is the
Kroneker's delta. This exact analytical form gives a criteria for
estimating the accuracy of our numerical method. The analytical
solution is compared with the numerical approximation. The table 1
shows the result for selective $l'$ and $m'$ with $\sigma=0.01$ and
$d=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. The fourth column shows
the average difference between the numerical approximation and the
expected theoretical value $e^{-l'(l'+1)\sigma}Y_{l'm'}$.

\begin{table}
\begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
$l'$ & $m'$ & integral & kernel smoothing\\
\hline

 1 & 1 & 0.9999 & $-3.8854 \cdot 10^{-17}$\\
10 & 5 & 0.9999 & $-3.7597 \cdot 10^{-17}$  \\
10     & 7 & 0.9999 &   $4.9109 \cdot 10^{-17}$  \\
15  & 10 & 0.9998 &  $-4.0601 \cdot 10^{-8}$  \\
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{Validation of heat kernel smoothing with test function
$Y_{l'm'}$.}
\end{table}
\section{Conclusions}
We have developed a mathematical framework for performing kernel
smoothing in a unit sphere. The heat kernel was constructed
analytically using the spherical harmonics. Based on the FWHM vs
$\sigma$ plot in Figure 2, we conclude that the heat kernel flattens
out more rapidly than the corresponding Gaussian kernel. Based on
the exact mathematical formulation, we developed the $O(n)$
smoothing algorithm. The algorithm is validated against the
analytical solution.



\bibliographystyle{IEEEbib}
%\bibliography{reference.OCT.2005}

\begin{thebibliography}{10}

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


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

\bibitem{chung.cvpr.2003}
M.K. Chung, K.J. Worsley, S.~Robbins, and A.C. Evans,
\newblock ``Tensor-based brain surface modeling and analysis,''
\newblock in {\em IEEE Conference on Computer Vision and Pattern Recognition
  (CVPR)}, 2003, vol.~I, pp. 467--473.

\bibitem{chung.ipmi.2005}
M.K. Chung, S.~Robbins, A.C. Evans,
\newblock ``Unified statistical approach to cortical thickness analysis,''
\newblock {\em Information Processing in Medical Imaging (IPMI)},
 in Lecture Notes in Computer Science (LNCS) 3565:627-638.
Springer-Verlag. 2005.



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

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

\bibitem{geric.2002}
M.S. Gudio~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{groemer} Groemer, H., {\it Geometric Applications of Fourier
Series and Spherical Harmonics,} Cambridge University Press, 1996.

\bibitem{hurdal.1999}
M.K., Hurdal, P.L., Bowers, et al.
\newblock ``Quasi-conformally flat mapping the human cerebellum,''
\newblock in {\em Medical Image Computing and Computer-Assisted
Intervention (MICCAI)}, 1999.

\bibitem{macdonald.2000}
J.D. MacDonald, N.~Kabani, D.~Avis, and A.C. Evans,
\newblock ``Automated 3-d extraction of inner and outer surfaces of cerebral
  cortex from mri,''
\newblock {\em NeuroImage}, vol. 12, pp. 340--356, 2000.

\bibitem{rosenberg}
S.~Rosenberg,
\newblock {\em The Laplacian on a Riemannian Manifold},
\newblock Cambridge University Press, 1997.


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

\end{thebibliography}
\end{document}
