\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 40\\
Image Registration 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}
This method is still under formulation. :)
\begin{enumerate}
\item {\em Function matching by iterative shifting.} We present a
curve matching methodology that can be extended in higher
dimension. Let $f(x)$ and $g(x)$ be two functions from
$\mathbb{R}$ to $\mathbb{R}$. Suppose that the two functions are
similar in shape such that
$$f(x) \approx g(x+u(x))$$
where $u$ is the displacement function. The reasonable stochastic
model would be
$$f(x) = g(x+u(x)) + \epsilon_0(x)$$
where $\epsilon_0$ is zero mean Gaussian fields. Linearize
$g(x+u(x))$ via Taylor expansion:
$$g(x+u(x))= g(x) + g'(x)u(x) + O(u^2(x))$$
We assume small displacement $u(x)$ otherwise algorithm based on
linear approximation may fail. We have
$$ f(x) - g(x+u(x))\approx f(x) - g(x) - g'(x)u(x) = \epsilon(x)$$
where $\epsilon$ is the stochastic noise and error of linear fit
combined together. At each point $x$, $u$ is estimated by \bq
u_0(x) = \frac{g(x)-f(x)}{g'(x)}.\eq To stabilize the algorithm,
we use $K_{\sigma}'(x) * g(x)$ in the denominator so that
$$ u_0(x) = \frac{K_{\sigma_0}*(g(x)-f(x))}{K_{\sigma_0}'(x) * g(x)}.$$
So our first estimation is $u_0(x)=
\frac{K_{\sigma}*(g(x)-f(x))}{K_{\sigma_0}'*g(x)}.$ Based on this
initial guess, we shift function $g(x)$ to $g(u(x)+u_0(x))$ and
perform the above procedure again and get the new estimate
$$ u_1(x) = \frac{K_{\sigma_1}*(g(x+u_0(x)) - f(x))}{K_{\sigma_1}'*g(x+u_0(x))}.$$
The total amount of displacement is $u_0 + u_1$. Following this
idea, we have iterative shifting algorithm
$$u_{n+1}(x)= \frac{K_{\sigma_n}*(g(x+u_0+u_1+\cdots u_n)-f(x))}{K_{\sigma_n}'*g(x+u_0+u_1+\cdots u_n)}.$$
We decrease the scale at each iteration $$\sigma_0 > \sigma_1 >
\sigma_2 \cdots$$ to have a coarse-to-fine approach. So the final
estimation is
$$\hat u = u_0 + u_1 + u_2 +\cdots.$$
In order for the algorithm to work, we need to show the
convergence of the
 iteration. If we have $|u_{n+1}| \leq c|u_n|$, $0<c<1$, the algorithm
 will converge due to contraction.

\item {\em Improvement via least squares estimation.} The
estimation (1) seems nice except it may not be very stable since
it is based on a single data point. A better estimation would be
to use a least squares estimation of $u$ requires a couple of
points around $x$. Suppose we use three points $x_{-}, x, x_{+}$
and set up a system of linear equations
$$ f(x_i) - g(x_i) - g'(x_i)u(x_i) = \epsilon(x_i).$$
We minimize the sum of squared residual of the 1st order Taylor
expansion:
$$u_0 = [g'(x)]^{-}[f(x) - g(x)]$$
where $[g'(x)] =(g'(x_), g'(x), g'(x_+))'$ denotes a vector and
similarly for $[f(x) - g(x)]$ and $^-$ is the generalized
Moore-Penrose inverse. It is trivial to show that for unit column
vector $V$, $V^-=V'$. Hence $[g'(x)]^- = [g'(x)]^T/\|g'(x)\|^2$.
This gives an weighted average estimation of the displacement.
Again the derivative $g'$ may not be well defined so we take the
derivative of kernel smoothing. So we estimate
$$u_0(x) = [K_{\sigma_0}*g'(x)]^-[K_{\sigma_0}(f(x) - g(x))].$$
Note that this is a weighted least squares estimation (0-th order
weighted polynomial regression with Gaussian weight).

Similar coarse-to-fine {\em multiscale} approaches have been used
in many image registration algorithms. It basically saves a lot of
computational time. Given data $Y$, multiscale data can be defined
as the collection of data
$$\{K_{\sigma}*Y:\sigma \in \mathbb{R}^+ \}.$$ Since the scale
$\sigma$ can vary continually, the multiscale representation of
image can contain global coarse-to-fine geometric features. For
instance $K_0*Y=Y$ is the original unsmoothed data while
$K_{\infty}*Y =\bar Y$ (assuming data is centered) is the mean
image. The advantages of the multiscale data representation are
the differentiability, compactness and reduction on computational
time.


\item {\em Derivative matching}. It is also possible to formulate
curve matching by matching higher order derivatives. It may give
more smooth displacement. This idea has been implemented by Thomas
Hoffmann via the dynamic time warping algorithm (Wang, K. and
Gaser, Y. Alignment of Curves by Dynamic Time Warping. Annals of
Statistics, 25:1251-1276, 1997).

Following the previous method for matching the curve functions, $$
f'(x)=g'(x) + g''(x)u(t)+ \epsilon(t)$$ Afterwards it follows
similarly.

The similarity measure we have used is based on function
difference. In the case of higher dimensional image data, it is
intensity difference. If we use more complex similarity measure
like mutual information, the problem becomes nonlinear
optimization. Other people used basis function expansion approach
such as polynomial basis: $u(x) = \sum_{i=1}^p \beta_i x^i$ (AIR)
or cosine basis (SPM'99).


\end{enumerate}
\end{document}


\item {\em Jacobian fields.} For a given displacement vector field
$U(t)$ from $O_1$ onto $O_2$, the actual deformation in the fixed
coordinate is given by
$$t \to t + U(t).$$
Then we are interested in local characteristic of the deformation
in the neighbourhood of a point $t \in O_1$. The Jacobian $J(t)=
det(I_N + \dot U)$ of this transformation measures the volume of
the unit cube at the transformed coordinate. The actual volume
change per unit volume is given by
$$K(t) = det(I_N + \dot U) - 1.$$
This quantity $K(t)$ is called {\it dilitation} in continuum
mechanics and its first order approximation is
$$K(t) \approx tr \dot U$$
for small deformation and we may simply let $K(t) = tr \dot U$
from now on. For a mean zero isotropic Gaussian random field
$U(t)$, $tr \dot U$ is again a mean zero isotropic Gaussian random
field with a different covariance function.



 \item{\em Multiscale data.}



\item {\em Optimal bandwidth.} On manifolds $\partial \Omega$, we
propose linear additive model
$$Y(x) = \mu(x) + \epsilon(x)$$
where $\epsilon$ is a mean zero Gaussian random field and $Y$ is
centered data. We can always make data centered by substracting
the sample mean $\bar Y$. We estimate $\mu$ by the usual kernel
smoother so that $\widehat{\mu} = K_{\sigma}*Y$. Let us denote the
residual by $$r(x) = Y(x) - K_{\sigma}*Y(x).$$

The performance of kernel smoother is measuredby the mean
integrated squred error (MISE)
$$\mathbb{E} \int_{\partial \Omega} [Y(x) - K_{\sigma}*Y(x)]^2 \;dx$$



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





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