\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 27\\
Laplacian in regular grid.}
\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}
The lecture is based on {\bf Problem 27.}
\begin{enumerate}
\item {\em 4 neighbor Laplacian.} Consider a regular grid system
$\delta l\mathbb{Z}^2$. It was shown that the 4 neighbor Laplacian
weights are given as
$$(w_{ij}) = \left(%
\begin{array}{ccc}
  0 & 1 & 0 \\
  1 & -4 & 1 \\
  0 & 1 & 0 \\
\end{array}%
\right).$$ Consider continues function $f$ evaluated in the
neighborhood of $(x_0,y_0)$. The Laplacian is
$$\Delta f = \frac{1}{\delta l^2}\sum_{i,j=-1,0,1} w_{ij}f(x_i,y_i)$$
and the iterative scheme is

$$f(x_0,y_0) \leftarrow f(x_0,y_0) + \lambda \sum_{i,j}
w_{ij}f(x_i,y_i)$$ where $\lambda = \delta t/\delta l^2$ is the
diffusion coefficient. This can be written as
$$f(x_0,y_0) \leftarrow \sum_{i,j} W_{ij}f(x_i,y_i)$$
where
$$(W_{ij}) = \left(%
\begin{array}{ccc}
  0 & \lambda & 0 \\
  \lambda & 1-4\lambda & \lambda \\
  0 & \lambda & 0 \\
\end{array}%
\right).$$ Note that $\sum_{ij} W_{ij}=1$.



\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.24]{ima27-1.eps}
\includegraphics[scale=0.24]{ima27-2.eps}
\includegraphics[scale=0.24]{ima27-3.eps}
\includegraphics[scale=0.24]{ima27-4.eps}
\caption{Top left: original image, Top right: 4 neighbor Laplacian
of the image based on built in {\tt MATLAB} function. Bottom left:
4 neighbor Laplacian of the image based on {\tt Laplacian4}.
Bottom right: 8 neighbor Laplacian of the image based on {\tt
Laplacian8}.}\end{figure}


\item {\em 8 neighbor Laplacian.} Since $K_{\sqrt{2\delta t}}*f$
is the solution of the diffusion equation after time $\delta t$,
we let
$$f(x_0,y_0) + \lambda \sum_{i,j} w_{ij}f(x_i,y_i) \approx
K_{\sqrt{2\delta t}}*f(x_0,y_0).$$ Similar formulation has been
used to estimates the weights in terms of exponential functions in
Belkin, M. and Niyogi, P., Laplacian eigenmaps for dimensionality
reduction and data representation (2002).

Due to symmetry, we only need to evaluate the kernel at
$(0,0),(\delta l, 0), (\delta l, \delta l)$. $K(0,0) \propto
\frac{1}{\delta t}$, $K(\delta l,0) \propto \frac{1}{\delta
t}e^{-\lambda/4} \approx \frac{1}{\delta t}(1-\frac{\lambda}{4})$
and $K(\delta l, \delta l) \propto \frac{1}{\delta
t}e^{-\lambda/2} \approx \frac{1}{\delta t} (1-
\frac{\lambda}{2})$. Since $\delta l$ is small, the integral can
be estimated by the finite sum so that weights are proportional to
$1, 1-\lambda/4, 1- \lambda/2$. Now the kernel weights are given
by
$$(W_{ij}) = \frac{1}{9 -3\lambda}\left(%
\begin{array}{ccc}
  1-\frac{\lambda}{2} & 1-\frac{\lambda}{4} &   1-\frac{\lambda}{2} \\
   1-\frac{\lambda}{4} & 1 &   1-\frac{\lambda}{4}\\
   1-\frac{\lambda}{2} & 1-\frac{\lambda}{4} &   1-\frac{\lambda}{2} \\
\end{array}%
\right).$$ The weights are normalized so that $\sum_{ij} W_{ij} =
1$. Then the 8 neighbor Laplacian weight is
$$(w_{ij}) = \frac{1}{9 -3\lambda}\left(%
\begin{array}{ccc}
  1-\frac{\lambda}{2} & 1-\frac{\lambda}{4} &   1-\frac{\lambda}{2} \\
   1-\frac{\lambda}{4} & 3\lambda-8 &   1-\frac{\lambda}{4}\\
   1-\frac{\lambda}{2} & 1-\frac{\lambda}{4} &   1-\frac{\lambda}{2} \\
\end{array}\right)%
.$$ Since $\lambda$ is very small, we have 8 neighbor Laplacian
weights
$$(w_{ij}) = \frac{1}{9}\left(%
\begin{array}{ccc}
  1 & 1 & 1 \\
  1 & -8 & 1 \\
  1 & 1 & 1 \\
\end{array}%
\right).$$
\begin{verbatim}
load trees
imagesc(X);colormap('bone');colorbar
L = del2(X)%built-in Laplacian.

Laplacian4=[0 1 0
          1 -4 1
          0  1 0]
L4X=conv2(X,Laplacian4,'same');

Laplacian8=[1 1 1
          1 -8 1
          1  1 1]
L8X=conv2(X,Laplacian8,'same');
\end{verbatim}

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.24]{ima27-5.eps}
\includegraphics[scale=0.24]{ima27-6.eps}

\caption{Left: smoothed image. Right: Laplacian of smoothed image.
Comparing with the previous figures, it gives more crisp edges due
to the suppression of noise.}
\end{figure}

{\bf Problem 40.} Based on the above 4 neighbor liberalization of
the above Laplacian, find function $f$ that minimize
$$\sum_{ij} (Y(p_{ij}) - f(p_{ij}))^2 + \lambda \int_{\Omega}
\Delta f(p)f(p) \;d\mu(p).$$ The integral can be approximated as
$\sum_{ij} \Delta f(p_{ij})f(p_{ij}).$ $f$ gives a smoothing of
data $Y$ defined on $\Omega$. Apply it to the sagittal image
$\Omega$ and smooth the image intensity $Y(p_{ij})$ this way.
Hint: Read Belkin and Niyogi's Semi-Supervised Learning on
Riemannian Manifolds (2003). I think Tim Ramsay (2000) has done
something similar and called it Laplacian smoothing.

\item {\em Edge detection.} In the place where the image intensity
changes rapidly, the Laplacian will have extreme values. This is
how the Laplacian can be used to detect edges in images without
much efforts; however the Laplacian is not well defined in images
due to the discreetness of image intensity. To make it more
smooth, kernel smoothing is usually performed before computing the
Laplacian, i.e. $\Delta K_{\sigma}*Y$.

\begin{verbatim}
Gaussian=[0 1/2 0
          1/2 1 1/2
          0  1/2  0]/3
GX=X;
for i=1:100
  GX=conv2(GX,Gaussian,'same');
end;
GL4X=conv2(GX,Laplacian4,'same');
\end{verbatim}

An alternate much simpler way to detect edge may be to use local
variance map. For image intensity $Y$, let ${\bf Var}_d(Y)$ to be
the variance map of $Y$ where the sample variance is computed
locally in square of size $d \times d$. In the regions of
homogenous image intensity, intensity variance doesn't change much
while in the region of rapidly changing image intensity, ${\bf
Var}_d(Y)$ will have large values. So the local intensity variance
can be used to extract edge information. In practice, ${\bf
Var}_d(Y)$ is not a continuous function. To make it more smooth,
possibly for modeling the covariance function, we smooth with
Gaussian kernel before the variance computation, i.e. ${\bf
Var}_d(K_{\sigma}*Y)$.

\begin{verbatim}
for i=2:247
  for j=2:348
    SX(i,j)=std(reshape(X((i-1)...
        :(i+1),(j-1):(j+1)),9,1));
  end
end;
\end{verbatim}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.24]{ima27-8.eps}
\includegraphics[scale=0.24]{ima27-9.eps}

\caption{Left: local variance map. Right: local variance map of
smoothed image. It is the variance map of Figure 2 left image.}
\end{figure}
\end{enumerate}

{\bf Problem 41.} Estimate the covariance function $R(t,s)$ of the
smoothed tree image (Figure 2 Left). $z=R(t,s)$ gives a 2D surface
with respect to variables $t$ and $s$. Plot your covariance
function surface. You may need to use {\tt mesh} in {\tt MATLAB}.


\end{document}