\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 38\\
Image Registration 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 Similarity measure.} Consider two geometric objects
$\Omega_1$ and $\Omega_2$ in images. We are interested in testing
if the shape of two objects $\Omega_1$ and $\Omega_2$ are
statistically different. This requires quantifying the shape
variations via the concept of curve and image registration and
template. Let $I_i$ be the image intensity that defines the
geometric object $\Omega_i$. In order to compare image intensities
across images, we need similarity metric that measures the
closeness of two images. One very simple similarity metric would
be the local sample intensity correlation $\rho$ defined as
$$\rho(I_1(x),I_2(x)) = \frac{{\bf Cov} (I_1(B_x), I_2(B_x))}{{\bf Var}^{1/2} I_1(B_x) {\bf
Var}^{1/2} I_2 (B_x)},$$ where $B_x$ is a small ball centered
around $x$. Other similarity measure such as the {\em mutual
information} have been used. We have computed sample intensity
correlation for both $5 \times 5$ and $10 \times 10$ windows
(Figure 3). Some regions show extremely lower correlation
indicating those regions are not matching anatomically. Then the
question is if we can increase the anatomical matching by
deforming image. This is the basic idea of intensity-based image
registration

\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{ima38-1.eps}
\caption{Sulcal patterns of brain at age 14 (left) and at age 19
(right). The global patterns are similar but there is a pattern
variation present. The question is if these pattern difference is
statistically significant.}
\end{figure}

\item {\em Image deformation.} Given image intensity $I_1(x)$, $x
\in \Omega_1$, the deformation of image via displacement field
$u=(u_1,\cdots,u_n)$ to be $I_(x+ u(x)), x \in \Omega_1$. Then we
will increase the matching of two images by increasing the total
similarity measure,
$$\int \rho^2(I_1(x+u(x)),I_2(x)) \;dx.$$
One may use $L_2$ norm such that the total similarity measure
would be
$$\int [I_1(x+u(x))-I_2(x)]^2\;dx.$$
Since planer curves can be characterized by the curvature
function, Thomas Hoffmann use the curvature difference as the
similarity measure in the curve registration problem, i.e.
$$\int [\kappa_1(x+u(x))-\kappa_2(x)]^2 \;dx.$$
However the minimization of the above similarity will result in
not very smooth displacement field. So we need to introduce a
measure of roughness. The usual measure of roughness is to use the
first derivative of the displacement
$$\|\partial u\| = \sum_{i,k} \Big|\frac{\partial u_k}
{\partial x_i}\Big|^2$$ or the second derivatives
$$\|\partial^2 u\| = \sum_{i,j,k} \Big|\frac{\partial^2 u_k}
{\partial x_i \partial x_j}\Big|^2.$$ Then we maximize the
similarity while minimizing the roughness, i.e.
$$\hat u=\arg\max_u \int \rho(I_1(x+u),I_2) \;dx - \lambda \int \|\partial^2 u \|^2 \;dx$$
where $\lambda >0$. Most of research in image registration have
been concentrated on choosing proper similarity and roughness
measures. See Brain Warping by Toga and Functional Data Analysis
by Ramsay and Silverman for an overview of image and curve
registration problem.


\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.5]{ima38-2.eps}
\caption{Brain cross-sections of two different subjects}
\end{figure}

\item {\em Shape inference.} Under the null assumption of no shape
difference between two geometric objects, the displacement fields
would be small. So it would be reasonable to model the
displacement field stochastically:
$$u(x) = \mu(x) + \Sigma^{1/2}(x)\epsilon(x).$$
where $\epsilon$ is the zero mean random vector field whose
components are i.i.d. and $\Sigma$ is the covariance matrix
accounting for the spatially varying correlated noise. See
\begin{verbatim}
http://www.math.mcgill.ca/chung/
deformation/ni_deformation.pdf
\end{verbatim}for detail. This model has been applied to the
detection of brain tissue growth in a group of 28 children.

The hypothesis of interest is
$$H_0: \mu(x) =0 \mbox{ for all } x$$
vs. $$H_1: \mu(x) \neq 0 \mbox{ for some } x.$$

The reasonable test statistic would be
$$\int \epsilon'(x)\Sigma(x)\epsilon(x) \; dx.$$
It can be shown that the test statistic is approximately a
$c\chi^2$-distribution (Chung, 2001) where the degrees of freedom
and proportionally $d$ is estimated by the method of moment
matching. This is a global inference. For a local inference, we
take $\sup_{x \in \Omega_1} \epsilon'(x)\Sigma(x)\epsilon(x)$ to
be the test statistic.

If there are two groups of geometric objects: $n$ objects
$\Omega_i$ embedded in image $I_i$ and $m$ objects $\Omega_i'$
embedded in image $I_i'$. We can model the displacement of the
first group
$$u_i(x) = \mu_1(x) + \Sigma_1^{1/2}(x)\epsilon_i(x)$$
and those of the second group
$$u_j'(x) = \mu_2(x) + \Sigma_1^{1/2}(x)\epsilon_j(x).$$
where $\epsilon_i$ would be zero mean i.i.d. Gaussian random
vector fields and $u_i$ and $u_j'$ are the displacement to a
template. If we use one subject as a template, effectively we are
losing one degree of freedom. So for small sample, it would be
better to register to the independent template. Then we would be
interested in testing if the group means are identical:
$$H_0: \mu_1(x) = \mu_2(x) \mbox{ for all } x.$$



\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.5]{ima38-3.eps}
\caption{local intensity correlation in $5 \times 5$ window (left)
and $10 \times 10$ window (right). By image registration, it is
needed to increase the correlation uniformly in the image.}
\end{figure}








\end{enumerate}



\end{document}


\subsection*{\bf One Dimensional Version of Multi-Step Intensity Matching}
Let $f(x)$ and $g(x)$ be two functions from $R$ to $R$. Suppose
that the two functions are similar in shape such that
$$f(x) = g(x-u(x)) + \epsilon(x),$$
where $\epsilon(x)$ is a random Gaussian noise. Then, by
linearizing $g(x-u(x))$ with respect to $u(x)$,

$$ f(x) - g(x-u(x))= f(x) - g(x) + g'(x)u(x) +\epsilon$$

$$u(x) = \frac{g(x)-f(x)}{g'(x)}$$

Hence, by solving the recurrence
$$u_{n+1}(x)= \frac{g(x)-g(x-u_n(x))}{g'(x)}$$
with the initial condition $$u_1(x)= \frac{g(x)-f(x)}{g'(x)}$$


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.

Let $f(x)$ and $g(x)$ be two functions of similar shape such that
$$f(x)=g(x-u) + \epsilon(x),$$
where $\epsilon(x)$ is a Gaussian random noise. By linearizing the
function $g(x-u)$,
$$f(x)= g(x) -g'(x)u + \epsilon(x),$$
We can then estimate $u$ as a least square estimator $\hat u_1$.
Since the equation is translation invariant, we have
$$g(x)= f(x+u) - \epsilon(x).$$
Again, by linearizing the function $g(x+u)$,
$$g(x) = f(x) + f'(x)u -\epsilon(x),$$
Then subtracting these two linearized equations, we get
$$f(x)-g(x)= -\frac{f'(x)+g'(x)}{2}u + \epsilon(x).$$
The least square estimator $\hat u$ based on the last equation is
better than $\hat u_1$ (I didn't write down the proof yet but it
looks like easy one.)


Let $I_1$ and $I_2$ be the intensities of brain images $\Omega_1$
and $\Omega_2$ respectively. We want to find displacement vector
field $u$ from $\Omega_1$ to $\Omega_2$. Assuming the deformation
to be small, the intensity and the intensity gradient should match
within a certain error. Let $t=(t_1,t_2,t_3)$.
\\
1. Intensity Matching\\
\bq
I_1(t) = I_2(t+u(t)) + \epsilon(t)\\
       \simeq  I_2(t)+ I_2'(t)u(t) + \epsilon(t)\eq
,where $I_2'(t)$ is the gradient of the intensity


2. Intensity Gradient Matching\\
$
I_1'(t)=I_2'(t + u(t)) + \epsilon(t) \\
\simeq (I_2(t)+ u'(t))^{T}I_2'(t+u(t)) +\epsilon(t)\\
\simeq (I_2(t)+u'(t))^{T}(I_2'(t)+I_2''(t)u(t)) + \epsilon(t)\\
=I_2'(t)+I_2''(t)u(t)+ u'(t)^{T}I_2'(t)+u'(t)^{T}I_2''(t)u(t)
+\epsilon(t)
$\\
Neglecting higer order terms, we have \bq I_1'(t)=I_2'(t) +
I_2''(t)u(t)+ \epsilon(t) \eq Assuming the 3 by 3 matrix
$I_2''(t)$ to be invertible,
the displacement vector field can be expressed as\\
\begin{equation}
u(t)= I_2''(t)^{-1}(I_1'(t)-I_2'(t))
\end{equation}

\subsection*{\bf Estimating a Shift Function}
If $u$ is not a parameter but a function of $x$, then
$$f(x) = g(x-u(x)) + \epsilon(x).$$
By linearizing $g$,
$$ f(x) - g(x-u(x))= f(x) - g(x) + g'(x)u(x) +\epsilon$$
Then the first approximation is given by
$$u_1(x)= \frac{g(x)-f(x)}{g'(x)}$$
The function $g(x)$ will be deformed to $g(x-u_1(x))$. Match the
deformed function $g(x-u_1(x))$ to $f(x)$ using the same
linearization. Iteratively,
$$u_{n+1}(x)= \frac{g_n(x)-f(x))}{g_n'(x)}$$
$$g_n(x) = g(x-u_n(x)-\ldots -u_1(x))$$
The estimated deformation of $u(x)$ after the $n$-th iteration is
$$\hat u(x) = u_1(x) + \ldots u_n(x)$$
The derivative $g(x)$ is not well defined when $g(x)$ is not a
$C^1$ function. If $g(x)$ is a discrete function defined on finite
number of points, then we may use the expected first derivative in
the neighborhood of $x$.

 \item{\em Multiscale data.} 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. One of the advantage of the multiscale data representation
is the differentiability.



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