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

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

\title{DIFFUSION SMOOTHING ON BRAIN SURFACE VIA FINITE ELEMENT METHOD}

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

\twoauthors {Moo K. Chung}
         {University of Wisconsin\\
         Department of Statistics\\
         Madison, WI 53706\\
         $\tt{mchung@stat.wisc.edu}$}
 {Jonathan Taylor}
         {Stanford University\\
         Department of Statistics\\
         Stanford, CA 94305\\
         $\tt{jonathan.taylor@stanford.edu}$}

\begin{document}

\maketitle
%

\begin{abstract}
Surface data such as the segmented cortical surface of the human
brain plays an important role in medical imaging. To increase the
signal-to-noise ratio for data residing on the brain surface, the
data is usually diffused. Most of diffusion equation approach for
triangulated mesh data is based on the finite element method and a
system of linear equations are iteratively solved without the
explicit representation of the Laplace-Beltrami operator. Such
implicit formulation requires inverting large sparse matrix that
is required in most finite element methods.

The novelty of our paper is in the explicit representation of the
Laplace-Beltrami operator derived from the finite element method
itself. The Laplace-Beltrami operator is expressed as a weighted
averaging operator where the weights are expressed in terms of the
interior angles and the area of triangles. The weights are
computed by actually solving the system of linear equations in the
finite element method and inverting a matrix in a computational
algebra system $\tt{MAPLE}$. Afterwards the diffusion equation is
solved via a simple finite difference scheme that speed up
computation.

\end{abstract}

\section{Introduction}
\label{sec:intro}
 When solving a diffusion equation on a curved surface as a way to
increase the signal-to-noise ratio (SNR), the Laplacian somehow
has to incorporate the geometry of the curved surface. The
extension of the Euclidean Laplacian to an arbitrary Riemannian
manifold is called the {\em Laplace-Beltrami operator}
\cite{Kreyszig}. Although diffusion equations are widely used to
smooth out data in imaging analysis
\cite{Malladi,Perona,Sochen,Tang:99,Taubin}, there is a very few
paper that uses diffusion to smooth out data residing on the brain
surface \cite{Andrade, Chung:03a}. Many previous work involving
diffusion are the problem of surface fairing \cite{Meyer, Taubin}
or anisotropic smoothing of intensity images where the image
intensities are taken as surfaces to be smoothed \cite{Sochen}.

In this paper, we show how to represent the Laplace-Beltrami
operator explicitly for an arbitrary triangular surface mesh using
the finite element method (FEM). Afterwards the explicit
representation of the Laplacian is used in the finite difference
method (FDM) for solving a diffusion equation iteratively on the
human brain surface for the cortical thickness and curvature
measurements \cite{Chung:03a}. A similar discrete representation
that is based on geometric arguments can be found in
\cite{Meyer,Schneider}.

\section{Laplace-Beltrami Operator}

Suppose we have an orientable smooth twice-differentiable $(C^2)$
2-dimensional surface $S$ in ${\mathbb{R}}^3$. Then we have a
parameterization of $S$:
$$X(u) = \{ x_1(u),x_2(u),x_3(u): u=(u^1,u^2) \in D \}$$
for some planar domain $D$. Let $T_pS$ be a tangent space at any
$p =X(u) \in S$ 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$. Any vector $d\xi \in T_pS$ can be written as
$d\xi=du^1X_1 + du^2X_2$ for some constants $du^1$ and $du^2$ and
the length of the vector $d\xi$ is \bqn d\xi^2 \equiv \langle
d\xi,d\xi \rangle = \sum_{i,j}g_{ij}du^i du^j,
\label{eq:bilinear}\eqn where the inner products $g_{ij} = \langle
X_i, X_j \rangle$ are called the {\it Riemannian metric tensor}.
Then the Laplace-Beltrami operator $\Delta_X$ corresponding to the
surface parameterization $X$ is defined as
 \bqn \Delta_X 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}\eqn Since the Laplace-Beltrami operator is
{\em self adjoint}  with respect to the $L^2$ norm for any $C^2$
functions $F$ and $G$ on $S$, it can be written as \bqn \int_S
G\Delta F dS= -\int_S \langle \nabla F,\nabla G \rangle \; dS =
\int_S F\Delta G dS. \label{eq:adjoint}\eqn

\section{Finite Element Method}
Let $N_T$ be the number of triangles in the triangular mesh $\hat
S$ that is the discrete estimation of true surface $S$. We seek an
approximate solution $F_i$ in triangle $T_i$ such that the
solution $F_i(x,t)$ is continuous across neighboring triangles.
The approximate solution $F$ for the whole region is then \bqn
F(x,t) \doteq \sum_{i=1}^{N_T} F_i(x,t)\label{eq:FEM}.\eqn A
slightly different formulation of FEM for the surface flattening
problem is given in \cite{Angenent}. 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(x,t) = \sum_{j=1}^3
\xi_{i_j}(x)F(p_{i_j},t),\label{eq:interpolation}\eqn where
nonnegative $\xi_{i_k}$ are given by the {\em barycentric
coordinates} \cite{Sadiku}. In the barycentric coordinates, any
point $x \in T_i$ is uniquely determined by two conditions: $$x
=\sum_{k=1}^3 \xi_{i_k}(x)p_{i_k}, \sum_{k=1}^3 \xi_{i_k}(x)= 1.$$
\subsection{Discrete Diffusion Equation}
Let $G$ be an arbitrary piecewise linear function given by
$$G(x) = \sum_{i=1}^{N_T} \xi_{i_1}(x)G_{i_1}+\xi_{i_2}(x)G_{i_2}+\xi_{i_3}(x)G_{i_3},$$
where $G_{i_k}=G(p_{i_k})$ are the values of function $G$
evaluated at vertex $p_{i_k}$ of $T_i$. Another piecewise linear
function $F(x,t)$ is given similarly so that $F_{i_k} =
F(p_{i_k},t)$.

Then from equation (\ref{eq:adjoint}), the integral version of
diffusion equation $\partial_t F = \Delta F$ can be written as
 \bqn
\int_{T_i} G\partial_t F \; dT = -\int_{T_i} \langle \nabla F,
\nabla G \rangle \; dT \label{eq:L2norm}.\eqn The left-hand term
in (\ref{eq:L2norm}) is \bq \int_{T_i} G
\partial_t F \; dT &=& \sum_{k,l=1}^3 G_{i_k}\partial_t F(p_{i_l},t)
\int_{T_i}
\xi_{i_k}\xi_{i_l} \;  dT \\
&=&[G_i]'[A^i] \frac{d}{dt} [F_i], \label{eq:Lhand}\eq where
$[G_i] =(G_{i_1},G_{i_2},G_{i_3})'$,
$[F_i]=(F_{i_1},F_{i_2},F_{i_3})'$ and $3 \times 3$ matrix
$[A^i]=(A_{kl}^i), A_{kl}^i = \int_{T_i} \xi_{i_k}\xi_{i_l} \;
dT$. It can be shown that
$$[A^i] = \frac{|T_i|}{12}\begin{pmatrix}
  2 & 1 & 1 \\
  1 & 2 & 1 \\
  1 & 1 & 2 \end{pmatrix},$$
where $|T_i|$ is the area of the triangular element $T_i$
\cite{Sadiku}. Similarly the right-hand term in (\ref{eq:L2norm})
can be written as \bq \int_{T_i}  \langle \nabla F, \nabla G
\rangle dT &=& \sum_{k,l=1}^3 G_{i_k}F(p_{i_l},t)\int_{T_i}
\langle \nabla \xi_{i_k},\nabla \xi_{i_l} \rangle dT\\
& =& [G_i]'[C^{i}][F_i], \label{eq:Rhand} \eq where
$[C^i]=(C_{kl}^i), C_{kl}^i = \int_{T_i} \langle \nabla
\xi_{i_k},\nabla \xi_{i_l} \rangle \; dT.$  Since $T_i$ is planar,
the gradient $\nabla \xi_{i_k}$ becomes the standard planar
gradient. Following \cite{Sadiku,Silvester}, it can be shown that
$2[C^i]$ is the matrix of the following form:
\begin{figure}[t]
\centering \centerline{\epsfig{figure=Figure1.eps,width=6.5cm}}
\caption{A typical triangular elements.} \label{fig:elements}
\end{figure}
$$\begin{pmatrix}
\cot \theta_{i_2}+\cot \theta_{i_3} & -\cot \theta_{i_3} & -\cot \theta_{i_2}\\
  -\cot \theta_{i_3} & \cot \theta_{i_1} + \cot \theta_{i_3} & -\cot \theta_{i_1}\\
-\cot \theta_{i_2} & -\cot \theta_{i_1} & \cot \theta_{i_1}+\cot
\theta_{i_2}\end{pmatrix},$$ where $\theta_{i_k}$ is the interior
angle of vertex $p_{i_k}$. Equating the left and right hand sides,
we get \bqn [G_i]'[A^i]\frac{d}{dt}[F_i] = -
[G_i]'[G^i][F_i].\label{eq:G} \eqn Since the equation (\ref{eq:G})
should be satisfied for an arbitrary vector $[G_i]$, we have a
system of ordinary differential equations (ODE) given by \bqn
\frac{d[F_i]}{dt} = -[A^i]^{-1}[C^i][F_i] \mbox{ for all }
i\label{eq:triangle}.\eqn

\subsection{Assembling Elements}
Having discretized an element, the next step is to assemble all
such elements in $m$ incident triangles around vertex $p$. Let
$p_1, \cdots, p_m$ be the
 $m$ neighboring vertices around $p=p_0$ in the counter-clockwise direction.
 Let $p,p_i,p_{i+1}$ be the vertices of the element $T_i$ (Figure 1). Then
 from combining elements, we have
\bq \int_{T_1 \cup \cdots \cup T_m} \langle \nabla F, \nabla G
\rangle \; dT &=&\sum_{i=1}^{m} \int_{T_i} \langle \nabla F,
\nabla G \rangle \; dT\\&=& [G]'[C][F],\eq 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)]'$. The matrix $[C]=(C_{ij})$ is
called the {\em global coefficient matrix}, which is the
assemblage of individual element coefficients. The contribution to
$C_{ij}$ comes from all elements containing vertices $i$ and $j$.
In the case of a hexagonal triangulation in Figure 1, $[C]$ is
given by
\begin{figure}[htb]
  \centering
\centerline{\epsfig{figure=diffusion_coefficient.eps,width=8.5cm}}
\end{figure}\\
Similarly
 \bq \int_{T_1 \cup \cdots \cup T_m}
G\der{t}{F} \; dT &=&[G]'[A]\frac{d[F]}{dt}, \eq where
$[A]=(A_{ij})$ has the same structure as $[C]$, i.e. write $A_{01}
= A_{01}^1 + A_{01}^6$ instead of $C_{01} = C_{01}^1 + C_{01}^6$.
Equating the above equations we have \bqn [G]'[A]\frac{d[F]}{dt} =
-[G]'[C][F]. \label{eq:discrete}\eqn Since equation
(\ref{eq:discrete}) should be satisfied for an arbitrary piecewise
linear function $G$, we have a discrete diffusion equation on $m$
elements $T_1,\cdots,T_m$ given by \bqn \frac{d[F]}{dt} = -
[A]^{-1}[C][F].\label{eq:ODE2}\eqn At first glance, it seems like
we need to solve a system of linear equations iteratively. Note
that the first row of the simultaneous ODE (\ref{eq:ODE2}) gives
the discrete 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 the diffusion equation
$\partial_t F(p,t) = \Delta F(p,t)$, we can see that the
right-hand side of equation (\ref{eq:FEMlaplace}) is the discrete
estimation of Laplacian of function $F$ evaluated at vertex $p$.
Simplifying the matrix computation 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).


\begin{figure}[t]
\centering
\centerline{\epsfig{figure=diffusion_sphere.eps,width=8.5cm}}
\caption{Diffusion smoothing was applied to smooth out the mean
curvature of the brain cortex and projected onto a sphere to show
how 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 we are beginning to see much clear sulcal
pattern.} \label{fig:sphere}
\end{figure}

\section{Finite Difference Method}
The diffusion equation is solved by the finite difference scheme:
\bqn F(p,t_{n+1})  = F(p,t_n)+ (t_{n+1}-t_n) \widehat{\Delta}
F(p,t_n)\label{eq:difference} \eqn where $\widehat{\Delta}
F(p,t_n)$ is estimated by (\ref{eq:femestimate}). We may fix the
iteration step size $t_{n+1}-t_n=\delta t$. For the convergence of
the finite difference scheme, $\delta t$ is chosen to satisfy the
harmonic condition \bqn \min_i F(p_i,t_n) \leq F(p,t_n+\delta t)
\leq \max_i F(p_i,t_n).\label{eq:inequality} \eqn For more
detailed convergence condition, see \cite{Chung:03b}.

Note that the Laplace-Beltrami operator in the conformal
coordinate system \cite{Kreyszig} $(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}.$$ %Computing the
%linear weights for the Laplace-Beltrami operator takes a fair
%amount of time in $\tt{MATLAB}$ but once the weights are computed,
%it is applied through the whole iteration repeatedly and the
%actual finite difference scheme takes less than one minute for 100
%iterations in Pentium 4.

\section{Conclusions}
Based on the FEM, we discretized a diffusion equation in a
triangular mesh patch centered around a vertex and solved a system
of equation using a computational algebraic system $\tt{MAPLE}$.
It turns out that the Laplace-Beltrami operator can be represented
as a weighted averaging operation where the weights are given in
terms of the interior angles and the areas of triangles. Based on
the FEM estimate for the Laplace-Beltrami operator, we iteratively
run the finite difference with a temporal step size that satisfies
a convergence criterion. This diffusion smoothing would be highly
useful in smoothing fMRI data \cite{Andrade} and anatomical data
\cite{Chung:03a} that are residing on the fixed brain surface and
would be useful in statistical inference based on the random
fields theory.

\section{Acknowledgements}
Most of work was done when authors were students of Keith Worsley
at McGill University. Authors wish to thank the Montreal
neurological institute for providing the MRI 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 3.

\begin{figure}[t]
\centering
\centerline{\epsfig{figure=diffusion_thickness.eps,width=8.5cm}}
\caption{Left: cortical thickness computed at the posterior right
hemisphere of the autistic brain. Right: Diffusion smoothing
estimation of the cortical thickness metric. We can see huge noise
reduction. Gyri have thicker gray matter compared to thinner
sulci.} \label{fig:thickness}
\end{figure}

\bibliographystyle{IEEEbib}
%\bibliography{surface}
\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}, vol. 12, pp. 79--93, 2001.


\bibitem{Angenent}
S.~Angenent, S.~Hacker, A.~Tannenbaum, and R.~Kikinis,
\newblock ``On the laplace-beltrami operator and brain surface flattening,''
\newblock {\em IEEE Transactions on Medical Imaging}, vol. 18, pp. 700--711,
  1999.

\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}, vol. 18, pp. 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{Kreyszig}
E.~Kreyszig,
\newblock {\em Differential Geometry},
\newblock University of Toronto Press, Toronto, 1959.


\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{Schneider}
R.~Schneider and L.~Kobbelt,
\newblock ``Geometric fairing of irregular meshes for freeform surface
  design,''
\newblock {\em Computer Aided Geometric Design}, vol. 18, 2001.


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

\bibitem{Silvester}
P.P. Silvester and R.L. Ferrari,
\newblock {\em Finite Elements for Electrical Engineers},
\newblock Cambridge: Cambridge Univ. Press, 1983.

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

\bibitem{Taubin}
G.~Taubin,
\newblock ``Curve and surface smoothing without shrinkage,''
\newblock in {\em The Proceedings of the Fifth International Conference on
  Computer Vision}, vol. 852-857, p. 1995.


\end{thebibliography}

\end{document}
