\documentstyle[12pt,fullpage,psfig]{article}

\newcommand{\lecturetitle}[4]{
  \noindent
  CS 280, Computer Vision (Prof. Jitendra Malik) \hfill #3\\
  Lecture \##1 (#2) \hfill {\bf DRAFT} Notes by #4\\[-.05in]
  \mbox{}\hrulefill \mbox{}\\}

\begin{document}

% MODIFY THE FOLLOWING TO INCLUDE THE LECTURE #, TOPIC, DATE, and YOUR NAME
\lecturetitle{8}{Optical Flow Estimation}{February 17, 1999}{Mark A.
Paskin}

% INCLUDE AN OUTLINE OF THE CLASS TOPICS
\begin{itemize} 
  \item Review: The Optical Flow Problem
  \item Computing the Error of an Estimate
  \item Computing the Estimate with Minimum Error
  \item An Aside: Geometry of Matrices
\end{itemize}

% THE NOTES WILL GO HERE

\section{Review: The Optical Flow Problem}

Last time, the optical flow equation was introduced:
\begin{eqnarray}
{\bf \dot{q}} &=& \lambda {\bf \widehat{q}}^2 {\bf t} + {\bf \widehat{q}} {\bf \omega} \label{optflow}
\end{eqnarray}
where ${\bf q}$ is the direction, $\lambda$ is the nearness ($\frac{1}{||{\bf x}||}$), ${\bf t}$ is the translation, and ${\bf \omega}$ is the rotation.  Given a set of $N$ point-velocity pairs 
$$
\{ ({\bf q}_1, {\bf \dot{q}}_1), ({\bf q}_2, {\bf \dot{q}}_2), \ldots , ({\bf q}_N, {\bf \dot{q}}_N) \}
$$
(which might be computed from two sets of image points and a correspondence mapping), we would like to ask the question, ``What is the most likely translation ${\bf t}$, rotation ${\bf \omega}$, and reduced nearnesses $\mu_i$ for the points ${\bf q}_i$ which would cause the observed velocities ${\bf \dot{q}}_i$?''\footnote{Recall that the reduced nearness $\mu_i$ is the value of the nearness $\lambda_i$ when ${\bf t}$ is normalized.}  In the previous lecture, we saw that with five such pairs, the ${\bf t}$, ${\bf \omega}$, and $\mu_i$ are exactly determined; however, in the presence of noise, an exact solution is likely to be incorrect.  To insure that our estimation of the actual ${\bf t}$, ${\bf \omega}$, and $\mu_i$ is good, we can estimate the ``best'' solution to the overdetermined system arising from more than five pairs.  To do this, we will use least-squares optimization.

\section{Computing the Error of an Estimate}

Define ${\bf p}_i$ as ${\bf \dot{q}}_i$, the measured parallax, and define ${\bf p}_i'$ as the value of ${\bf \dot{q}}$ in equation (\ref{optflow}), the predicted parallax given an estimate of $\left< {\bf t}, {\bf \omega}, \mu_i \right>$.  Then we can define the error of our estimate of $\left< {\bf t}, {\bf \omega}, \{ \mu_i \, | \, 1 \le i \le N \} \right>$ as follows:
\begin{eqnarray}
E &=& \frac{1}{N} \sum_{i = 1}^{N} || {\bf p}_i - {\bf p}_i' ||^2
\end{eqnarray}

\section{Computing the Estimate with Minimum Error}
Since $E$ is a function of our estimate $\left< {\bf t}, {\bf \omega}, \{ \mu_i \, | \, 1 \le i \le N \} \right>$, we can find its global minima by locating the points where its gradient is zero.  To do this, we must find the partial derivatives of $E$ with respect to our parameters:
\begin{eqnarray}
\frac{\partial E}{\partial {\bf t}} &=& \frac{1}{N} \sum_{i = 1}^N 2 \mu_i \widehat{{\bf d}}_i^2 ({\bf p}_i - {\bf p}_i') \label{par-t} \\
\frac{\partial E}{\partial \mu_i} &=& 2 (\widehat{{\bf d}}_i^2 {\bf t})^T ({\bf p}_i - {\bf p}_i') \label{par-m} \\
\frac{\partial E}{\partial {\bf \omega}} &=& \frac{1}{N} \sum_{i = 1}^N 2 \widehat{{\bf d}}_i ({\bf p}_i - {\bf p}_i') \label{par-o}
\end{eqnarray}
We want to find an estimate of $\left< {\bf t}, {\bf \omega}, \{ \mu_i \, | \, 1 \le i \le N \} \right>$ such that
$$
\frac{\partial E}{\partial {\bf t}} = \frac{\partial E}{\partial \mu_i} = \frac{\partial E}{\partial {\bf \omega}} = 0
$$
where ${\bf d}_i = {\bf q}_i$.  We must be careful to add the constraint that $||{\bf t}|| = 1$ (since we are using reduced nearness values).

Now, all the ${\bf d}_i$ and ${\bf p}_i$ are known; they are the $({\bf q}_i, {\bf \dot{q}}_i)$ pairs given to us.  The ${\bf p}_i'$ are unknown, as in them are the hidden unknowns ${\bf t}$, ${\bf \omega}$, and $\mu_i$.  Equations (\ref{par-t}) and (\ref{par-m}) are nonlinear in ${\bf d}_i$, so linear optimization techniques cannot be applied.  Instead we will use a nonlinear iterative technique; we will begin with a guess for $\left< {\bf t}, {\bf \omega}, \{ \mu_i \, | \, 1 \le i \le N \} \right>$ and iterate until the error $E$ stays within some preset tolerance.

\subsection{The Iteration Step}

Substituting the RHS of equation (\ref{optflow}) for ${\bf p}_i'$ in equation (\ref{par-t}) and setting the partial derivative equal to zero, we get
\begin{eqnarray*}
\frac{1}{N} \sum_{i = 1}^N \left[ -\mu_i{\bf p}_i - \mu_i(\mu_i{\bf t} - \mu_i({\bf d}_i^T{\bf t}){\bf d}_i + \widehat{{\bf \omega}}{\bf d}_i \right] = 0
\end{eqnarray*}
Rearranging the terms, we have
\begin{eqnarray*}
\frac{1}{N} \sum_{i = 1}^N  -\mu_i^2{\bf t} &=& -\frac{1}{N} \sum_{i = 1}^N \mu_i{\bf p}_i + \frac{1}{N} \sum_{i = 1}^N \mu_i^2({\bf d}_i^T{\bf t}){\bf d}_i - \frac{1}{N} \sum_{i = 1}^N \widehat{{\bf \omega}}\mu_i{\bf d}_i
\end{eqnarray*}
Recognizing each of the summations as an average over all $N$ pairs, we can rewrite this as
\begin{eqnarray*}
\left< \mu^2 \right> {\bf t} &=& -\left<\mu{\bf p}\right> + \left<\mu^2({\bf d}^T{\bf t}){\bf d} \right> - \widehat{{\bf \omega}} \left<\mu{\bf d}\right>
\end{eqnarray*}
where $\left< \alpha_i \right>$ is shorthand for the average of all $\alpha_i$, $1 \le i \le N$.  If we let $\xi = 1/\left< \mu^2 \right>$, then we have an equation which can be used to ``predict'' a new value of ${\bf t}$, given a current estimate:
\begin{eqnarray}
{\bf t}_\textrm{\tiny new} &=& \xi \left( -\left<\mu{\bf p}\right> + \left<\mu^2({\bf d}^T{\bf t}_\textrm{\tiny old}){\bf d} \right> - \widehat{{\bf \omega}} \left<\mu{\bf d}\right> \right) \label{t-predictor}
\end{eqnarray}
In a similar fashion, we can derive predictors of ${\bf \omega}$ and the $\mu_i$:
\begin{eqnarray}
{\bf \mu_i}_\textrm{\tiny new} &=& \frac{-{\bf t}^T({\bf p}_i + \widehat{{\bf \omega}}{\bf d}_i)}{1 - ({\bf d}_i^T{\bf t})^2} \label{mu-predictor}\\
{\bf \omega}_\textrm{\tiny new} &=& \frac{\left<\widehat{{\bf p}}{\bf d} \right> + \widehat{{\bf t}}\left< \mu \widehat{{\bf d}} \right>}{\left< I - {\bf d}{\bf d}^T \right>} \label{omega-predictor}
\end{eqnarray}
Equations (\ref{t-predictor}), (\ref{mu-predictor}), and (\ref{omega-predictor}) give us the means to obtain better and better approximations to the estimate $\left< {\bf t}, {\bf \omega}, \{ \mu_i \, | \, 1 \le i \le N \} \right>$ which minimizes our error function $E$.  Given some initial estimate, we can evaluate each of equations (\ref{t-predictor}), (\ref{mu-predictor}), and (\ref{omega-predictor}) to obtain a better estimate, repeating the process until the error function decreases by less than some preset tolerance.  The only remaining question is how to obtain the initial estimate to start the procedure.

\subsection{Initialization}

First we observe that once we fix ${\bf t}$ and ${\bf \omega}$, all of the $\mu_i$ are determined by equation (\ref{optflow}).  One approach to estimating ${\bf t}$ and ${\bf \omega}$ is to pick five of the point-velocity pairs and solve explicitly for ${\bf t}$ and ${\bf \omega}$.  However, the error function $E$ may have local minima, and it is possible this method will not find the global minima.  In practice, it is better to choose ${\bf t}$ and ${\bf \omega}$ by assuming one of the two is zero; i.e., estimate the translation component from the data as if there were no rotation, and vice-versa for the rotational component.

To initialize ${\bf \omega}$ in this way, we can simply set ${\bf t} = {\bf 0}$ in equation (\ref{omega-predictor}) and solve for ${\bf \omega}$.  This is possible since the predicted value of ${\bf \omega}$ does not rely on a previous prediction of its value.  However, the predicted value of ${\bf t}$ in equation (\ref{t-predictor}) does depend upon a previous prediction of ${\bf t}$; so we must find another way to initialize the value of ${\bf t}$.

Note that when ${\bf \omega} = {\bf 0}$, the parallax vectors ${\bf p}_i$ arise purely from translation.  Consider a single parallax vector ${\bf p}_i$ on the tangent plane to the viewing sphere, and its corresponding distance vector ${\bf d}_i$.  As we are assuming there is no rotational component to the parallax ${\bf p}_i$, it must be true that ${\bf t}$, ${\bf p}_i$ and ${\bf d}_i$ all lie in a single plane.  Thus, if we contruct ${\bf v}_i$ as a normal to this plane, i.e., ${\bf v}_i = {\bf p}_i \land {\bf d}_i$, then ${\bf t}$ will have no projection on ${\bf v}_i$, i.e., ${\bf v}_i \land {\bf t} = 0$.  Since this will only be true when ${\bf \omega} = {\bf 0}$, we can use
\begin{eqnarray*}
E_{\bf t} &=& \sum_{i = 1}^N \left({\bf t} \cdot {\bf v}_i \right)^2 \label{e_f}
\end{eqnarray*}
as an error function, and minimize it.  Manipulating this equation, we see
\begin{eqnarray*}
E_{\bf t} &=& \sum_{i = 1}^N \left({\bf t}^T{\bf v}_i \right)\left({\bf v}_i^T{\bf t} \right) \\
&=& \sum_{i = 1}^N {\bf t}^T\left({\bf v}_i{\bf v}_i^T \right){\bf t} \\
&=& {\bf t}^T\left(\sum_{i = 1}^N {\bf v}_i{\bf v}_i^T \right){\bf t}
\end{eqnarray*}
The ${\bf t}$ which minimizes $E_{\bf t}$ is the eigenvector of $\sum_{i = 1}^N {\bf v}_i{\bf v}_i^T$ with the smallest eigenvalue!  Thus, to compute our initial estimate of ${\bf t}$, we can form the matrix $\sum_{i = 1}^N {\bf v}_i{\bf v}_i^T$ and find the eigenvector with smallest eigenvalue.  (Note that this matrix is symmetric positive-definite, and so all its eigenvalues will be real.)

\section{An Aside: The Geometry of Matrices}

It is often useful to visualize the geometric transforms represented by two- and three-dimensional matrices.  This can most easily be accomplished by using a {\em gauge figure}, in which a familiar shape, such as the unit circle, is subjected to the transform to illuminate the transform's effects.

For example, the identity matrix $I_2$ maps the unit circle to itself (as one would expect).  Scalings of the identity matrix $\alpha I_2$ are isotropic magnifications, e.g., the matrix
$$
\left[ \begin{array}{cc}
2 & 0 \\
0 & 2
\end{array} \right]
$$
magnifies the unit circle to a circle of radius two.  The matrix
$$
\left[ \begin{array}{cc}
2 & 0 \\
0 & 5
\end{array} \right]
$$
is an anisotropic magnification, since it expands the $y$-axis more than the $x$-axis.  The matrix
$$
\left[ \begin{array}{cc}
0 & -1 \\
1 & 0
\end{array} \right] 
=
\left[ \begin{array}{cc}
\cos(\pi/2) & -\sin(\pi/2) \\
\sin(\pi/2) & \cos(\pi/2)
\end{array} \right] 
$$
rotates the $xy$-plane by 90 degrees, since
$$
\left[ \begin{array}{cc}
\cos(\theta) & -\sin(\theta) \\
\sin(\theta) & \cos(\theta)
\end{array} \right] 
$$
is a rotation by $\theta$.

\end{document}
