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


% 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}{-3cm} \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 2004 IEEE International Symposium on Biomedical
Imaging Paper No. 562\\
\huge Diffusion Smoothing on Brain Surface via Finite Element
Method}
\author{Moo K. Chung\\University of Wisconsin\\
       Department of Statistics, Biostatistics and Medical Informatics\\
        Madison, WI 53706\\
     $\tt{mchung@stat.wisc.edu}$
      \and  Jonathan E. Taylor\\Stanford University\\
        Department of Statistics\\
        Stanford, CA 94305\\
      $\tt{jonathan.taylor@stanford.edu}$}
\maketitle \thispagestyle{empty}


%\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}
\subsection*{Introduction} In order to perform statistical
analysis on the brain cortical surface, data are usually diffused
to increase the signal-to-noise ratio and Gaussianness
\cite{Andrade, Chung:03a}. Most of the diffusion equation
approaches for polygonal surfaces are based on the finite element
method (FEM) without the explicit representation of the
Laplace-Beltrami operator.

We present two new methods for solving the diffusion equations on
the brain manifolds. The first method uses the explicit
representation of the Laplace-Beltrami operator from and the
diffusion equation is solved via a simple finite difference scheme
that speed up computation. The second methods uses the heat kernel
on the brain manifolds and the diffusion equations are solved via
iterative kernel smoothing.

\subsection*{Laplace-Beltrami Operator}
Consider orientable smooth twice-differentiable 2-dimensional
surface $\partial \Omega \subset {\mathbb{R}}^3$. Then we have a
parameterization of $\partial \Omega$: $X: D \to
\partial \Omega$ for some planar domain $D \subset \mathbb{R}^2$. Let $T_p(\partial \Omega)$
be a tangent space at any $p =X(u) \in \partial \Omega$ such that
partial derivatives
$$X_1(u)=\partial_{u^1} X(u), X_2(u)= \partial_{u^2} X(u)$$ form a
basis in $T_p(\partial \Omega)$. The inner products $g_{ij} =
\langle X_i, X_j \rangle$ are the Riemannian metric tensors. Then
the isotropic diffusion equation on $\partial \Omega$ is given by
\bqn \partial_t F = \Delta F\label{eq:diffusion}\eqn with initial
condition $F(p,0)=f(p)$ and the {\em Laplace-Beltrami operator}
\bq \Delta F = \frac{1}{|g|^{1/2}}\sum_{i,j=1}^2
\der{u^i}{}\Big(|g|^{1/2}g^{ij}\der{u^j}{F}\Big).
\label{eq:Beltrami}\eq Since the Laplace-Beltrami operator is {\em
self adjoint}  with respect to the $L^2(\partial \Omega)$ norm
$\langle F,G\rangle_{\partial \Omega} = \int_{\partial \Omega}
F(p)F(p)\;d\mu(p),$ we have \bqn \langle G,\Delta F
\rangle_{\partial \Omega}= -\int_{\partial \Omega} \langle \nabla
F,\nabla G \rangle \; d\mu(p) = \langle F,\Delta G
\rangle_{\partial \Omega}. \label{eq:adjoint}\eqn

\begin{figure}
\includegraphics[height=140mm]{Figure1.eps}
\caption{A typical triangular elements.} \label{fig:elements}
\end{figure}

\subsection*{Finite Element Method} Let $N_T$ be the number of
triangles in the triangular mesh $S$ that is the discrete
realization of true manifold $\partial \Omega$ (Figure
\ref{fig:thickness}). We seek solution $F_i$ to the diffusion
equation in each triangle $T_i$ such that the solution $F_i(p,t)$
is continuous across neighboring triangles, i.e. piecewise linear
function. The solution $F$ for $S$ is then \bqn F(p,t) \doteq
\sum_{i=1}^{N_T} F_i(p,t)\label{eq:FEM}.\eqn  Let
$p_{i_1},p_{i_2},p_{i_3}$ be the vertices of element $T_i$. In
$T_i$, we interpolate $F_i$ linearly by \bqn F_i(p,t) =
\sum_{j=1}^3 \xi_{i_j}(p)F(p_{i_j},t),\label{eq:interpolation}\eqn
where nonnegative $\xi_{i_k}$ are given by the {\em barycentric
coordinates} \cite{Sadiku}. Let $G$ and $F(x,t)$ be arbitrary
piecewise linear functions given by the barycentric coordinate
representation with $G_{i_k} = G(p_{i_k})$ and
$F_{i_k}=F(p_{i_k},t)$. Then from (\ref{eq:adjoint}), the
diffusion equation can be written as
 \bqn \langle G, \partial_t F \rangle_{T_i} =
 -\int_{T_i} \langle \nabla F, \nabla G \rangle \; d\mu(p)\label{eq:L2norm}\eqn
The above integral equation can be written in matrix form: \bqn
[G_i]'[A^i]\frac{d}{dt}[F_i] = - [G_i]'[C^i][F_i].\label{eq:G}
\eqn where $[G_i] =(G_{i_1},G_{i_2},G_{i_3})'$,
$[F_i]=(F_{i_1},F_{i_2},F_{i_3})'$, $[A^i]=(A_{kl}^i), A_{kl}^i =
\int_{T_i} \xi_{i_k}\xi_{i_l} \; d\mu(p)$, $[C^i]=(C_{kl}^i),
C_{kl}^i = \int_{T_i} \langle \nabla \xi_{i_k},\nabla \xi_{i_l}
\rangle \; d\mu(p).$ Since equation (\ref{eq:G}) should be true
for all $[G_i]$, we have a system of ordinary differential
equations (ODE) \bqn \frac{d[F_i]}{dt} = -[A^i]^{-1}[C^i][F_i]
\mbox{ for all } i\label{eq:triangle}.\eqn for each element $T_i$.
Having discretized an element, the next step is to assemble all
such elements in $m$ incident triangles around vertex $p$ (Figure
1). Combining $m$ elements,
 \bqn \frac{d[F]}{dt} = -
[A]^{-1}[C][F].\label{eq:ODE2}\eqn
  where $[F]=[F(p,t),
F(p_1,t),\cdots,F(p_m,t)]'$ and
$[G]=[G(p),G(p_1),\cdots,G(p_m)]'$. Matrix $[C]=(C_{ij})$ is
called the global coefficient matrix and comes from all elements
containing vertices $i$ and $j$. It seems we need to solve a huge
system of linear equations iteratively; however the first row of
the simultaneous ODE (\ref{eq:ODE2}) gives the diffusion equation
at the vertex $p=p_0$: \bqn \frac{dF(p,t)}{dt} = - \sum_{i,k=0}^m
A_{0k}^{-1}C_{ki}F(p_i,t), \label{eq:FEMlaplace}\eqn where
$A^{-1}_{0k}$ is the $0k$-th element of $A^{-1}$. Comparing this
with equation (\ref{eq:diffusion}), we can see the right-hand side
of equation (\ref{eq:FEMlaplace}) is the discrete estimation of
the Laplace-Beltrami operator at vertex $p$. Simplifying the
matrix inversion using the computational algebraic system
$\tt{MAPLE}$, we have the FEM estimation for the Laplace-Beltrami
operator given by \bqn \widehat \Delta F(p) = \sum_{i=1}^m
w_i\big(F(p_i)-F(p)\big)\label{eq:femestimate}\eqn with the
weights $w_i= (\cot \theta_i + \cot \phi_i)/|T|,$
 where  $\theta_i$ and $\phi_i$ are the two angles opposite to
 the edge $p_i-p$ and $|T|=\sum_{i=1}^m |T_i|$ is the sum of the
 areas of the incident triangles (Figure 1). A similar discrete representation that is based on geometric
arguments can be found in \cite{Meyer}.  The diffusion equation is
then solved by the finite difference scheme: \bqn F(p,t_{n+1})  =
F(p,t_n)+ \delta t \widehat{\Delta} F(p,t_n)\label{eq:difference}
\eqn with the iteration step size $\delta t =t_{n+1}-t_n$. The
convergence condition can be found in \cite{Chung:03b}. Note that
the Laplace-Beltrami operator in the conformal coordinate system
$(u^1,u^2)$ can be written as $\Delta = \dder{(u^1)}{} +
\dder{(u^2)}{}.$ So we can define the FWHM of diffusion smoothing
locally as the FWHM of the corresponding Gaussian kernel in the
conformal coordinate system. Then diffusion smoothing with $N$
iterations and the step size $\delta t$ would be equivalent to
Gaussian kernel smoothing with
$$\mbox{FWHM} =4(\ln2)^{1/2}\sqrt{N\delta t}.$$

\begin{figure}
\centering
\includegraphics[width=230mm]{diffusion_sphere.eps}
\caption{Diffusion smoothing was applied to smooth out the mean
curvature and projected onto a sphere to show how the hidden
sulcal pattern can be enhanced over time. (a) initial mean
curvature. (b) after 20 iterations with $\delta t =0.2$. (c) after
100 iterations.} \label{fig:sphere}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=240mm]{thickness_statistics5.eps}
\caption{Left, middle: Quantitle-Quantile(QQ) plot of cortical
thickness dilatation measurements before and after diffusions
smoothing for 28 subjects \cite{Chung:03a}. The horizontal axis
displays the quantiles of Gaussian distribution while the vertical
axis displays the quantiles of an empirical distribution. How
closely the blue dots lie along the straight red line gives an
idea if the underlying empirical distribution follows Gaussian .
Right: Lilliefors statistic measures the maximum difference
between the empirical and a theoretical Gaussian distributions.
Most of cortex shows value less than the cutoff value 0.19
indicating that the data have been smoothed to follow Gaussian
after.} \label{fig:stat}
\end{figure}
\subsection*{Iterated Heat Kernel Smoothing}
The drawback of the FEM approach is for larger $\delta t$, finite
difference scheme (\ref{eq:difference}) may diverge
\cite{Chung:03b}. We have developed a completely different method
based on the heat kernel construction that avoids this problem.
The solution to diffusion equation (\ref{eq:diffusion}) when
$t=\sigma^2/2$ is given by the convolution
$$F(p,\sigma) = K_{\sigma} *f (p)=\int_{\partial \Omega} K_{\sigma}(p,q)f(q)\mu(q)$$
where heat kernel $K_{\sigma}$ is given by the {\em parametrix
expansion} \cite{Rosenberg}: \bq K_{\sigma}(p,q) &=&
\frac{1}{(2\pi\sigma)^{1/2}}\exp
\big[-\frac{d^2(p,q)}{2\sigma^2}\big]\big[u_0(p,q) +
O(\sigma^2)\big]\eq and $d(p,q)$ is the geodesic distance between
$x$ and $y$ and $u_0(p,q) \to 1$ as $p \to q$. Under some
regularity condition, $K_{\sigma}$ is a probability distribution
on $\partial \Omega$. Note that $K_{\sigma}* (K_{\sigma} * f)$ is
the diffusion of signal $f$ after time $\sigma^2$ so that
$K_{\sigma}* (K_{\sigma} * f) = K_{\sqrt 2 \sigma} * f.$ Arguing
inductively, we have an iterated kernel smoothing formula
$$K_{\sigma}^{(m)} * f =\underbrace{K_{\sigma} *\cdots * K_{\sigma}}_m *f = K_{\sqrt m \sigma} * f.$$
So the kernel smoothing with large bandwidth can be performed
iteratively with a smaller bandwidth. In practice, we use
truncated and normalized kernel
$$\widetilde{K}_{\sigma}(p,q) =
\frac{\exp\big[-\frac{d^2(p,q)}{2\sigma^2}\big]{\bf 1}_B(p,q)}
{\int_{B} \exp\big[-\frac{d^2(p,q)}{2\sigma^2}\big]
\;d\mu(p)d\mu(q)}$$ for some small region $B \subset \partial
\Omega$ and indicator function ${\bf 1}_B$, i.e. $B=T_1 \cup
\cdots \cup T_m$. Then our iterated heat kernel smoothing of data
$f$ is given by
$$\widetilde{K_{\sigma}}^{(m)}*f=\underbrace{\widetilde{K_{\sigma}} *\cdots * \widetilde{K_{\sigma}}}_m
*f$$ which is the integral version of the Nadaraya-Watson kernel
estimator \cite{Chung:04e}.
\begin{figure}
\centering
\includegraphics[width=240mm]{diffusion_thickness2.eps}
\caption{Top: cortical thickness computed at the posterior right
hemisphere of the autistic brain and its iterated heat kernel
smoothing with $\sigma = 1$ and $m=100$ iterations. Bottom:
QQ-plot of a single subject 108,588 thickness measurements showing
increased Gaussianness after smoothing.} \label{fig:thickness}
\end{figure}
\subsection*{Conclusions} Based on the FEM, we discretized a
diffusion equation in a triangular mesh patch centered around a
vertex and solved a system of linear equation. It turns out that
the Laplace-Beltrami operator can be represented as a weighted
averaging where the weights are given in terms of the geometry of
triangular mesh elements. Then the diffusion equation is solved
via the finite difference with a temporal step size that satisfies
a convergence criterion. An alternate method is to iteratively
convolve with the heat kernel of the Laplace-Beltrami operator
with small bandwidth. Our diffusion and much simpler iterated heat
kernel smoothing would be highly useful in smoothing fMRI data
\cite{Andrade} and anatomical data \cite{Chung:03a}. Afterwards
statistical inference based can be based on the random fields
theory.
\subsection*{Acknowledgements} Authors wish to thank
Keith Worsley of McGill University for the motivation of the
problem. Authors also wish to thank Montreal neurological
institute for providing the mesh used in Figure 2 and the Keck
laboratory for functional imaging and behavior of the University
of Wisconsin-Madison for the MRI used in Figure 4. The mesh used
in Figure 4 was generated using FreeSurfer \cite{Dale}.

\tiny
\begin{thebibliography}{10}
\bibitem{Andrade}
A.~Andrade, F.~Kherif, J.~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}, {\bf 12}:79--93, 2001.

\bibitem{Chung:04e} M.K. Chung,
\newblock ``Iterated heat kernel smoothing on manifolds,''
\newblock Department of Statistics, University of
Wisconsin-Madison. Technical Report, 2004.

\bibitem{Chung:03a}
M.K. Chung, K.J. Worsley, S.~Robbins, T.~Paus, J.~Taylor,
J.~Giedd, J.L.
  Rapoport, and A.C. Evans,
\newblock ``Deformation-based surface morphometry applied to gray matter
  deformation,''
\newblock {\em NeuroImage}, {\bf 18}:198--213, 2003.

\bibitem{Chung:03b}
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},
  2003, vol.~I, pp. 467--473.

\bibitem{Dale}
A.M. Dale, B. Fischl, M.I. Sereno, \newblock ``Cortical
surface-based analysis I: segmentation and surface
reconstruction,'' \
\newblock {\em NeuroImage} {\bf 9}:179--194, 1999.

\bibitem{Meyer}
M.~Meyer, M.~Desbrun, P.~Schröder, and A.~H. Barr,
\newblock ``Discrete differential-geometry operators for triangulated
  2-manifolds,''
\newblock in {\em VisMath}, 2002.

%\bibitem{Malladi}
%R~Malladi and I.~Ravve,
%\newblock ``Fast difference schemes for edge enhancing beltrami flow,''
%\newblock in {\em Computer Vision-ECCV}, 2002, vol.~I, pp. 343--357.

%\bibitem{Perona}
%P.~Perona and J.~Malik,
%\newblock ``Scale-space and edge detection using anisotropic diffusion,''
%\newblock {\em IEEE Trans. Pattern Analysis and Machine Intelligence}, vol. 12,
%  pp. 629--639, 1990.

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

\bibitem{Sadiku}
M.N.O. Sadiku,
\newblock {\em Numerical Techniques in Electromagnetics},
\newblock CRC Press, 1992.

%\bibitem{Sochen}
%N.~Sochen, R.~Kimmel, and R.~Malladi,
%\newblock ``A general framework for low level vision,''
%\newblock {\em IEEE Transactions on Image Processing}, vol. 7, pp. 310--318,
%  1998.

%\bibitem{Tang:99}
%B.~Tang, G.~Sapiro, and V~Caselles,
%\newblock ``Direction diffusion,''
%\newblock in {\em ICCV}, 1999, pp. 1245--1452.

\end{thebibliography}

\end{multicols}
\end{document}
