\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 37\\
Anisotropic kernel smoothing II.}
\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 Image gradient.} In Lecture 34,
 we smoothed data $Y(x)$ via iterated kernel smoothing of form
$$K_{H}*K_{H}*\cdots*K_{H}*Y(x)$$
where $H$ is determined by matching it to the Riemannian metric
tensor $G$ of data, i.e. $H=2tG$. The parameter $t$ measures the
over all smoothness and by leaving it this way, we are leaving a
room for choosing it automatically as a optimal filtering problem.

In diffusion tensor imaging (DTI), the diffusion coefficient gives
the the natural metric. We generalize this idea to gray scale
images where the metric tensor is not given. The problem is
figuring out the most natural metric tensor of the image. Suppose
image $I$ can be modelled additively as
$$I(x) = \mu(x) + K_{M}*W(x)$$
where $W$ is white noise such that ${\bf Var} (K_M*W) =1$. Then it
is shown that the smoothness of image can be measured by the
variance of the image derivatives, i.e. ${\bf Var} (\partial I) =
(MM')^{-1}/2$. So this gives a motivation for using the image
gradient as the image metric that measures. Unfortunately the
image gradient is not well defined. To make it numerically well
defined, we differentiate after smoothing. This is the usual way
in computer vision area. Also see Fan's nonparametric book where
the derivative is estimated this way. Thomas Hoffmann estimated
higher order derivatives after spline smoothing. See
http://www.stat.wisc.edu/~mchung/papers/HBM2004\\
/HBM2004thomas.html. Yuefeng Lu is working on derivative
estimation problem using polynomial regression. So we estimate the
image gradient as
$$\widehat{\partial I} = \partial (K_{\sigma} *I) =
(\partial K_{\sigma}) *I.$$ Figure 2 left image shows the image
gradient estimated this way. Now if we take this as the image
metric and match it proportionally to the covariance matrix of
kernel, i.e. $HH'=2tG= 2t \widehat{\partial I}$, metric will be
larger near edges so the kernel smoothing will smooth more across
the edges which is not what we want to do.
\begin{verbatim}
file='dooly.jpg'; I=imread(file);
dooly=I(:,:,1);
[xdim, ydim]=size(dooly)
ndooly=dooly+normrnd(0,10,xdim,ydim);
sndooly=gaussblur(ndooly,9);
[Vx,Vy]=gradient(sndooly);
\end{verbatim}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.5]{ima37-1.eps}
\includegraphics[scale=0.5]{ima37-2.eps}
\caption{Left: a cartoon image, intensity scale is between 0 and
255. Middle: Gaussian white noise $N(0,10)$ added. Right:
Isotropic kernel smoothing. Noise has been removed inside the big
eyes but the boundary has been also blurred.}
\end{figure}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima37-3.eps}
\includegraphics[scale=0.25]{ima37-4.eps}
\caption{Left: Gradient fields $V=(\partial K_{\sigma})*I$. If we
use this as the metric, we will smooth more across the edges. We
need to come up with the vector fields that will not smooth near
edges. Right: fields $1/(1+|V_i|)$. These vector fields are not
smooth enough.}


\end{figure}
\item {\em Image metric.} What we need is a smaller metric near
the edges. We can construct it in by inverting the gradient but to
make it numerically stable, we let
$$G = \left(%
\begin{array}{cc}
  1+ |\widehat{\partial_{x_1} I}| & 0 \\
  0 & 1 + |\widehat{\partial_{x_2} I|}\\
\end{array}%
\right)^{-1}.$$ We have taken the absolute value of the gradient
to avoid singularity when $\widehat{\partial_{x_i} I} +1=0$. This
is fairly reasonable metric. Similar numerical technique has been
also used in Sethian's level set method. There are different kind
of metric to investigate. Hessian matrix of the image intensity,
i.e. $\partial_{x_i}\partial_{x_j} I$ or the Fisher information
matrix, local sample covariance matrix. Unfortunately the above
metric will not be smooth enough due to the inverse operator. So
we smooth out the vector fields so that our final metric would be
$$G = K_{\sigma}*\left(%
\begin{array}{cc}
  1+ |\widehat{\partial_{x_1} I}| & 0 \\
  0 & 1 + |\widehat{\partial_{x_2} I|}\\
\end{array}%
\right)^{-1}.$$

\begin{verbatim}
Vx=1./(1+abs(Vx));
Vy=1./(1+abs(Vy));
Vx=gaussblur(Vx,9);
Vy=gaussblur(Vy,9);
\end{verbatim}

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.5]{ima37-5.eps}
%\includegraphics[scale=0.25]{ima37-4.eps}
\caption{ Smoothed vector fields $K_{\sigma}*1/(1+|V_i|)$.}
\end{figure}

Then the covariance matrix of the kernel is $HH'= 2tG$ and we
perform iterated kernel smoothing with $K_H$. We performed 100
iterations with $t=0.5$ which gives the right images in Figure 4.
Notice how it clears noise but preserves all the geometric
features.

\begin{figure}[t] \centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.5]{ima37-6.eps}
%\includegraphics[scale=0.25]{ima37-4.eps}
\caption{Left: original image, middle: noisy image, right:
anisotropic smoothing of noisy image. See how it reduces noise
dramatically without blurring the geometric features.}
\end{figure}

\item {\em M-estimation.} There is  a slightly different appraoch
to edge preserving smoothing. Given additive model
$$I(x) = \mu_(x;\sigma) + \epsilon(x)$$ for image intensity, the $M$-estimator of
unknown parameter $\sigma$ is defined as a minimizer $$\hat \sigma
= \arg \min_{\sigma} \int \rho \big(I(x) - \mu(x;\sigma)\big)
\;dx$$ In the case $\rho(z) = z^2$, we have the usual least
squares estimation (LSE). The motivation for M-estimation is the
(LSE) is not robust enough in the presence of outliers. The usual
choice for the score function is that $\rho(z)/z^2 < \infty$ as $z
\to \infty$.

For edge detecting problem, zero-one score function $\rho(z) =
\mathbf{1}_{|z|>k}(z)$ has been used (Chen {\em et al.}, 2001,
CVPR). Chu {\em et al.} (1998. JASA) used $\rho(z) = -\exp
(-z^2/2\sigma^2)/\sigma$ and $\mu(x,\sigma) = \mu(x)$ which
downweights image intensity far from the estimate. Huber in Robust
statistics 1981 used $\rho(z) = z^2$ if $|z| \leq k$ and
$2k|z|+k^2$ if $|z| > k$.

The main problem is formulating kernel smoothing in the context of
the M-estimation problem which would give a statistical
justification for this adaptive kernel smoothing. This requires
figuring out the variational version of the kernel smoothing and
diffusion smoothing. See Cai's Weighted Nadaraya-Watson Regression
Estimation for reference and Fan's book. Consider the following
discrete minimization problem
$$\hat \beta(x) = \arg\min_{\beta_j} \sum_{i=1}^n
K_{\sigma}(x-x_i)\big[Y(x_i) -\sum_{j=0}^p \beta_j(x-x_i)^j\big]^2
$$ It can be shown that the Nadaraya-Watson estimator
$\hat \mu_{\sigma} = \widetilde{K_{\sigma}}*Y$ is the minimizer
when $p=0$. The continues version would be obviously
$$K_{\sigma}*Y(x) = \arg\min_{\beta_0} \int
K_{\sigma}(x-y)\big[Y(y) -\beta_0]^2 \;dy$$ This should be
trivially true. Then the anisotropic version should be also
correct.\\ {\bf Theorem.} $$K_{H}*Y(x) = \arg\min_{\beta} \int
K_{H}(x-y)\big[Y(y) -\beta]^2 \;dy.$$ Then it follows that we can
extend this result to manifolds so that we have the following most
general theorem.

\begin{figure}[t] \centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.5]{ima37-7.eps}
%\includegraphics[scale=0.25]{ima37-4.eps}
\caption{Left: original image, Right: noisy image}
\end{figure}

{\bf Theorem.} For $C^2$ manifolds $\partial \Omega$ with the
Laplace-Beltrami operator $\Delta$ with eigenvalue $\lambda_j$ and
eigenfunction $\psi_j$, the solution to diffusion equation
$$\partial_t f = \Delta f$$
with initial value $f(x,0)=Y(x)$ is the minimizer
$$K_{\sqrt{2t}}*Y(x) = \arg\min_{\beta} \int_{\partial \Omega}
K_{\sqrt{2t}}(x,y) \big[Y(y) - \beta\big]^2\; d\mu(y).$$ where
heat kernel is given by $$K_{\sigma}(x,y)
=\sum_{j=0}^{\infty}e^{-\lambda_j\sigma} \psi_j(x)\psi_j(y).$$
Hence the diffusion and kernel smoothing on manifold can be viewed
as the local polynomial regression. :)



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