\documentclass[portrait,a0]{a0poster}
%120cm tall by 240cm wide (4ft. tall by 8ft. wide) (landscape
%orientation).
% 33 x 47 in  = 83cm  x 119 cm

\usepackage{color}
\usepackage{times}
\usepackage{a0poster}
\usepackage{graphicx}
\usepackage{multicol}
%\usepackage{pstricks,pst-grad}
%\usepackage[dvips]{graphicx}
\usepackage{times,amsmath,amsfonts}

%\usepackage{spconf,amsmath,epsfig}
\usepackage{amsfonts,amssymb}

% This is the sideways spacing between the columns of  text
\columnsep=70pt

% This is the thickness of the black line between the columns of text
\columnseprule=2pt

%\usepackage{pstricks}

%\psset{unit=1cm}

% Gives you the times font
%\usepackage{times}

% Alllows you to include postscript files. Instructions in
% /usr/local/teTeX/texmf/doc/latex/graphics/epslatex.ps
%\usepackage{graphicx}


%\setlength{\parindent}{0.0cm}

\addtolength{\textwidth}{-4cm} \addtolength{\oddsidemargin}{1cm}

\def\x{{\mathbf x}}
\def\L{{\cal L}}

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

\newcommand {\diff}[1] {\frac{\partial}{\partial #1}}
\newcommand{\jacob}[3]{\frac{\partial^2 #3}{\partial #1 \partial #2}}
\newcommand{\der}[2]{\frac{\partial #2}{\partial #1}}
\newcommand{\dder}[2]{\frac{\partial^2 #2}{\partial #1^2}}


\begin{document}




\title{\small 2006 IEEE International Symposium on Biomedical
Imaging (ISBI) Paper No. 1430\\
\huge \bf Heat Kernel Smoothing on Unit Sphere}
\author{ Moo K. Chung\\
       Department of Statistics, Biostatistics and Medical Informatics\\
       Waisman Laboratory for Brain Imaging and Behavior\\
       University of Wisconsin-Madison\\
     \tt{http://www.stat.wisc.edu/}$\sim$\tt{mchung}} \maketitle \thispagestyle{empty}

%\date{April 8, 2006}
%\author{Jonathan L. Marchini and Brian D. Ripley}
%\address{Department of Statistics, University of Oxford, UK$^\dagger$}
%\email{\{marchini,ripley\}@stats.ox.ac.uk}



\begin{multicols}{3}
 \large \subsection*{Introduction}

In brain imaging, cortical data such as cortical thickness, surface
curvatures and surface coordinates have been mapped to a unit sphere
for visualization, surface registration and data analysis.

To increase the signal-to-noise ratio, the cortical data have been
mainly 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 these problems, we propose a new technique that construct
the heat kernel analytically using the spherical harmonics.


\subsection*{Spherical Harmonics} Given the following parametrization of unit sphere $S^2$
$$p = (\sin \theta \cos \varphi, \sin \theta \sin \varphi,
\cos \theta),$$ the corresponding spherical Laplacian is given by
$$\Delta = \frac{1}{\sin \theta}\frac{\partial }{\partial
\theta}\Big(\sin \theta \frac{\partial}{\partial
 \theta}\Big) +\frac{1}{\sin^2 \theta} \frac{\partial^2 }{\partial^2
 \varphi}.$$
For $f,h \in L^2(S^2)$, the space of square integrable functions in
$S^2$, the inner product is
$$\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
orthonormal basis functions of $L^2(S^2)$ are given by the {\em
spherical harmonic} of degree $l$ and order $m$, denoted by
$Y_{lm}$:
 \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 $P_l^m$ is the associated Legendre polynomials of order $m$.
So any $f \in L^2(S^2)$ can be expressed as \bq
f=\sum_{l=0}^{\infty} \sum_{m=-l}^{l} f_{lm} Y_{lm},
\label{eq:fourier}\eq where the Fourier coefficient $f_{lm} =
\langle f, Y_{lm} \rangle$. This is the basis of the spherical
harmonic (SPHARM) representation used in computational neuroanatomy
\cite{geric.2002,shen04}.

\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=1\linewidth]{harmonics.eps}
\caption{Spherical harmonics of for degree $l=5, 30, 45$.}
\end{figure}





\subsection*{Heat Kernel}

The {\em heat kernel} or Gauss-Weistrass kernel
\cite{chung.ipmi.2005} is defined as $$ 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:hk}.$$ It directly
generalize the Gaussian kernel in the Euclidean space to $S^2$
\cite{chung.ipmi.2005}. The heat kernel can be written in more
compact form as $$ 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}$$



The shape of heat kernel is shown in Figure \ref{fig:heat.kernel}
for varying $\theta = \cos^{-1}(p\cdot q)$ with different bandwidths
$\sigma =$ 0.01, 0.02, 0.05, 0.10, 0.50.

\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=0.6\linewidth]{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 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}

The full width at the half maximum (FWHM) of kernel has been widely
used as a unit for measuring the amount of 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
$$\mbox{FWHM } =\sqrt{8\ln 2} \sigma.$$ For the heat kernel in
$S^2$, since we can not find FWHM analytically in a close form, we
estimate it numerically  by solving for $\theta$ in
$$ \sum_{l=0}^{d} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}P_l^0(\theta)
= \frac{1}{2}\sum_{l=0}^{d} \frac{2l+1}{4\pi}e^{-l(l+1)\sigma}$$
where $d$ is chosen before hand (Figure \ref{fig:FWHM}).

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


\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=0.6\linewidth]{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.}
\end{figure}

\subsection*{Heat Kernel Smoothing}

We define heat kernel smoothing of data $f$ to be \bq
K_{\sigma}*f(p) &=&
\int_{S^2} K_{\sigma}(p,q)f(q) \;d\mu(q)\\
&=&  \sum_{l=0}^{\infty} \sum_{m=-l}^l e^{-l(l+1)\sigma} f_{lm}
Y_{lm}(p) \label{eq:convolution}.\eq

Note that this is the solution to the isotropic diffusion equation
$$ \frac{\partial g}{\partial \sigma} = \Delta
g, g(\sigma=0,p)=f(p)\label{eq:iso.diff}.$$




\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=1\linewidth]{thickness.smoothing.eps}
\caption{\label{fig:thick} 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}



\subsection*{Numerical Implementation} The $S^2$ surface is realized
as a triangle mesh. 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}.
The Fourier coefficients are estimated using the numerical
integration technique. The total area of the mesh is 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. As an illustration,
we mapped the cortical thickness data \cite{chung.ipmi.2005}
obtained from MRI onto $S^2$ mesh. We performed heat kernel
smoothing with various bandwidths $\sigma$ on the thickness data
(Figure \ref{fig:thick}).

\subsection*{Applications}
As an application, we show how to estimate the cortical thickness by
reconstructing the cortex using the heat kernel smoothing technique.
The Cartesian coordinates of the both outer and inner surfaces are
mapped onto $S^2$ via the deformable surface algorithm (Figure
\ref{fig:xyz}).
\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=0.8\linewidth]{xyz3.eps}
\caption{\label{fig:xyz} Heat kernel smoothing of coordinates}
\end{figure}

Then the surfaces are reconstructed with varying $\sigma$ and up to
$l=80$ degrees of spherical harmonics (Figure \ref{fig:org}).
\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=1\linewidth]{cortex.eps}
\caption{\label{fig:org} Original cortex and its reconstruction at
different $\sigma=0, 0.0001, 0.001, 0.01$ with up to $l=80$ degree
harmonics. When $\sigma=0$, the heat kernel smoothing gives the
traditional SPHARM \cite{geric.2002, shen04}.}
\end{figure}

By averaging the Fourier coefficients of heat kernel smoothing, we
can construct the average cortical surface for 12 normal subjects
(Figure \ref{fig:avg}).
\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=1\linewidth]{average-fourier.eps}
\caption{\label{fig:avg} Average template constructed by averaging
the coefficients of heat kernel smoothing.}
\end{figure}

The cortical thickness is then estimated by taking the mean square
of the Fourier coefficients of heat kernel smoothing (Figure
\ref{fig:thick2}).
\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[width=1\linewidth]{thickness3.eps}
\caption{\label{fig:thick2} Cortical thickness estimation for
different $\sigma=0, 0.0001, 0.001, 0.01$.}
\end{figure}

\subsection*{Conclusions} We have developed a theoretical framework
for performing heat kernel smoothing on a unit sphere. The heat
kernel was constructed analytically using the spherical harmonics.


%\subsection*{Acknowledgements}

\tiny
\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{multicols}
\end{document}
