

% Make two column format for LaTex 2e.
%\documentclass[10pt,twocolumn]{article}
\documentclass{llncs}
\usepackage{llncsdoc}

\usepackage{epsfig,epsf,psfrag}
\usepackage{times,amsmath,amsfonts}
%\usepackage[dvips]{graphicx,graphics,color}
%\usepackage{times,amsmath,amsfonts}


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

%        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}
% 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}
\newcommand{\bq}{\begin{eqnarray*}}
\newcommand{\eq}{\end{eqnarray*}}


\begin{document}


\hspace{-0.6cm}DEPARTMENT OF STATISTICS\\
University of Wisconsin\\
1210 West Dayton St.\\
Madison, WI 53706

\vspace{2cm}
\begin{center}
\Large
{TECHNICAL REPORT NO. 1104}\\
\large June 20, 2005

\vspace{0.5cm} {\Large Anisotropic Kernel Smoothing in Diffusion\\
Tensor Imaging: Theoretical Framework}
% Make title bold and 14 pt font (Latex default is non-bold, 16pt)

\vspace{1cm}

$^{1,2,3}$Moo K. Chung, $^{3,4}$Jee Eun Lee, $^{3,4}$Andrew L.
Alexender \vspace{0.5cm}

$^1$Department of Statistics\\
 $^2$Department of Biostatistics and Medical Informatics\\
$^3$Waisman Laboratory for Brain Imaging and Behavior\\
$^4$Department of Medical Physics\\
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
%\date{}
\title{\Large \bf Anisotropic Kernel Smoothing in Diffusion\\
Tensor Imaging: Theoretical Framework}

% For single author (just remove % characters)
%\author{I. M. Anonymous \\
%  My Department \\
%  My Institute \\
%  My City, STATE, zip}

% For two authors (default example)
\author{$^{1,2,3}$Moo K. Chung, $^{3,4}$Jee Eun Lee, $^{3,4}$Andrew L. Alexender}
\institute{$^1$Department of Statistics, $^2$Department of Biostatistics and Medical Informatics\\
$^3$Waisman Laboratory for Brain Imaging and Behavior\\
$^4$Department of Medical Physics\\
University of Wisconsin-Madison\\
Madison, WI 53706}
%\date{November 10, 2003}
\maketitle

%\large
\section*{\centering Abstract}
{\em We present a theoretical framework for smoothing diffusion
tensor images. We formulate smoothing data along the white fiber
tracks as iterated anisotropic Gaussian kernel smoothing. The
formulation is derived from the construction of the transitional
probability of diffusion process. The smoothing is performed in such
a way that it smooth data more along the most likely direction of
water molecule diffusion direction. Statistical properties of our
smoothing are also presented.}

\section{Introduction}

Diffusion tensor imaging (DTI) is a new magnetic resonance imaging
technique technique that provides the directional information of
water diffusion in the white matter of the brain. The directional
information is usually represented as a symmetric positive definite
$3 \times 3$ matrix which is usually termed as the {\em diffusion
tensor} $D=(d_{ij})$. The water diffusion is assumed to follow
diffusion equation \bqn
 \frac{\partial C}{\partial t} = \nabla \cdot
(D \nabla C)\label{eq:diffusion1},\eqn where $C$ is the
concentration of water molecules \cite{kang.2004}.
%The diffusion
%tensor can be used to construct the kernel shape in kernel smoothing
%formulation.

The previous approaches in smoothing DTI data are usually formulated
as the anisotropic diffusion or the Laplace-Beltrami flow
\cite{Batchelor,Parker:2000,Spira,Tang,Weinstein}. In anisotropic
diffusion framework, the white fiber connectivity metric is given as
a solution to equation (\ref{eq:diffusion1}) with the initial
condition where every voxel is zero except a seed region where it is
given the value one \cite{Batchelor}. Then the value 1 at the seed
region is diffused though the brain and the numerical value that
should be between 0 and 1 is taken as a probability of connection.
Mathematically it is equivalent as the Monte-Carlo random walk
simulation with no restriction \cite{Koch}.

Our kernel smoothing approach is different from \cite{Batchelor} and
\cite{Koch}. In our kernel-based approach, we construct anisotropic
kernel that is the transition probability density of a continues
diffusion process and perform iterative anisotropic Gaussian kernel
smoothing. This formulation is simpler than solving equation
(\ref{eq:diffusion1}) numerically. Further, unlike the Monte-Carlo
random walk simulation approach, it always gives the deterministic
results that does not change from one run to another. Due to the
simplistic nature of our smoothing formulation, it is also easier to
derive its statistical properties.

As applications of our technique, we will present two examples. In
the first example, the transitional probability from a region in the
corpus callosum will be computed. In the second example, the
fractional anisotropy (FA) map will be smoothed out to show that our
smoothing will preserve FA while increasing the signal-to-noise
ratio (SNR).


\section{Transition probability}

Let $P_t({\bf p}, {\bf q})$ be the {\em transition probability
density} of a particle going from ${\bf p}$ to ${\bf q}$ under
anisotropic diffusion process.  This is the probability density of
the particle hitting ${\bf q}$ at time $t$ when the particle is at
${\bf p}$ at time 0. So we have $\int_{\mathbb{R}^n}P_t({\bf p},{\bf
x}) \; d{\bf x}= 1$ and $P_t({\bf p}, {\bf q}) = P_t({\bf q}, {\bf
p}).$ The transition probability serves as a mathematical analogue
of the water molecule concentration $C$ in equation
(\ref{eq:diffusion1}). The solution to equation
(\ref{eq:diffusion1}) with the initial condition $p(0,{\bf x}) =
\delta({\bf x})$, the Dirac delta function, has been used as the
probabilistic representation of white fiber track connectivity in
\cite{Batchelor}.


When $D$ is constant, it can be shown that $P_t({\bf p}, {\bf q}) =
K_t({\bf p} - {\bf q})$, the $n$-dimensional anisotropic Gaussian
kernel defined as follows. First, let ${\bf x}=(x_1,\cdots,x_n)' \in
\mathbb{R}^n$. The $n$-dimensional isotropic Gaussian kernel is
defined as $ K({\bf x})=(2\pi)^{-n/2}e^{-{\bf x}'{\bf x}/2}$ which
is the joint multivariate normal probability density function of $n$
independent standard normal random variables. Note that
$\int_{\mathbb{R}^n} K({\bf x}) \;d{\bf x}=1$. Then we define
anisotropic kernel $K_t$ via transformation $h: {\bf x} \to
(2\sqrt{t}D)^{-1/2}{\bf x}$ such that
$$K_t({\bf x}) = {\det}^{-1}\Big(\frac{\partial h}{\partial {\bf x}}\Big)K \circ h({\bf x})
= \frac{1}{(4\pi t)^{n/2} {\det}^{1/2} D} \exp \Big(-\frac{{\bf
x}D^{-1}{\bf x}}{4t}\Big).$$ The transformation matrix
$2\sqrt{t}D^{-1/2}$ is called the bandwidth matrix and controls the
concentration of the anisotropic kernel weights. $K_t$ is a
multivariate normal density with the mean zero and the covariance
matrix $2tD$.

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.1]{transition.prob3.eps}
\caption{\label{fig:transition} The transition probabilities from
the center voxel to neighboring voxels. They form the anisotropic
kernel weights. Left: transition probability without restriction.
Right: transition probability within  FA $> 0.4$.}
\end{figure}

\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.17]{transition.prob.eps}
\caption{\label{fig:trans.prob} The long distance transition
probability from the seed voxel in the corpus callosum with $\Delta
t = 0.1$ and $k=40,80, 120, 160$ iterations. The scale in the top
images is $10^{-3}$. The bottom images are in natural log scale.}
\end{figure}


For long distance transition probability, the computation is based
on the Chapman-Kolmogorov equation \cite{Paul}. The transition
probability density of a particle going from ${\bf p}$ to ${\bf q}$
is the total sum of the probabilities of going from ${\bf p}$ to
${\bf q}$ through all possible intermediate points ${\bf x} \in
\mathbb{R}^n$: \bqn P_t({\bf p}, {\bf q}) = \int_{\mathbb{R}^n}
P_s({\bf p},{\bf x})P_{t-s}({\bf x}, {\bf q}) \; d{\bf
x}\label{eq:kolmogorov}\eqn for any $0< s < t$. The equation still
hold in the case when $s$ is either 0 or $t$, since in that case one
of the probability in the integral will become the Dirac-delta
function and in turn the integral collapse to the probability on the
left side. For ${\bf x} \in B_{\bf p}$, small neighborhood centered
around ${\bf p}$, $D$ can be considered constant and hence $P_t({\bf
p},{\bf q}) \doteq K_t({\bf p} -{\bf q})$. Therefore, the
Chapman-Kolmogorov equation (\ref{eq:kolmogorov}) can be written as
\bqn P_t({\bf p}, {\bf q}) \doteq \frac{\int_{B_{\bf p}} K_s({\bf
p}-{\bf x})
 P_{t-s}({\bf x}, {\bf q}) \; d{\bf x}}
 {\int_{B_{\bf p}} K_{s}({\bf p}-{\bf x}) \; d{\bf x}}.
 \label{eq:transition} \eqn
The denominator is introduced to correct for the truncation error.
Note that when $s \to 0$ or $B_p$ becomes large, the equation
(\ref{eq:transition}) converges to equation (\ref{eq:kolmogorov}).
Let us define a truncated Gaussian kernel $$\tilde K_{t}({\bf
p}-{\bf x}) = \frac{K_{t}({\bf p} - {\bf x}){\bf 1}_{B_{\bf p}}({\bf
x})} {\int_{B_{\bf p}} K_t({\bf p} - {\bf x}) \;d{\bf x}}$$ where
${\bf 1}_{B_{\bf p}}({\bf x})$ is the index function taking value 1
if ${\bf x} \in B_{\bf p}$ and 0 otherwise. Rewriting the
Chapman-Kolmogorov equation, for each fixed ${\bf p}$, we have
$$P_t({\bf p}, {\bf q}) \doteq \tilde K_{s}*P_{t-s}({\bf p},{\bf q}).$$
The convolution is with respect to the first argument ${\bf p}$. By
letting $t=k\Delta t$ and $s=\Delta t$, we have iterative
relationship: $P_{k\Delta t}({\bf p},{\bf q}) = \tilde K_{\Delta t}
* P_{(k-1)\Delta t}.$ Solving this we have $k$-fold convolution
\bqn P_{k\Delta t}({\bf p},{\bf q}) = \underbrace{\tilde K_{\Delta
t}
* \tilde K_{\Delta t} * \cdots * \tilde K_{\Delta t}}_{k \mbox{
times}} *P_{0}({\bf p},{\bf q}), \label{eq:transition2}\eqn where
$P_{0}({\bf p},{\bf q}) = \delta({\bf p} - {\bf q})$. For numerical
implementation, we take the Kroneker-delta function rather than the
Dirac-delta function. We will denoting the $k$-fold convolution as
$\tilde K_{\Delta t}^{(k)}$. Note that
$$\int_{\mathbb{R}^n} P_{k \Delta t}({\bf p},{\bf q}) \; d{\bf q}=1$$
for all $k$.

The computation for the long distance transition probability can be
viewed as smoothing the initial probability distribution $P_0({\bf
p} - {\bf q}) = \delta({\bf p} -{\bf q})$. In the numerical
implementation, the initial probability 1 at the seed will be
smoothed over all voxels.


Figure \ref{fig:transition} shows the transition probability from
one voxel to the neighboring voxels computed from the anisotropic
Gaussian kernel. In the left figure, the transition probability is
computed without constrain while the right figure shows the
transition probability computed within FA $> 0.4$. Figure
\ref{fig:trans.prob} shows the long distance transition probability.
The bottom figure shows the transition probability represented in
log-scale, i.e. $\rho=\log P_t({\bf p},{\bf q})$. $\Delta t = 0.1$
with $k=40, 80, 120, 150$ iterations are used.






\section{Anisotropic kernel smoothing}


%A kernel-based method for solving
%(\ref{eq:diffusion2}) can be found in \cite{Spira}.

Instead of smoothing the initial probability, the same approach can
be used to smooth noisy data in such a way that we smooth more along
the direction the white matter tracks and less across the tracks.
Consider observation $Y$ which will be assumed to be
$$Y({\bf x}) = f({\bf x}) + \epsilon({\bf x}),$$
where $f$ is the mean function and $\epsilon$ is a zero mean
Gaussian random field. We assume that $\epsilon$ is an anisotropic
field. We are interested in estimating $f$ via kernel smoothing. If
the signal is defined along the white matter tracks, it is necessary
to incorporate the directional information of the white fiber into
the shape of kernel. In computing the transitional probability, we
smoothed the initial probability distribution by the iterated kernel
smoothing framework. Instead smoothing probability, we will smooth
observation $Y$ along the white fiber tracks. We formulate {\em
anisotropic Gaussian kernel smoothing} as the convolution of
observation $Y$ and kernel $K_t$: \bqn \hat f({\bf x}) = K_t
* Y({\bf x}) = \int_{\mathbb{R}^n} K_t({\bf x}-{\bf y})Y({\bf y})\;d{\bf
y} \label{eq:kernel1}. \eqn It can be shown that the above
convolution is the unique solution to (\ref{eq:diffusion1}) with the
initial condition $C({\bf x},0) = Y({\bf x})$ after time $t$ for
fixed $D$.


The anisotropic kernel smoothing can be viewed as a local
regression. To see this, consider partial differential operator
$\mathcal{L} \psi =\nabla \cdot (D \nabla \psi)$ and its eigenvalue
problem $\mathcal{L} \psi = \lambda \psi$. Since $D$ is positive
definite, we have ordered positive eigenvalues $0 = \lambda_0 \leq
\lambda_1 \leq \lambda_2 \leq \cdots$ and corresponding
eigenfunctions $\psi_0,\psi_1,\cdots$. The eigenfunctions $\psi_j$
form an orthonormal basis of $L^2(\mathbb{R}^n)$, the $L^2$ space of
functions defined on $\mathbb{R}^n$. %Any function $ h \in
%L^2(\mathbb{R}^n)$ can be expanded as $h = \sum_{j=0}^{\infty}
%h_j^*\psi_j$ where $h_j^* = \int  h({\bf x})\psi_j ({\bf x}) \;d{\bf
%x}.$
%Consider the finite number of basis $\psi_0, \psi_1,\cdots,
%\psi_l$ which spans subspace $\mathcal{H}_l \subset
%L^2(\mathbb{R}^n)$.
Then it can be shown that
%the anisotropic kernel $K_t$ can be
%written as \bqn K_{t}({\bf x}-{\bf y})
%=\sum_{j=0}^{\infty}e^{-\lambda_j\sqrt{2t}} \psi_j({\bf
%x})\psi_j({\bf y})\label{eq:heatkernel}.\eqn

%\cite{berline}
%\cite{rosenberg}.


\begin{theorem}
$$K_t*Y({\bf x}) =\arg\min_{h \in L^2(\mathbb{R}^n)}\int K_{t}({\bf x} -{\bf y})\Big[Y({\bf y}) - h\Big]^2 \; d\mu({\bf
y}).$$\label{theorem:finite.polynomial}
\end{theorem}
By expanding the integral as a quadratic in $h$: $ h^2 - 2h K_{t}*Y
+K_{t}*Y^2$ and differentiating with respect to $h$, we obtains the
minimum when $h({\bf x}) = K_{t}*Y({\bf x})$. This theorem is also
true for the $l$-dimensional finite subspace $\mathcal{H}_l \subset
L^2(\mathbb{R}^n)$ that is spanned by the finite basis functions
$\psi_0,\psi_1,\cdots,\psi_l$.

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.18]{transition.prob2.eps}
\caption{\label{fig:trans.prob2} Transition probability from a
single point in the corpus callosum with $t = 0.1$ and $k=8,16, 24,
28$ iterations. The scale in the top images is $10^{-3}$. The bottom
images are in natural log scale.}
\end{figure}

The following theorem shows the total amount of measurements is
invariant under smoothing.
\begin{theorem}[Conservation of total signal]
$$\int_{\mathbb{R}^n} K_t^{(k)} *Y \;d{\bf x}= \int_{\mathbb{R}^n} Y
\;d{\bf x}.\label{theorem:conservation}$$
\end{theorem}
It can be easily seen by interchanging the order of integral
$$\int_{\mathbb{R}^n} K_t *Y \;d{\bf x} =\int_{\mathbb{R}^n}\int_{\mathbb{R}^n} K_t({\bf
x} - {\bf y})Y({\bf y}) \;d{\bf y} \;d{\bf x}$$ and noting $K_t$ is
a probability distribution. Then the statement follows inductively.
The bias of the heat kernel smoothing  $\mathbb{E} \hat f - f =
K_t*Y -f \to 0$ as $t \to 0$. The total bias always vanishes:
$\int_{\mathbb{R}^n} [\mathbb{E} \hat f ({\bf x}) -f({\bf x})] \;
d{\bf x} = \int_{\mathbb{R}^n} [K_t*f({\bf x}) -f({\bf x})] \; d{\bf
x} = 0$. This is the consequence of Theorem
\ref{theorem:conservation}.

We have applied our smoothing method on each component of diffusion
tensor $D=(d_{ij})$ for increasing the signal-to-noise ratio in the
FA measure. The following theorem guarantee the positive
definiteness is preserved under smoothing.


\begin{figure}[t]
%figure 3
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{FA-smooth.eps}
\caption{\label{fig:FA-smooth} Left: FA image. Middle: isotropic
Gaussian kernel smoothing with 12mm FWHM. Right: anisotropic
Gaussian kernel smoothing with the equivalent filter size showing
the preservation of directional information. The advantage of the
anisotropic smoothing is that it increase the signal-to-noise ratio
without sacrificing much of the directional characteristic of DTI.}
\end{figure}

\begin{theorem}[Conservation of positive-definiteness]\\ If $D=(d_{ij})$ is positive definite, $K_H*D =
(K_t*d_{ij})$ is also positive-definite.
\end{theorem}
Note that $D$ is positive-definite if  ${\bf x}'D{\bf x} = \sum_{ij}
x_i x_jd_{ij} > 0$ for all ${\bf x}=(x_1,\cdots,x_n) \in
\mathbb{R}^n$. Then ${\bf x}'(K_H*D){\bf x}= \sum_{ij} x_ix_j
K_t*d_{ij} = K_t*\sum_{ij} x_i x_j d_{ij} >0$ proving the statement.

%The important property of the kernel smoothing is \bqn |K_H*f({\bf
%x})| \leq |f({\bf x})|,\eqn

Figure \ref{fig:FA-smooth} shows the FA images before and after
smoothing. For comparison, we performed both isotropic and
anisotropic Gaussian kernel smoothing with the similar bandwidth.
The standard unit of isotropic smoothing bandwidth is the full width
at the half maximum (FWHM) of the Gaussian kernel. For anisotropic
smoothing, we have taken the average full width at the half maximum
along the three principal axis of the kernel ellipsoid.


%In the case $V_j=0$ while $|{\bf V}| \neq 0$ in (\ref{eq:vector}),
%the kernel diverges. To avoid the divergence in numerical
%implementation, we can introduce very small $\epsilon$ and use
%${\bf V}+\epsilon$ instead of ${\bf V}$ in computing kernel. In
%the numerical computation, we used $\epsilon = 10^{-10}$.

%Further the above kernel has a problem in the
%case of constant field ${\bf V}=(1,\cdots, 1)/\sqrt{n}$. In this
%particular case the kernel is isotropic. Although it is not a
%problem in almost all DTI images, a way to avoid this pathological
%example is to introduce the cross-product terms into the matric
%tensor. Metric like $g_{ij}=V_iV_j$ would give such structure. It
%In this case the metric structure is given by
%We wish to find the covariance
%matrix $HH'$ of
%the kernel that gives ${\bf V}$ as a eigenvector. Without the full
%set of eigenvalues and eigenvectors, the exact determination of
%the covariance matrix is not possible. So we estimate the
%components of $HH'$ by minimizing the sum of the errors, i.e. the
%least-squares estimation. Solving $HH'{\bf V}=\lambda{\bf V}$, we
%have an estimate $\widehat{HH'} = \lambda{\bf V}{\bf V}^{-}$,
%where ${\bf V}^{-}$ is the Moore-Penrose generalized inverse of
%${\bf V}$. It is trivial to show that for unit column vector ${\bf
%V}^{-}={\bf V}'$. However ${\bf V}{\bf V}^-$ singular. To get the
%positive definite symmetric estimate $\widehat{HH'}$, we add small
%epsilon in the diagonal term:

%$$(HH'-\epsilon I){\bf V}=  (\lambda -\epsilon){\bf V}$$
%for some $0< \epsilon < \lambda$. Then taking the Moore-Penrose inverse, we get
%$$\widehat{HH'} =\epsilon I + (\lambda-\epsilon){\bf V}{\bf V}^{-}.$$
%When $\epsilon$ is relatively small compared to $\lambda$, i.e. $\lambda - \epsilon \approx \lambda$ so
%\bqn \widehat{HH'} = \epsilon I + \lambda {\bf V}{\bf V}'. \label{eq:penrose}\eqn
%Note that this basically corresponds to matching $HH'$ to the Riemannian metric tensor
%of the form
%$$d\psi^2= (\epsilon\delta_{ij} + \lambda V_iV_j) dx_idx_j$$
%for small $\epsilon$ which makes the metric positive definite. The
%corresponding kernel would be approximately
%$$K_{\epsilon}({\bf x}) = \alpha \frac{(\epsilon\lambda)^{1/2}}{(2\pi \epsilon)^{n/2}} \exp \Big( \frac{{\bf x}'(I-{\bf v}{\bf v}'){\bf x}}{2\epsilon}\Big),$$
%where $\alpha$ is a normalizing constant to make
%$K_{\epsilon}({\bf x})$ a probability density [blinded]. It can be
%computed numerically at each voxel. $\epsilon$ can be used to
%control the amount of smoothing.








%\begin{figure}[t]
%%figure 5
%\centering
%\includegraphics[scale=0.17]{DTIinv.eps}
%\caption{\label{fig:maxeig} Top: the original diffusion
%tensors $d_{11},d_{22}$ and $d_{12}$. Bottom: the inverse of
%the diffusion tensors, $d^{11},d^{22}$ and $d^{12}$, which are
%need for computing the anisotropic kernel. Note that $d^{ij}$ have
%been smoothed via the Cholesky factorization and normalized by
%the trace of the inverse.}
%\end{figure}



%Suppose that $S_j = ({\bf x}_1^j,{\bf x}_2^j, \cdots, {\bf x}_{m_j}^j)$ is the $j$-th path connecting ${\bf p}$ and ${\bf q}$ in a DTI. Even thought there would be infinite number of paths connecting from ${\bf p}$ to ${\bf q}$, we will only consider paths that are
%Since we know the transition probability of going from voxel ${\bf x_i^j}$  to adjacent voxel ${\bf x_{i+1}^j$ to be $P({\bf x_i}^j, {\bf x_{i+1}^j}) = K_t({\bf x}_{i+1}^j -{\bf x}_i^j)$, the above equation can be written as
%\bqn P({\bf p}, {\bf q}) &=& \sum_{S_j} P({\bf p}, {\bf x}_1^j)P({\bf p}, {\bf x}_2^j)\cdots P({\bf x}_{m_j}^j, {\bf q})\\
%&=&  \sum_{S_j}  K_t({\bf p}-{\bf x}_1^j) K_t({\bf x}_1^j -{\bf x}_2^j) \cdots K_t({\bf x}_m^j -{\bf q}),
%\eqn
%where the summation is taken over all possible path that connects from ${\bf p}$ to ${\bf q}$. Although this formulation is computationally feasible, the numerical implementation would be messy. Basically we need to sum over all possible %paths


%\begin{figure}[t]
%figure
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.20]{logp.eps}
%\caption{\label{fig:logp} Log-transition probability at 20, 40,
%60, 100,180 and 200 iterations respectively showing the
%anisotropic propagation of the initial probability 1. $t=0.1$ has
%been used for smoothing with kernel type (\ref{eq:vector}).}
%\end{figure}



%In discretizing the integral version of kernel smoothing
%(\ref{eq:kernel2}), we assume that the center of voxels in DTI to
%be on the integer lattice $\mathbb{Z}^n$. The affine
%transformation that maps DTI to the lattice $\mathbb{Z}^n$ without
%changing the maximal connectivity measure is used to simplify the
%computation.
%We took $B_{\bf x}$ to be the collection of voxels incident to
%voxel ${\bf x}$ and ${\bf x}$ itself. There are $3^n$ voxels in
%$B_{\bf x}$. The smaller the bandwidth parameter $\Delta t$, the
%more concentrated the weighting is to these incident voxels. The
%smaller support $B_{\bf x}$ is used to speed up the computational
%time.
%Let's illustrate how the method works in the
%case of isotropic kernel. Anisotropic kernel case is similar
%except that the kernel changes from voxel to voxel. The 2D
%isotropic kernel when $H=I$ for $\Delta t=0.1$ is given by
%$$\left(\begin{array}{ccc}
%    0.0054 &   0.0653 &   0.0054\\
%    0.0653 &   0.7958 &   0.0653\\
%    0.0054 &   0.0653 &   0.0054\\
%\end{array}\right).$$
%Contributions outside 9 voxels are negligible (less than $0.0001$)
%so we use $\Delta t \leq 0.1$ and $9$ voxel neighborhood for our
%computation. If a more accurate connection probability is desired,
%one may tempted to use more neighborhood; however, that is not
%necessary if $\Delta t$ is very small.
%Note that $\sum_{{\bf x}_j
%\in B_{\bf x}} K_{\Delta t}({\bf x}_j) = 1.0785$ while $\sum_{{\bf
%x}_j \in \mathbb{Z}^2} K_{\Delta t}({\bf x}_j) = 1.0787$, about a
%$0.0002\%$ difference. So even after 500 iterations using the
%above truncated kernel, there will be only $0.1\%$ difference to
%true value in the discrete lattice. Since the dimension of DTI is
%less than $500^3$, the probability of connection can be computed
%in less than 500 iterations. Normalizing the above kernel, we get
%$$\tilde K_{\Delta t}=\left(\begin{array}{ccc}
%    0.0050  &  0.0606  &  0.0050\\
%    0.0606  &  0.7378  &  0.0606\\
%    0.0050  &  0.0606  &  0.0050\\
%\end{array}\right).$$




\section{Conclusions and Discussions}
G\"{o}ssl {\em et al.} \cite{Gossl} applied an isotropic Gaussian
kernel smoothing with small filter size to get the smooth,
continuous representation of DTI data using a linear state space
model. However, the isotropic smoothing tends to smooth out the
directional characteristic of the DTI data. Motivated by this, we
have introduced a new method for smoothing data along the white
fiber tracks to get the continuous and smooth representation of the
data while preserving the directional characteristic of DTI by
performing spatially adaptive iterative anisotropic kernel
smoothing. The kernel is constructed using the concept of the
transitional probability of the diffusion process that directly
corresponds to the water molecule diffusion in DTI.


Compared to the previous approaches of a Monte-Carlo simulation
\cite{Koch} or solving a diffusion equation
\cite{Conturo,Parker:2000,Weinstein}, our kernel method is simpler
to implement numerically and statistical properties can be
investigated easily. To speed up the kernel smoothing, one may use
the decomposition scheme of Geusebroek \cite{Geusebroek}.

An additional advantage of our smoothing formulation is that it can
be interpreted differential geometrically. Consider vector field
${\bf V} =(V_1,\cdots,V_n)'$ that is the principal eigenvector of
$D$. Suppose we wish to smooth data along the field in such a way
that we smooth more along the direction of the field. Let the flow
${\bf x}=\psi(t)$ corresponding to the field is given by
$$\frac{d \psi}{d t} = {\bf V} \circ \psi(t).$$
This system of ordinary differential equations give a family of
integral curves whose tangent vector is ${\bf V}$. The line element
is
$$d\psi^2 = V_1^2\; dx_1^2 + \cdots + V_n^2 dx_n^2$$
and the Riemannian metric tensor is given by $G=(g_{ij}),
g_{ij}=V_{i}^2\delta_{ij}$. By matching the bandwidth matrix of a
general anisotropic kernel to the Riemannian metric tensor
proportionally, we get

$$ K_t({\bf x}) = (4\pi t)^{-n/2}\prod_{j=1}^n \frac{1}{|V_j|}\exp
\Big(-\frac{x_j^2}{4tV_j^2}\Big).\label{eq:vector}$$

By introducing the scaling parameter $t$, we left a room for
adjusting the amount of smoothing. The advantage of using only the
principal eigenvectors would be the simplicity of the implementation
while the drawback is that the method does not completely utilize
the full diffusion tensor $D$. To use $D$ directly, we note that $D$
gives a natural Riemannian metric tensor in the diffusion process,
i.e. $G=D$ \cite{DeLara}. This results in the same anisotropic
kernel we have introduced in computing the transition probability. A
similar Riemannian metric tensor approach in connection with DTI is
given in \cite{Chefdhotel,ODonnell}. In particular \cite{ODonnell}
matched the diffusion tensor $D$ to the inverse of the metric tensor
$G$ and applied it in computing the tensor-warped distances in the
white fibers.


%\begin{figure}
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.45]{DTIyy.eps}
%\caption{\label{fig:DTIyy} Isotropic kernel simulation with $H=I$, $\Delta t = 0.1$ for different number of iterations.}
%\end{figure}


%Authors wish to thank Keith Worsley of McGill University and John Ashburner of Functional Imaging Laboratory, London for useful
%discussions about simulating diffusion tensors. Also thanks is to Sigurd Angenent of the department of mathematics for the discussion of the Riemannian metric tensor that has the cross-product terms.


\begin{thebibliography}{9}
\small  % Use 9 point text.


%\bibitem{Alexander} Alexander, A.L., Hasan, H.M., Lazar, M.,
%Tsuruda, J.S., Parker, D.L. Analysis of Partial Volume Effects in
%Diffusion-tensor MRI. {\em Magn. Res. Med} {\bf 45}:770-780. 2001.

%\bibitem{Anderson} Anderson, T.W. {\em An introduction to multivariate
%statistical analysis}, Wiley, New York. 1984.

%\bibitem{Ashburner} Ashburner, J. and Friston, K.J. Voxel-Based
%Morphometry: The Methods. {\em NeuroImage} {\bf 11}:805-821. 2000.

\bibitem{Basser:96} Basser, P.J., and Pierpaoli, C.  Microstructural and
physiological features of tissues elucidated by
quantitative-diffusion-tensor MRI. {\em J. Magn. Reson. B} {\bf
111:}209-219. 1996.

\bibitem{Basser:00} Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J.,
Aldroubi, A. In Vivo Tractography using DT-MRI data. {\em Magn.
Res. Med.} {\bf 44}:625-632. 2000.

\bibitem{Batchelor} Batchelor, P.G., Hill, D.L.G., Atkinson, D., Calamanten,
F., Atkinson, D. Study of Connectivity in the Brain Using the Full
Diffusion Tensor from MRI, IPMI 2001, {\em Lecture Notes in
Computer Science} {\bf 2082}:121-133. 2001.

%\bibitem{Betounes} Betounes, D. {\em Partial Differential Equations for
%computational Science.}
% Springer-Verlag New York, Inc. 1998.

\bibitem{Chefdhotel} Chefd'hotel, C., Tschumperl\'{e}, D., Deriche, R., Faugeras,
O. Constrained flows of matrix-valued functions: Application to
diffusion tensor regularization. {\em Lecture Notes in Computer
Science} {\bf 2350}:251-265, Springer-Verlag, Berlin. 2002.

\bibitem{Conturo} Conturo, T.E., Lori, N.F., Cull, T.S., Akbudak, E.,
Snyder, A.Z., Shimony, J.S., McKinstry, R.C., Burton, H., Raichle,
M.E. Tracking Neuronal Fiber Pathways in the Living Human Brain,
{\em Proc. Natl. Acad. Sci.} {\bf 96}:10422-10427. 1999.

\bibitem{DeLara} De Lara, M.C. Geometric and Symmetric Properties of a
Nondegenerate Diffusion Process, {\em Annals of Probability}. {\bf
23}:1557-1604. 1995

%\bibitem{Dougherty} Dougherty, E.R. {\em Random processes for image and signal processing}. IEEE Press.
%1998.

\bibitem{Hagmann} Hagmann, P., Thiran, J.-P., Vandergheynst, P., Clarke, S.,
Meuli, R. Statistical fiber tracking on DT-MRI data as a potential
tool for morphological brain studies. {\em ISMRM Workshop on
Diffusion MRI: Biophysical Issues}. 2000.

%\bibitem{Harville} Harville, D.A. {\em Matrix algebra from a
%statistician's perspective}, Springer-Verlag New York. 1997.

\bibitem{Gossl}  G\"{o}ssl, C., Fahrmeir, L., P\"{u}tz, B., Auer, L.M., and
Auer, D.P. Fiber tracking from DTI using linear state space
models: delectability of the pyramidal tract. {\em NeuroImage}
{\bf 16:}378-388. 2002.

\bibitem{Geusebroek} Geusebroek, J.-M., Smeulders, A. W. M.,
van de Weijer, J. Fast Anisotropic Gauss Filtering, {\em Lecture
Notes in Computer Science} {\bf 2350}:99-112, Springer-Verlag,
Berlin. 2000.

%\bibitem{Kimmel} Kimmel, R., Sochen, N., Orientation diffusion or
%how to comb a porcupine?

\bibitem{kang.2004} Kang, N. and Zhang, J. and Carlson, E.S., and
Gembris, D. White matter fiber tractography via anisotropic
diffusion simulation in the human brain. Technical report. 410-04.
Department of Computer Science. University of Kentucky. Lexington.

\bibitem{Koch} Koch, M.A., Norris, D.G., and Hund-Georgiadis, M. An
investigation of functional and anatomical connectivity using
magnetic resonance imaging. {\em NeuroImage} {\bf 16}:241-250.
2002.

%\bibitem{Lazar} Lazar, M., Weinstein, D.M., Tsuruda, J.S., Hasan, K.M., Arfanakis,
%K., Meyerand, M.E., Badie, B., Rowley, H.A., Haughton, V., Field,
%A., Alexander. L. White Matter Tractography using Diffusion Tensor
%Deflection. {\em Human Brain Mapping} {\bf 18}:306-321. 2003

\bibitem{ODonnell} O'Donnell, L., Haker, S., Westin, C.-F., New Approaches to
Estimation of White Matter Connectivity in Diffusion Tensor MRI:
Elliptic PDEs and Geodesics in a Tensor-Warped Space. Fifth
International Conference on Medical Image Computing and
Computer-Assisted Intervention - MICCAI 02. 2002.


\bibitem{Parker:2000} Parker, G.J.M., Schnabel, J.A., Symms, M.R.,
Werring, D.J, Barker, G.J. Nonlinear Smoothing for Reduction of
Systematic and Random Errors in Diffusion Tensor Imaging. {\em
Journal of magnetic Resoance Imaging.} {\bf 11}:702-710. 2000.

\bibitem{Paul} Paul, W., Baschnagel, J., {\em Stochastic Processes from
Physics to Finance}. Springer-Verlag Berlin. 2000.

%\bibitem{Perona}
%Perona, P., Malik, J. \newblock Scale-space and Edge Detection
%Using Anisotropic Diffusion. \newblock {\em IEEE Trans. Pattern
%Analysis and Machine Intelligence}, 12:629-639, 1990.

\bibitem{Mori} Mori, S.,
Crain, B., Chacko, V.P., van Zijl, P.C.M. Three Dimensional
Tracking of Axonal Projections in the Brain by Magnetic Resonance
Imaging. {\em Ann. Neurol.} {\bf 45}:265-269. 1999.

\bibitem{Stevens} Stevens, C.F.
{\em The six core theories of modern physics}, The MIT Press. 1995.

\bibitem{Spira} Spira, A., Kimmel, R., Sochen, N.
Efficient Beltrami Flow using a Short Time Kernel. {\em Lecture
Notes in Computer Science.} {\bf 2695}: 511-522. 2003.

\bibitem{Tang} Tang, B., Sapiro, G., Caselles, V. Direction Diffusion. ICCV
1245-1252. 1999.

\bibitem{Tench} Tench, C.R., Morgan, P.S., Blumhardt, L.D., Constantinescu,
C. Improved White Matter Fiber Tracking Using Stochastic
Labelling, {\em Magn. Res. Med.} {\bf 48}:677-683. 2002.

\bibitem{Wang} Wang, Z., Vemuri, B.C., Chen Y, Mareci, T., Simultaneous
Smoothing and Estimation of the Tensor Field from Diffusion Tensor
MRI. CVPR. 2003.

\bibitem{Weinstein} Weinstein, D.M., Kindlmann, G.,Lundberg, E.
 Tensorlines: Advection-Diffusion Based Propagation Through Diffusion
 Tensor Fields. {\em IEEE Visualization} 249-253. 1999.\\
 %http://www.sci.utah.edu/publications/vis99$\_$tensorlines-paper.pdf

%\bibitem{Worsley} Worsley, K.J., Marrett, S., Neelin, P., Vandal, A.C.,
%Friston, K.J. and Evans, A.C. A unified statistical appraoch for
%determining significant signals in images of cerebral activation.
%{\em Human Brain Mapping} {\bf 4}:58-73. 1996.

\end{thebibliography}
\end{document}
