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



\setlength{\textheight}{10.5in} \setlength{\textwidth}{7.3in}
\setlength{\columnsep}{0.3in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.7in} \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 06\\
Iterated kernel smoothing.}
\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 Iterated kernel smoothing.} From the property of a
diffusion smoother, we have \bq K_{m \sigma}
=\underbrace{K_{\sigma}*\cdots *K_{\sigma}}_{m^2 \mbox{
times}}.\eq Based on this identity, we implement kernel smoothing
numerically. Let $\sigma=0.4$ and apply iterated kernel smoothing.
With small window size $5 \times 5$, Gaussian kernel smoothing
with $\sigma =28$ will chop off the tail regions too much and this
is one of the motivation for doing iterated kernel smoothing.
\begin{verbatim}
w=K(0.4,dx,dy)/sum(sum(K(0.4,dx,dy)))
smooth_Y=Y;
for i=1:250
  smooth_Y=conv2(smooth_Y, w,'same');
  ISE(i)=sum(sum((mu-smooth_Y).^2));
end;
\end{verbatim}


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima06-0.eps}
\includegraphics[scale=0.25]{ima06-15.eps}
\includegraphics[scale=0.25]{ima06-2.eps}
%\includegraphics[scale=0.2]{ima06-3.eps}
\includegraphics[scale=0.25]{ima06-4.eps}
\includegraphics[scale=0.25]{ima06-5.eps}
%\includegraphics[scale=0.25]{ima06-6.eps}
\caption{Original unknown signal $\mu$ and kernel smoothing
estimations corresponding to $\sigma=0.8 (m=2),\sigma=4
(m=10),\sigma=28 (m=70)$. $\sigma=28$ is obviously the case of
over smoothing while $\sigma=0.8$ is under smoothing. The bottom
left image is ISE plotted over the number of iterations.
\label{fig:lec6-1}}
\end{figure}

{\bf Problem 12.} Come up with a similar identity like equation
(1) for anisotropic kernel smoothing with bandwidth matrix $H$.

\item {\em Integrated squared error.} When $\sigma=28$ with
$490$-iterations, Figure \ref{fig:lec6-1} shows obvious over
smoothing while $\sigma=0.8$ with $4$-iterations are under
smoothing. So there must be an optimal bandwidth that estimate the
signal best. Most automatic bandwidth selection methods are based
on i.i.d. error assumption so the results do not usually work in
imaging where image intensities have highly correlated geometric
structure.  A reasonable way to chose optimal bandwidth that
minimizes the fit of the estimation in $L_2$ norm, i.e.
$$\min_{\sigma} \int \big[\mu(x) - K_{\sigma}*Y(x)\big]^2 \; dx$$
which is called the integrated squared error (ISE). From our
simulation the optimal bandwidth is when $\sigma=\sqrt{18} \cdot
0.4 = 1.70$. Unfortunately this gives under smoothing. See J.S.
Simonoff's Smoothing Methods in Statistics (1996) for optimal
kernel smoothers. With i.i.d. error assumptions, you may get rid
of outliers easily and run statistical analysis again but in
imaging data we can not just delete outliers since it will destroy
the geometry of the image itself.

{\bf Problem 13.}\footnote{{\bf Project 4.} Investigation of CV in
2D image. Most CV results in the statistical literature that
assume i.i.d. error will not work here.} Check if cross-validation
(CV) method can produce a reasonable optimal bandwidth. For
example, out of $101 \times 101$ image data $Y$, you may decompose
images into training image $\tt{ Y_L=Y(1:2:101,1:2:101)}$ and
validation image $\tt{Y_R=Y(2:2:101,2:2:101)}$ and try to minimize
the integrated squared error in the validation image using the
model fit of the training image. You may try other partitioning
schemes. CV is a standard technique in computer vision and
statistics. See Ruppert {\em et al.} (2003) Semiparametric
Regression.


\item {\em Mesh data format.} As an application of iterated
smoothing approach, we will show iterated kernel smoothing for
surface data. Download all files
http://tezpur.keck.waisman.wisc.edu/$\sim$chung/teaching/
data/MNImesh/

\begin{verbatim}
[tri,coord,nbr,normal]=MNIgetmesh
                       ('gray.obj');
load thickness.vv
trisurf(tri,coord(1,:),coord(2,:),
         coord(3,:));

\end{verbatim}

$\tt{tri}$ gives the list of 3 indices that form a triangle.
$\tt{coord}$ gives the coordinates of vertices $\tt{nbr}$ gives
the list of adjacent neighbors. $\tt{normal}$ gives unit normal
vectors at each vertex. $\tt{trisurf}$ will display the surface
mesh in {\tt MATLAB}. Some graph theoretic concept is useful in
describing and implementing kernel smoothing on the mesh. The
degree of a vertex is the number of edges that connects to the
vertex. Read $\tt{MNIgetmesh.m}$ for understanding how one gets
it.

\begin{verbatim}
>> coord(:,1:4)
ans =
   49.9831   44.0817  -54.7131  -40.8825
   30.4634  -51.8008   32.6077  -58.8814
    7.4515   19.1203   17.1390   19.4011

>> tri(1:5,:)
ans =
           3       10243       10245
       10243        2563       10244
       10245       10243       10244
       10245       10244        2565
        2563       10246       10248

>> normal(:,1:4)
ans =
    0.4151   -0.1340   -0.9101   -0.2198
   -0.1868    0.7082    0.1541    0.9744
    0.8904   -0.6932    0.3847    0.0482

>> nbr(1:3,:)
ans =
       15473       15474       16484       28770       28772           0
       30322       30321       30324       31844       38497           0
       10243       10245       11827       19554       20082           0
\end{verbatim}

\item {\em Topological invariant.} We will only deal with closed
surface mesh that is topologically invariant to a sphere. One way
to check it is to compute the Euler characteristic $\chi(M)$of the
triangular mesh $M$ which is defined as
$$\chi(M) = V -E + F$$
where $V,E$ and $F$ are the number of vertices, edges and faces
respectively. In our data $V= 40962$ and $F= 81920$. For a mesh
topologically equivalent to a sphere, $\chi(M)=2$. If we count the
number of edges $E$ by counting the number of faces it is $3F$ but
we are counting them twice so $3F=2E$. So simplifying  the
relation we have \bq 2V - F = 4. \label{eq:EC}\eq Equation
(\ref{eq:EC}) gives a way to check the topology of mesh. A
particular note of interest is that from the {\em Gauss-Bonnet
theorem}, the integral of the Gaussian curvature over $M$ is
$2\pi\chi(M)$. So when we integrate the Gaussian curvature
$\kappa_G$ over the brain surface $M$, we should get
$$\int_M \kappa_G(x) \; dx = 4\pi.$$


\end{enumerate}
{\bf Homework 2.} Submit solutions to problems 4-12 by Feb 16
11:00am. You only need to do two thirds of the total problems I
pose in the course. For computational problems, please do not just
dump me codes. It is your responsibility that you document your
codes and proving the correctness of your algorithm.


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima06-8.eps}
\includegraphics[scale=0.2]{ima06-7.eps}
\caption{The left image is noisy thickness measure on the surface
of brain. The right image shows an image processing artifact that
will cause serious problem in statistical analysis later. For
example in the estimation of curvature \label{fig:lec6-2}}
\end{figure}

\end{document}