\documentclass[11pt,twocolumn]{article}
\usepackage{times,amsmath,amsfonts}
\usepackage[dvips]{graphicx,graphics}

\setlength{\textheight}{10in} \setlength{\textwidth}{7.1in}
\setlength{\columnsep}{0.3in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.3in} \setlength{\headsep}{0in}
\setlength{\oddsidemargin}{-0.6in}  % Centers text.

% Add the period after section numbers.  Adjust spacing.
%\newcommand{\Section}[1]{\vspace{-5pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
%\newcommand{\SubSection}[1]{\vspace{-2pt}\subsection{\hskip -1em.~~#1}
%        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray*}}
\newcommand{\eqn}{\end{eqnarray*}}
\newcommand{\bq}{\begin{eqnarray}}
\newcommand{\eq}{\end{eqnarray}}

\begin{document}

\title{Stat 992: Lecture 34\\
Anisotropic kernel smoothing I.}
\author{Moo K. Chung \tt{mchung@stat.wisc.edu}}
\maketitle \thispagestyle{empty}

%\begin{figure}
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.8]{homework05-1.eps}
%\end{figure}

\begin{enumerate}
\item {\em Motivation.} An image may consists of two different
classes and each class may have different variability so that when
we want to smooth image we may want to apply different bandwidth
to different classes. This requires smoothing along the class
boundaries and the smoothed data should not cross over the
boundary. Data may be defined along vector fields so that we need
to smooth data along the tensor fields (Figure 1).

\item {\em Diffusion tensor data.} Anisotropic kernel smoothing is
best illustrated by doing a specific example. Let
$D(x)=(d_{ij}(x)), x \in \mathbb{R}^n$ be raw diffusion tensor
images given in Figure 1. For fixed $x$, $d_{ij}(x)$ is a non
negative definite tensor which can be checked by computing the
determinant and trace. For simplicity we normalize $D$ by
$D/\max_x D(x)$. In this way two positive eigenvalues of $D$ are
constrained so that $1 \geq \lambda_1(x) \geq \lambda_2(x) \geq
0$. After downloading data from
http://www.stat.wisc.edu/$\sim$mchung/\\
teaching/datatensor/DTI.mat run following codes
\begin{verbatim}
load DTI
vecxx=vec(xx); vecxy=vec(xy);
vecyy=vec(yy);
gmax=max([vecxx;vecxy ;vecyy]);
xx=xx/gmax; xy=xy/gmax; yy=yy/gmax;
imagesc([xx xy yy]);colorbar;
\end{verbatim}

\item Consider principal eigenvector ${\bf V} =(V_1,\cdots,V_n)'$
of $D$. For simplicity, we will only consider smoothing along
${\bf V}$ for now. Running the following code, we generate ${\bf
V}$ and $\lambda_1$ (Figure 2).
\begin{verbatim}
xdim=170;ydim=203;
for i=1:xdim*ydim
    A=[vecxx(i) vecxy(i)
       vecxy(i) vecyy(i)];
    [v,d]=eig(A);
    V(i,:)=v(:,2)';
    lambda1(i)=d(2,2);
end;
\end{verbatim}
The eigenvectors are normalized so that $\|{\bf V}\|=1$. For
simplicity, we may let ${\bf V} \leftarrow \lambda_1{\bf V}$. The
stream line ${\bf x}=\psi(t)$ corresponding to the vector field is
given by
$$\frac{d \psi}{d t} = {\bf V} \circ \psi(t) + \epsilon(t)$$
where $\epsilon$ is random noise vector. This ordinary
differential equation (ODE) gives a family of integral curves
whose tangent vector is ${\bf V}$ (Betounes, 1998). The line
element is
$$d\psi^2 = V_1^2\; dx_1^2 + \cdots + V_n^2 dx_n^2.$$
which gives Riemannian metric tensor $g_{ij}=V_{i}^2\delta_{ij}$.
If we smooth some data with anisotropic Gaussian kernel $K_{H}$
such that we smooth more along the larger metric distance, we may
let $HH' = 2tG, G=(g_{ij}),$ i.e.
$$HH'= 2t
\left(\begin{array}{cccc}
  V_1^2 & 0 & \cdots & 0 \\
  0 & V_2^2 & \cdots & 0\\
  \cdots & \cdots & \cdots& \cdots\\
  0 & 0 & \cdots & V_n^3\\
\end{array}
\right).$$
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{ima34-1.eps}
\caption{Diffusion tensor images $D_{xx},D_{xy}, D_{yy}$.}
\end{figure}

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima34-4.eps}
\includegraphics[scale=0.25]{ima34-3.eps}
\caption{Left: Principal eigenvalues $\lambda_1$, Right: Principal
eigenvectors ${\bf V}$ displayed using {\tt quiver} command in
{\tt MATLAB}.}
%\caption{\label{fig:dti_quiver} Left: normalized FA-map. Middle:
%isotropic Gaussian kernel smoothing with 12mm FWHM. Right:
%anisotropic Gaussian kernel smoothing with generalized 12mm FWHM.}
\end{figure}

By introducing the scaling parameter $t$, we left a room for
adjusting the amount of smoothing possibly by an optimal
bandwidth. For the above choice of the covariance matrix,
anisotropic kernel is given by \bqn K_t({\bf x}) = \frac{1}{(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}\eqn In the case
$V_j=0$ while $|{\bf V}| \neq 0$ in (\ref{eq:vector}), the kernel
is undefined. To avoid the divergence in numerical implementation,
we can introduce very small $e$ and let $HH'= 2tG +eI$ instead of
${\bf V}$ in computing kernel.
\begin{verbatim}
%kernel computation
[dx,dy]=meshgrid(-1:1)
K=inline('exp(-(x.^2+y.^2)/(4*t))
/sum(sum(exp(-(x.^2+y.^2)/(4*t))))')
e=0.0001;
ker=zeros(ksize,ksize);
for i=2:(xdim-1)
    for j=2:(ydim-1)
      wx=  Vx(i-1:i+1,j-1:j+1)+e;
      wy=  Vy(i-1:i+1,j-1:j+1)+e;
      weight = K(0.5,dx./wx,dy./wy);
      kernel(:,:,i,j)=weight;
  end;
end;
\end{verbatim}
\begin{verbatim}
%iterated smoothing
f=lambda1; g=zeros(xdim,ydim);
for k=1:20
  k
  for i=2:(xdim-1)
    for j=2:(ydim-1)
      temp=f((i-1):(i+1),(j-1):(j+1));
      g(i,j)=sum(sum(temp.*kernel
             (:,:,i,j)));
    end;
  end;
  f=g;
end;
\end{verbatim}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima34-5.eps}
\includegraphics[scale=0.25]{ima34-7.eps}
\includegraphics[scale=0.25]{ima34-6.eps}
\includegraphics[scale=0.25]{ima34-8.eps}
\caption{Top left: 9-pixel neighboring kernel, Top right:
Spatially varying filter size. Filter size is large where the
directional information is strongest. Smoothing of the principal
eigenvalues. Bottom left: anisotropic smoothing with $m=20,t=0.5$
which gives maximum filter size of 3.1. Bottom right: Isotropic
kernel smoothing with the uniform filter size of 3.1.}
%\caption{\label{fig:dti_quiver} Left: normalized FA-map. Middle:
%isotropic Gaussian kernel smoothing with 12mm FWHM. Right:
%anisotropic Gaussian kernel smoothing with generalized 12mm FWHM.}
\end{figure}
 \item {\em Discussion.} This is a spatially adaptive filter
which has the property of reducing noise while preserving edge
information. Because the bandwidth changes spatially, in order to
be able to compare the result with isotropic kernel, we need to
come up with the amount of smoothing in average sense. Note that
at each position ${\bf x}$, the amount of smoothing is $2tV_j^2$
along the $j$-th coordinate. The shape of kernel is an ellipsoid
$\frac{x_1^2}{2tV_1^2} + \cdots + \frac{x_n^2}{2tV_n^2}=1.$ The
mean distance from the origin to the surface of the ellipsoid is
$\sigma({\bf x})=\sqrt{2t}\sum_{j=1}^n |V_j|/n$. Then we perform
iterated kernel smoothing with $m$ iterations, the total amount of
smoothing would be $\sigma_{total}({\bf x}) = \sqrt{m}\sigma({\bf
x})$. Figure 3 shows spatially varying filter size. The advantage
of using only the principal eigenvectors would be the simplicity
of the implementation while the drawback is that the computation
of the principal eigenvectors takes a fair amount of time in
$\tt{MATLAB}$ as well as the fact that the method does not
completely utilize full diffusion tensor information. 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. A way to avoid this is to introduce the cross-product
terms into the matric tensor. Metric like $g_{ij}=V_iV_j$ would
give such structure. However such metric tensor is not positive
definite.
\end{enumerate}
\end{document}

Consider following random field model on manifolds $Y_i(x) =
\mu_i(x) +\epsilon(x), x \in \partial \Omega$ for subject $j$. We
assume i.i.d. zero mean Gaussian random field $\epsilon_j$ with
isotropic covariance $R(x,y) = \sigma^2f(d(x,y))$ where $d$ is a
geodesic distance and $f$ some function.

Suppose there are other measurement $Z_j \sim N(\mu_Z
,\sigma_Z^2)$. The sample correlation between two measurements
$Y_j$ and $Z_j$ at fixed point $x$ is given by
$$\rho(x) = \frac{\mbox{Cov}(Y(x),Z)}{$$


\item {\em Covariance function estimation.} On manifolds $\partial
\Omega$, we have model
$$Y(x) = \mu(x) + \epsilon(x)$$
where $\epsilon$ is unknown random field. We may estimate $\mu$ by
kernel smoothing so that $\widehat{\mu} = K_{\sigma}*Y$. Now we
estimate the noise with the residual $\widehat{\epsilon} = Y -
K_{\sigma}*Y$.