%* Last edited: Apr 14 12:39 1997 (gideon)
% this version includes comments from Eric after the paper was shipped to
% CVPR review.
% Based on direct-method-cvpr97.ver5.tex but formated for AI memo
% This is the camera ready copy submitted to cvpr97 proceedings based
% on the AI-MEMO but with many changes and new images.
% New version for journal paper based on CVPR97 final paper
% New version for journal paper with stuff added from the thesis. 

\documentstyle[IEEEtran-gideon,epsf]{article}
%\documentstyle[12pt,epsf]{article}
%\documentstyle[times,art12,twocolumn,epsf]{article}
% \documentstyle[twocolumn,epsf]{article}
% \pagestyle{empty}
%set dimensions of columns, gap between columns, and paragraph indent 
\setlength{\textheight}{9.5in}
\setlength{\columnsep}{2.0pc}
\setlength{\textwidth}{7.2in}
\setlength{\footheight}{0.0in}
\setlength{\topmargin}{-0.25in}
\setlength{\headheight}{0.0in}
\setlength{\headsep}{0.0in}
\setlength{\oddsidemargin}{-.35in}
\setlength{\parindent}{1pc}

\makeatother

\def\single{\def\baselinestretch{1.0}\large\normalsize}
\def\double{\def\baselinestretch{1.5}\large\normalsize}
\def\triple{\def\baselinestretch{2.0}\large\normalsize}

%\newcommand{\gideonHome}[0]{/home/ai/gideon}
\newcommand{\gideonHome}[0]{.}
\newcommand{\bibDir}{\gideonHome /tex/bib/}
\newcommand{\PlotsDir}[0]{plots/}
\newcommand{\TbcImagesDir}[0]{images/}
\newcommand{\TbcHandheldImagesDir}[0]{images/handheld/}
\newcommand{\TbcJamesImagesDir}[0]{images/james1/}
\newcommand{\TbcBneiImagesDir}[0]{images/bnei1/}


\newtheorem{proposition}{Proposition}
\newtheorem{corollary}{Corollary}
\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{result}{Result}
\newtheorem{definition}{Definition}
\newenvironment{proof}{{\em Proof:}}{\hfill\rule{2mm}{2mm}}
\def\eop {{\noindent\framebox[0.5em]{\rule[0.25ex]{0em}{0.75ex}}}}
\def\be {\begin{equation}}
\def\ee {\end{equation}} 
\def\bea {\begin{eqnarray}}
\def\eea {\end{eqnarray}} 


\newcommand{\mbf}[1]{\mbox{\boldmath $#1$}}
\newcommand{\rank}{\mbox{rank}}
\newcommand{\diag}{\mbox{diag}}
\newcommand{\sign}{\mbox{sign}}

\begin{document}
\title{\Large\bf Model-based Brightness Constraints: on Direct Estimation of Structure and 
Motion\thanks{We acknowledge US-IS BSF contract 94-00120
and the European ACTS project AC074. General support for G.S. comes
from DARPA contracts N00014-94-01-0994 and
95009-5381.}
}
%%for single author (just remove % characters) 
%%\author{Author \\
%%  Department \\
%%  Institute \\
%%  City, STATE, zip}
%%for two authors (this is what is printed) 

\author{\begin{tabular}[t]{cc} 
G. P. Stein  & A. Shashua\\
 \\
        Artificial Intelligence Laboratory & Institute of Computer Science\\
        MIT & Hebrew University of Jerusalem\\
        Cambridge, MA~~02139 & Jerusalem 91904, Israel\\
        gideon@ai.mit.edu & http://www.cs.huji.ac.il/$\sim\rm shashua$/
\end{tabular}}

%%don't want date printed but will use it to add keywords
\date{}
%\date{Keywords: Shape Representation and Recovery, 3D Recovery from 2D,
%  Shape from Motion, Image Sequence Analysis, Algebraic and Projective
%  Geometry }


\maketitle
%I don't know why I have to reset thispagestyle, but otherwise get page numbers 
\thispagestyle{empty}

\subsection*{\centering Abstract}
%IEEE allows italicized abstract
%{\abstract
{\em
  We describe a new direct method for estimating structure and motion from
  image intensities of multiple views. We extend the direct methods of
  \cite{HornWeldon:88} to three
  views. Adding the third view enables us to solve for motion, and
  compute a dense depth map of the scene, directly from
  image spatio-temporal derivatives in a linear manner without first
  having to find point 
  correspondences or compute optical flow.

  We describe the advantages and limitations of this method which are then
  verified with experiments using real images. 
}

\single

\section{Introduction}\label{sec:introduction}

The geometry of multiple views of a 3D scene is well understood. There
exist geometric constraints, which relate corresponding points and
lines in multiple views to the camera geometry. These constraints take
the form of the trilinear tensor equations for three views, and the
epipolar constraints for two views. Given a set of feature
correspondences we can recover the camera geometry, and in particular
camera motion (or the relative position of multiple cameras). But
finding correspondences is a hard problem, and the features must be
recovered accurately in order to correctly recover camera motion.

Feature correspondence, whether optical flow or discrete features, is
based in some form or another, on the constant brightness
assumption. That is, we assume the same feature in two views share
some photometric properties.  In it's strictest form, the {\em
constant brightness constraint}, assumes that the brightness of the
corresponding point does not change between views. Alternatively, we
might look at some function of the brightness, possibly a nonlinear
function such as  brightness edges. The constant brightness constraint is a
good approximation for many surfaces in the real world, especially if
the motion is small.  However, the constant brightness
constraint is not strong enough to give us true correspondences. Local
measurements, for example, cannot give 'optical flow' but only {\em
normal flow}, that is, the image flow estimates in the direction of
the image brightness gradient. This is a particular problem in scenes
with long, nearly straight, edges. We will show some examples later in
the paper.

In this paper we present  the {\em model-based brightness
constraints}, where we combine geometric motion models with the
brightness constraint.  This provides us with a stronger constraint
that can be used for direct estimation of structure and
motion. We can threreby bypass the correspondence problem, and recover the
camera motion directly from the image brightness values.  This results
in accurate motion estimates and good 3D reconstruction in very
challenging scenes.

We now provide the key ideas behind the new method. The {\em optical
flow constraint equation} \cite{Horn-Schunk81} provides a matching
constraint between a point in one image, and a line passing through
the corresponding point in the second image. In other words, for every
point in one image it gives us the equation of a line along which the
corresponding point must lie in the second image. By {\em point} we
refer to every pixel in the image which has a none zero brightness
derivative. The equation is given in terms of the spatio-temporal
derivatives of the image brightness, and the line is parallel to the
iso-brightness contour in the first image. The `optical flow
constraint equation' is a first order approximation and assumes that
the image motion (and hence typically, the camera motion) is 
small.

There are no geometric constraints on point-line correspondences
between two views. This can be seen using the following reasoning. A
point in the image defines a line in space. A line in the second image
defines a plane in space. A line and a plane always
intersect. Therefore, given a set of points in one image and a set of
corresponding lines in the second image, for every camera geometry
there exists a valid 3D interpretation and no constraint on the
geometry exists.  We
must therefore use three view geometry, where a point in one view, and
lines through the corresponding points in two other views, provide a
constraint, which can be written in the form of the {\em trilinear
tensor}
\cite{Shashua-tri-pami,Hartley:95,Spetsakis-Aloimonos-iuw90}. The 27
coefficients of the `trilinear tensor' encapsulate the camera motions
and the internal parameters of the camera such as the focal
length. The new method can therefore be viewed as a necessary
extension of the `direct methods' of Horn \& Weldon
\cite{HornWeldon:88} from two views (one motion) to three views (two
motions). These methods are dubbed 'direct
methods' because they do not require prior computation of optical flow.

By combining the optical flow constraint equation with the geometric
model of the `trilinear tensor' we obtain the {\em tensor brightness
constraint} \cite{Stein-Shashua-cvpr97} that describes the
relationship between the spatio-temporal brightness derivatives at
each pixel in the image, with the camera parameters modeled by the 27
coefficients of the trilinear tensor. This `tensor brightness
constraint' provides one linear equation per pixel in the image which
results in a highly over-constrained set of equations. The `tensor
brightness constraint' is valid for the most general case (projective)
where the cameras undergo general motion and we do not know the
internal camera parameters which might vary from frame to frame.

By combining the constant brightness constraint  \cite{Horn-Schunk81}
with the trilinear tensor
\cite{Shashua-tri-pami,Hartley:95,Spetsakis-Aloimonos-iuw90} we obtain
a model-based brightness constraint that describes the relationship
between the spatio-temporal derivatives at each pixel in the image
with the ego-motion parameters modeled by the 27 coefficients of the
trilinear tensor. Hence, the model-based brightness constraint
provides a linear constraint per pixel in the image which results in a
highly over-constrained set of equations.
We then proceed
through a hierarchy of reduced motion models first by assuming
calibrated cameras 
and then by assuming the Longuett-Higgins \& Prazdny small motion model
\cite{Higgins-Prazdny80} resulting in reduced model-based brightness
constraints for those motion models.

%  Many natural scenes such as tree branches or
%man made objects such as window frames, lamp posts and fences often
%give rise to these problems.

%Section (\ref{sec:details}) describes some of the implementation
%details needed to make the method work.  Section
%(\ref{sec:real.images}) shows the results with real images. The method
%is shown to produce useful results in both camera motion and depth
%estimation. Section (\ref{sec:future.work}) discusses some of the open
%questions and suggests possible solutions.

\subsection{Related Work}
\subsubsection{Correspondence Based Methods}

The standard approach to the problem of structure from motion is to
first compute correspondences. These might be either dense
correspondence in the form of optical flow
\cite{Horn-Schunk81,LucasKanade:81,Nagel-opticflow92} or feature correspondence
\cite{Deriche-Giraudon90,Weng-cvpr89,Faugeras93}. Then the
correspondences are used to compute the camera motion and scene
structure.  The advantage of our method over both optical flow methods
\cite{Longuet-Higgins-Prazdny:80,LucasKanade:81,Heeger-Jepson-iccv90,Reiger-Lawton85,Tian-et-al96}
and feature based methods 
\cite{Tomasi-Kanade92,Faugeras93,ShashuaWerman:95,Soatto-Perona96,Azarbayejani95}
is that no prior computation of correspondences 
is needed, a computation which in itself is error prone. Since we obtain a
linear system of equations that combines the information from all the
pixels in the image, we avoid the aperture problem without having to
apply a smoothness assumption. By avoiding the need
to explicitly detect feature points we can use information from areas
where gradient information is weak such as shadow edges. 

These advantages are highlighted in a scene such as in
Figure~\ref{fig:lab213.input.ps}. The scene contains a background of
straight bars, a plaster bust in the foreground, and an oblique line
cutting the bars.  The background bars provide an unreliable source of
information for point-to-point correspondence because of the aperture
problem. Even large regions of the image have no information which
constrains the flow in the $Y$ direction. Many feature points in the
image do not correspond to real feature points in the 3D world. For
example, the intersection of the bars and the line in the image do not
correspond to physical features in 3D since the line does not lie on
the background plane. Figure~\ref{fig:TBC.LK.compare.flow.only} shows
a detail of the optical flow field in the upper left of the scene
where these two problems occur. The flow was computed using a widely
available, `industrial strength', optical flow program
\cite{Bergen-Hingorani90,Anandan-eccv92}. (For more examples see
Section~(\ref{sec:TBC.optical.flow.compare}).) 

Finally, the texture on the plaster bust varies smoothly and is low
contrast thereby giving rise only to a relatively small number of
reliable features to track. Yet, our method produces a full 3D model
of the scene as demonstrated in Figure~\ref{fig:lab213.bump.ps}{\it
a,b}.  Many natural scenes such as tree branches or man made objects
such as window frames, lamp posts and fences, often give rise to these
problems.

\begin{figure}
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/lab212.s.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/lab213.s.ps}
              \hfill }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/lab210.s.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/lab213.depth.s.ps}
              \hfill }
                    \centerline{ 
              \hfill (c) \hfill \hfill (d) \hfill 
              } 
  \caption[The three input images and the estimated depth map]{\em 
    The three input images (a,b,c) and the estimated depth map
  (d). The motion between images (a) and (b) is horizontal. The motion
  is small but can be seen in the width of the rightmost vertical stripe
  and in the parallax between the cylinder on the right and the vertical
  stripes behind it. The motion of between images (a) and (c) was
  vertical. }\label{fig:lab213.input.ps}
\end{figure}


\begin{figure}
          \centerline{ 
               %\hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/flowDetail.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/lkFlowDetail.ps}
              \hfill }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
\protect\caption[Detail of optical flow computed for
Figure~\protect\ref{fig:lab213.input.ps}]{\em Detail of optical flow
computed for Figure~\protect\ref{fig:lab213.input.ps}. (a) Correct
optical flow computed using motion and depth map recovered by the
direct method. The camera motion was horizontal parallel to the image
plane. The flow vectors have a zero $Y$ component. (b) Optical flow
computed using Bergen \& Hingorani optical flow program. Note in the
upper left corner the flow vectors have a strong $Y$ component. This
is due to the aperture problem where even a large window (aperture)
will see image gradients in only one direction giving no constraint on
the $Y$ component of the flow.  Near the intersections (in the image)
of the diagonal line and the vertical bars the flow has a small $Y$
component. This error occurs because the program 'tracks' the
intersection point as if it were a real 'feature' point, but the lines
do not intersect in space. Although smaller in magnitude this is the
more significant error because the program will also give a high
confidence to this value. }\label{fig:TBC.LK.compare.flow.only}
\end{figure}

\begin{figure}
  \centerline{ 
        \epsfxsize=1.6in \epsffile{lab4-images/lab213.smooth.ps}
               \hfill 
        \epsfxsize=1.6in \epsffile{lab4-images/lab213.no-tex.ps}
        }
  \centerline{ 
              \hfill (a) \hfill \hfill (b)\hfill 
              } 
\protect\caption[3D rendering of the estimated depth map from images in
Figure~\protect\ref{fig:lab213.input.ps}]{\em 
    3D rendering of the estimated depth map from images in
Figure~\protect\ref{fig:lab213.input.ps}. These images show the
inverse depth, $k=\frac{1}{z}$, which is the natural value to
compute. In (b) the texture map was removed to show the detail and the
flaws. }\label{fig:lab213.bump.ps}
\end{figure}

\subsubsection{Direct Methods}

The `direct methods' were pioneered by Horn \& Weldon in
\cite{HornWeldon:88}. Using only a single image pair they ended up
with $N$ equations in $N+5$ unknowns, where $N$ is the number of
points in the image. The unknowns are the $N$ unknown depths and the
$3$ translation and $3$ rotation parameters with one unknown dropping
out because translation and depth can only be found up to a common
scale factor. The problem is therefore ill posed and additional
constraints are needed. Negahdaripour \& Horn \cite{Negahdaripour:87}
present a closed form solution assuming a planar or quadratic surface.
Szeliski \& Kang \cite{Szeliski-Kang-95} describe an iterative
solution using splines to enforce a smoothness constraint on the
depth.  McQuirk \cite{Ignacio:96} shows that in a pure translation
model the subset of the image points with a nonzero spatial derivative
but a zero time derivative gives the direction of motion, thus the
{\em focus of expansion} (FOE) is on a line perpendicular to the
gradient at these points. Since this uses only a subset of the image
points, one ends up discarding most of the image data.

Heel \cite{Heel:90} constructs a Kalman filter to build up a structure
model from more than one image pair but the core computation is
fundamentally the same single image pair computation. The basic idea
is that the Horn \& Weldon equation is linear in depth if motion is
known and linear in motion if depth is known. Heel initially assumes a
depth (constant depth surface) and estimates the motion and then uses
the motion to estimate the depth. The depth map is then warped using
the computed motion and used as the initial depth estimate for the
next image pair in the sequence. The question of convergence is not
answered and failures are not reported. Results are only shown for
pure translation. Heel limits the motion to image motions of less than
3 pixels which is one of the limiting factors in the accuracy of the
surface reconstructions.  Within a coarse-to-fine implementation
(Section~\ref{sec:TBC.coarse.to.fine}) our method can handle much
larger image motions averaging up to 50 pixels for $640\times 480$
resolution images, thereby increasing the dynamic range by an order of
magnitude..



Michaels \cite{Michaels-thesis92} uses three frames (two motions). Each
motion gives $N$ equations totaling $2N$ equations. The unknowns are
the $N$ unknown depths and the $2 \times 6$ translation and rotation
parameters. Again, one parameter drops out because of the translation
and depth scale ambiguity.  He solves the $2N$ equations with $N+11$
unknowns as a large nonlinear optimization problem using the
Levenberg-Marquart algorithm \cite{Minpack80}. Since it took a very
long time to converge results are shown for very low resolution images
only. He also shows theoretically that a large field of view ($\approx
120^o$) is required for accurate direct estimation of motion given
structure. Since estimation of motion given structure is a key stage
in the process he suggests that a wide field of view is required for
the whole process of estimation of structure and motion.  Kumar et al
\cite{Kumar-Anandan94-icpr} first align a dominant plane in the
images. The residual image motion is epipolar motion (i.e. pure
translation). The epipole of the residual motion is found using
iterative techniques starting with an initial guess. Others who have
explored planar alignment include
\cite{Irani-et-al94,Sawhney94-cvpr}. The theoretical framework for
planar alignment can also be found in \cite{Shashua-Navab94-cvpr}. 

More recently, Fermuller \& Aloimonos
\cite{Fermuller-Aloimonos-iccv95,Fermuller-Aloimonos-cvpr97} 
describe global geometric properties of the flow-field that give rise
to direct relationships between the measurement of normal flow and the
ego-motion parameters solved by means of search techniques. They apply
the {\em depth is positive} constraint and try to find an optimally
smooth surface.

The tensor brightness constraint was first presented in
\cite{ShashuaHanna:95} and a practical implementation with results was
first described in \cite{Stein-Shashua-cvpr97}. In this paper we
present new theoretical results which lead to a modification of the
algorithm. The modified algorithm gives accurate motion estimates even
in the presence of considerable camera rotation. We present new
quantitative and qualitative results including some with handhled
cameras and outdoor scenes. 

\section{Mathematical Background}
\label{sec:math}

\subsection{Notation}

We will use uppercase bold to denote matrices (e.g. ${\bf
A}$). Vectors representing 3D points which will be upper case and not
bold. Other vectors and scalars will be in lower case.

We will occasionally use tensorial notations.  We
use the covariant-contravariant summation convention: a point is an
object whose $n$ coordinates are specified with superscripts, i.e.,
$p^i=(p^1,p^2,...,p^n)$. These are called contravariant vectors. An
element in the dual space (representing hyperplanes --- lines in
${\cal P}^2$), is called a covariant vector and is represented by
subscripts, i.e., $s_j=(s_1,s_2,....,s_n)$. Indices repeated in covariant
and contravariant forms are summed over, i.e., $p^is_i = p^1s_1 +
p^2s_2+...+p^ns_n$. This is known as a contraction.
An outer-product of two 1-valence tensors (vectors), $a_ib^j$, is a
2-valence tensor (matrix) $c_i^j$ whose $i,j$ entries 
are $a_ib^j$ --- note that in matrix form $C=ba^\top$.

Matching image points across three views will be denoted by
$p,p',p''$; the homogeneous coordinates will be referred to as
$p^i,p'^j,p''^k$, 
or alternatively as non-homogeneous image coordinates
$(x,y),(x',y'),(x'',y'')$ --- hence, $p^i=(x,y,1)$, etc. 



\begin{figure}
  \centerline{ \epsfysize=2.5in \epsffile{3view.geom.ps}} \caption[The geometry of three views]{\em
The Geometry of Three Views: Image~1, Image~2 and Image~3 are three
views obtained from cameras centered at $o$, $o'$ and $o''$
respectively. Point $P$ in the scene projects to image points $p$,
$p'$ and $p''$. $s'$ and $s''$ are any image lines passing through the
points $p'$ and $p''$ respectively. }\label{fig:MATH.3view.geom}
\end{figure}

We will now consider three perspective views $\psi$, $\psi'$ and
$\psi''$ of a 3D scene.  Figure~\ref{fig:MATH.3view.geom} shows a 3D
point $P \in {\cal P}^3 $ and its image in the three views $p \in
\psi$, $p' \in \psi'$ and $p'' \in \psi''$.
Without
loss of generality we can align the 3D world coordinate system with
the coordinate system of the first camera:
\be
\begin{array}{lcccccr}
p & \cong &\left(\begin{array}{c}
  x\\
  y\\
  1
\end{array}
\right)&
\cong& \left[ {\bf I} ; 0 \right] P &&
\end{array}
\ee
Thus  $P \cong (x,y,1,\rho)^\top$. 
After we have set the first camera coordinate system, the other two
camera coordinate systems are in general given by:
\be\label{eq:MATH.A}
\begin{array}{lcccccccr}
p' & \cong &\left(\begin{array}{c}
  x'\\
  y'\\
  1
\end{array}
\right)&\cong& {\tilde {\bf A}} P&\cong& \left[ {\bf A} ; t' \right] P
&\cong &{\bf A} p + \rho t' 
\end{array}
\ee
\be\label{eq:MATH.B}
\begin{array}{lcccccccr}
p'' & \cong &\left(\begin{array}{c}
  x''\\
  y''\\
  1
\end{array}
\right)&\cong& {\tilde {\bf B}} P&\cong &\left[ {\bf B} ; t'' \right] P
&\cong &{\bf B} p + \rho t'' 
\end{array}
\ee
The matrices ${\bf A}$ and ${\bf B}$ are homography matrices from
Image~1 to Image~2 and to Image~3 respectively {\em due to the same
plane $\pi$}. The vectors $t'$ and $t''$ are the epipoles, the
projection center of camera 1 projected onto the image planes of the
2nd and 3rd cameras respectively. $\rho$ is the relative affine
depth\cite{Shashua-Navab-pami}. It is important to note that $\rho$ in
independent the second view point and thsu has he same value in
equations~(\ref{eq:MATH.A}) and~(\ref{eq:MATH.B}). In a calibrated
setting, (i.e. if the intrinsic parameters are known), the matrices
${\bf A}$ and ${\bf B}$ are rotations and the vectors $t'$ and $t''$
are the translations. Finally, in the calibrated setting $\rho$ is
replaced by $k = \frac{1}{z}$.

We will now derive the triliner tensor of \cite{Shashua-tri-pami}. 
Let $s'$ and $s''$ be lines through points $p'$ and $p''$
respectively:
\bea
\begin{array}{lccccccr}
s'^\top p'=0 &&&&&&&
s''^\top p''=0 
\end{array}
\eea
Premultiplying the left and right hand sides of~(\ref{eq:MATH.A}) by
$s'^\top$ and of~(\ref{eq:MATH.B}) by $s''^\top$ we get:
\be
\begin{array}{lcr}
s'^\top {\bf A} p + \rho s'^\top t' &=&0\\
s''^\top {\bf B} p + \rho s''^\top t''&=&0
\end{array}
\ee
%\noindent
Eliminating $\rho$ from the above equations results in the equation:
\be
s'^\top t's''^\top {\bf B} p - s''^\top t''s'^\top {\bf A} p=0
\ee
%\noindent
This can be written compactly using tensor notation:
\be
p^i s''_k s'_j  {\cal T}_i^{jk} = 0 \label{eq:MATH.tensor.equation}
\ee
where ${\cal T}$ is the
tensor representing a bilinear
function of the camera matrices:
\be
{\cal T}_i^{jk}=t'^jb_i^k - t''^ka_i^j\label{eq:MATH.tensor}
\ee

Equation~(\ref{eq:MATH.tensor.equation}) relates a point $p$ in Image 1
 and lines  $s'$ and $s''$ passing through the corresponding
points $p'$ and $p''$ in Image 2 and Image 3 respectively. 
It is important to note that the lines $s'$ and $s''$ do not have to
correspond to any physical line in space or  the image.  
These constraints first became prominent in
\cite{Shashua-tri-pami} and the underlying theory has
been studied intensively in
\cite{ShashuaWerman:95,Hartley:95,Shashua-Anandan-iuw96,Faugeras-Mourrain-iccv95,Triggs-iccv95,Heyden-iccv95,Shashua-Avidan-eccv96}.

%Let $\mbf x$ be a point in 3D space and its projection in
%a pair of images be $p$ and $p'$.  Then $p=[I;0]{\mbf x}$ and
%$p'\cong A{\mbf x}$. The left $3\times 3$ minor of $A$ stands for a 2D
%projective transformation of the chosen plane at infinity and the
%fourth column of $A$ stands for the epipole (the projection of the
%center of camera 0 on the image plane of camera 1). In particular, in
%a calibrated setting the 2D projective transformation is the
%rotational component of camera motion and the epipole is the
%translational component of camera motion.



%Three views, $p=[I;0]{\mbf x}, p'\cong A{\mbf x}$ and $p''\cong
%B{\mbf x}$, are known to produce four trilinear forms
%whose coefficients are arranged in a tensor representing a bilinear
%function of the camera matrices $A,B$:
%\be
%\alpha_i^{jk}=t'^jb_i^k - t''^ka_i^j
%\ee
%where $A=[a_i^j, t'^j]$ ($a_i^j$ is the $3\times 3$ left minor and $t'$
%is the fourth column of $A$) and $B=[b_i^k, t''^k]$. The tensor acts on
%a triplet of matching points in the following way:
%\begin{equation}
%p^is_j^\mu r_k^\rho \alpha_i^{jk} = 0 \nonumber
%\end{equation}
%where $s_j^\mu$ are any two lines ($s_j^1$ and $s_j^2$) intersecting at
%$p'$, and $r_k^\rho$ are any
%two lines intersecting $p''$. Since the free indices are $\mu,\rho$
%each in the range 1,2, we have 4 trilinear equations (unique
%up to linear combinations). 

\section{Model-based Brightness Constraints}

\subsection{Photometric Constraints}
Geometrically, a trilinear matching constraint is produced by
contracting the tensor with the point $p$ in Image 1, {\em any} line
coincident with $p'$ in Image 2, and {\em any} line coincident with
$p''$ in Image 3. In particular, we may use the tangent to the
iso-brightness contour at $p'$ and $p''$, respectively, and thus one
can recover in principle the camera matrices across three views in the
context of the ``aperture'' problem, as suggested by
\cite{Spetsakis-Aloimonos-iuw90}. However, there still remains the
problem of finding those matching tangents in the first place. This we
now solve. 

A first order approximation of the constant brightness constraint
leads to the optical flow constraint equation \cite{Horn-Schunk81}:
\begin{equation}\label{eq:TBC.constant.brightness}
u' I_x + v' I_y + I'_t = 0
\end{equation}
where $(u',v')$  
are the {\em optical flow} values at $(x,y)$ between Image 1 and Image
2 (i.e.  $u'=x' - x$
and $v' = y' - y$).  $(I_x, I_y, I'_t)$ are the spatial and temporal
derivatives at the coordinates $(x,y)$. In
practice $I'_t = I_2(x,y)-I_1(x,y)$.

The optical flow constraint
equation~(\ref{eq:TBC.constant.brightness}) can be rewritten in the
form:
$$\left(
I_x,
I_y,
I'_t
\right)^\top
\left(
u,v,1\right) =0.$$
The line $$s=\left(I_x,I_y,-x I_x - y I_y \right)^\top$$ in the projective
plane passes through the point $p=(x,y,1)^\top$ since:
$$
\left(
I_x,I_y,-x I_x - y I_y\right)
^\top
\left(
x,y,1\right) =0.$$ Combining those two equations together:
{\small
\begin{eqnarray}\label{eq:TBC.optical.projective.3}
\begin{array}{rcl}
\left(
\begin{array}{c}
I_x\\
I_y\\
-x I_x - y I_y
\end{array}\right)
^\top
\left(
\begin{array}{c}
x\\
y\\
1
\end{array}\right) + 
\left(
\begin{array}{c}
I_x\\
I_y\\
I'_t
\end{array}\right)
^\top
\left(
\begin{array}{c}
u\\
v\\
1
\end{array}\right) &=&\\
\left(
\begin{array}{c}
I_x\\
I_y\\
I_t'-x I_x - y I_y
\end{array}\right)
^\top
\left(
\begin{array}{c}
x'\\
y'\\
1
\end{array}\right) &=&0
\end{array}
\end{eqnarray}}
Where we have used $x'=x+u$ and $y'=y+v$.
Therefore, the line:
\begin{equation}\label{eq:TBC.optical.projective.s'}
s'=\left(\begin{array}{c}
I_x\\
I_y\\
I'_t -x I_x - y I_y \end{array}\right)^\top
\end{equation}
passes through the
point $p'=(x',y',1)^\top$. 

Thus the photometric constraints provide a matching constraint between a
point $p$ and a line $s'$ passing through the corresponding point
$p'$ in Image~2, and between a
point $p$ and a line $s''$ passing through the corresponding point
$p''$ in Image~3. 
The lines:
\be\label{eq:TBC.s1}
%\begin{array}{ccc}
 s'=\left(\begin{array}{c}
I_x\\
I_y\\
I'_t - xI_x - yI_y
\end{array}\right) 
\ee
and
\be\label{eq:TBC.s2}
 s''=\left(\begin{array}{c}
I_x\\
I_y\\
I''_t - xI_x - yI_y
\end{array}\right)
%\end{array}
\ee 
are lines coincident with $p'$ and $p''$ respectively and parallel
to the iso brightness contour at $(x,y)$. $I''_t$ is the temporal derivative between the Image 3 and Image 1. (i.e. $I''_t =
I_3(x,y)-I_1(x,y)$.) 

\subsection{The Projective Model: The Tensor Brightness Constraint}

Substituting equations~(\ref{eq:TBC.s1}) and~(\ref{eq:TBC.s2}) into the tensor
equation~(\ref{eq:MATH.tensor}) results in the {\em tensor brightness
  constraint}:
\be
\fbox{$s''_ks'_jp^i{\cal T}_i^{jk} = 0.$}\label{eq:TBC.tbc.proj} 
\ee 
We have one such equation for each point on the image where $s''_{k}$ and
$s'_{j}$ can be computed from the image gradients and $p^i=(x,y,1)$ are
the (projective) image coordinates of the point in Image 0. We 
solve for ${\cal T}_i^{jk}$ which combines the motion and camera parameters.
{\em The coordinates of the corresponding points $(x',y')$ and
  $(x'',y'')$ are not required}. 

Every pixel with a non-vanishing gradient contributes one linear
equation to the 27 unknown parameters comprising ${\cal T}_i^{jk}$.
However, a configuration of a point in the first image and two {\it
  parallel\/} lines in Images~2 and~3 is a particular instance of a
degenerate line configuration called a Linear Line Complex. The
general solution for the LLC case is explored
in~\cite{Stein-Shashua-eccv98}. In this particular case the system of
equations can provide a linear solution to 21 of the parameters and
the remaining 6 parameters of ${\cal T}_i^{jk}$ can be determined
using quadratic admissibility constraints. We provide here a proof for
this special case:

From equation~(\ref{eq:TBC.tbc.proj}) the coefficients of the terms
${\cal T}_i^{12}$ and the terms ${\cal T}_i^{21}$ are both $p^i I_x
I_y$. Therefore the linear equations cannot be solved for the six terms (${\cal
T}_i^{12}$ and ${\cal T}_i^{21}$) individually but only for the three sums
$${\cal T}_i^{12} + {\cal T}_i^{21}, i=1...3$$.

The sum ${\cal T}_i^{12} + {\cal T}_i^{21}$ provides one linear
equation in ${\cal T}_i^{12}$ and ${\cal T}_i^{21}$. We will use the
admissibility constraint on the Standard Correlation Slices, $T_i$.
The constraint states that the matrix $T_i$ is of $rank=2$. This leads
to a quadratic equation in ${\cal T}_i^{12}$ and ${\cal T}_i^{21}$.
The quadratic and linear equations together result in two solutions
for each pair ${\cal T}_i^{12}, {\cal T}_i^{21}, i=1...3$ for a total
of 8 discrete solutions.  A unique solution is obtainable by applying
further admissibility constraints as shown
in~\cite{Stein-Shashua-eccv98}.

This added complexity in finding the unique solution leads us to
search for simpler models. We next consider the case of a calibrated
camera with the rotation limited to small rotation angles. 

%\begin{theorem}
%Every point-line-line matching triplet across three images provides one
%linear constraint to $\alpha_i^{jk}$. In case all the matching triplets
%have the two lines parallel to each other, then $\alpha_i^{jk}$ can be
%determined up to an 8-fold ambiguity.
%\end{theorem}

%We will not prove this result here (see \cite{Stein-Shashua-eccv98}
%for proof).  Therefore, employing the constant brightness equation
%within a general three-frame sequence of images provides a linear
%system of equations but which is not fully constrained to solve for
%the tensor uniquely.

\subsection{The Small Rotation Model with Calibrated Cameras}

The next model is defined for small-angle rotations with calibrated
cameras. Assuming calibrated cameras equations~(\ref{eq:MATH.A})
and~(\ref{eq:MATH.B})  
become:
\bea\label{eq:calib.mot.model}
\begin{array}{lcccr}
p' & \cong& \left[ {\bf R}' ; t' \right] P
&\cong &{\bf R}' p + k t' \\
p'' & \cong& \left[ {\bf R}'' ; t'' \right] P
&\cong &{\bf R}'' p + k t'' 
\end{array}
\eea 
i.e., ${\tilde A}=[{\bf R'}; t']$ and ${\tilde B}=[{\bf R''};
t'']$ and $k = \frac{1}{z}$ replaces $\rho$ in
equation~(\ref{eq:MATH.A}) and~(\ref{eq:MATH.B}).

If we also assume small angle rotations and can therefore make the
approximations: 
$$
\begin{array}{lr}
\cos(\theta) \approx 1, & \sin(\theta) \approx \theta
\end{array}
$$
then the rotation matrices can be approximated as:
\bea
\begin{array}{lcccr}
R' &\approx & \left[ I + [w']_x \right] \\
R'' &\approx  & \left[  I + [w'']_x \right]
\end{array}
\eea
where $w',w''$ are the
angular velocity vectors and $[\cdot]_x$ is the skew-symmetric matrix
of vector products. 
Equations~(\ref{eq:MATH.A})
and~(\ref{eq:MATH.B}) now become:
\bea\label{eq:TBC.AB.calib}
\begin{array}{lcccr}
p' & \cong&  {\tilde {\bf A}} P&\cong& \left[ I + [w']_x ; t' \right] P\\
p'' & \cong&  {\tilde {\bf B}} P&\cong& \left[  I + [w'']_x  ; t'' \right] P.
\end{array}
\eea

We can now proceed along two paths. The different results will give us
some further insight into the problem. First we will substitute
${\tilde A}=[ I + [w']_x; t']$ and ${\tilde B}=[ I + [w'']_x; t'']$
directly into equation~(\ref{eq:MATH.tensor}) and then write out the
tensor explicitly:
\bea
T_1 = 
\left[
\begin{array}{ccc}
       t'_1-t''_1 &       t'_1 w''_3-t''_2 &      -t'_1 w''_2-t''_3\\
       t'_2-t''_1 w'_3 &  t'_2 w''_3-t''_2 w'_3& -t'_2 w''_2-t''_3 w'_3\\
       t'_3+t''_1 w'_2 & t'_3 w''_3+t''_2 w'_2& -t'_3 w''_2+t''_3 w'_2
\end{array}
\right]\\
T_2 = 
\left[
\begin{array}{ccc}
 -t'_1 w''_3+t''_1 w'_3 &      t'_1+t''_2 w'_3  &  t'_1 w''_1+t''_3 w'_3\\
      -t'_2 w''_3-t''_1 &            t'_2-t''_2 &       t'_2 w''_1-t''_3\\
 -t'_3 w''_3-t''_1 w'_1 &       t'_3-t''_2 w'_1 &  t'_3 w''_1-t''_3 w'_1
\end{array}
\right]\\
T_3 = 
\left[
\begin{array}{ccc}
  t'_1 w''_2-t''_1 w'_2 &  -t'_1 w''_1-t''_2 w'_2   &      t'_1-t''_3 w'_2\\
  t'_2 w''_2+t''_1 w'_1 &  -t'_2 w''_1+t''_2 w'_1    &     t'_2+t''_3 w'_1\\
       t'_3 w''_2-t''_1  &      -t'_3 w''_1-t''_2     &         t'_3-t''_3
\end{array}
\right]
\eea
where $T_1, T_2$ and $T_3$ are the 
Standard
Correlation Slices \cite{Stein-Shashua-eccv98}. Studying
the tensor terms we notice that the terms are not linearly independent
and the following seven linear equations hold:
\be
\begin{array}{l}
{\cal T}_1^{1,2}+{\cal T}_2^{1,1}+{\cal T}_1^{2,1}-{\cal T}_2^{2,2}=0\\
{\cal T}_2^{2,1}+{\cal T}_1^{2,2}+{\cal T}_2^{1,2}-{\cal T}_1^{1,1}=0\\
{\cal T}_3^{1,1}+{\cal T}_1^{1,3}+{\cal T}_1^{3,1}-{\cal T}_3^{3,3}=0\\
{\cal T}_1^{3,3}+{\cal T}_3^{3,1}+{\cal T}_3^{1,3}-{\cal T}_1^{1,1}=0\\
{\cal T}_2^{3,3}+{\cal T}_3^{3,2}+{\cal T}_3^{2,3}-{\cal T}_2^{2,2}=0\\
{\cal T}_3^{2,2}+{\cal T}_2^{2,3}+{\cal T}_2^{3,2}-{\cal T}_3^{3,3}=0\\
{\cal T}_2^{1,3}+{\cal T}_3^{2,1}+{\cal T}_1^{2,3}+{\cal T}_3^{1,2}+{\cal
T}_1^{3,2}+{\cal T}_2^{3,1}=0
\end{array}
\ee
Therefore, for the calibrated small rotation model, the 27-parameter
Tensor Brightness Constraint reduces to a {\em 20-parameter}
constraint equation. Each point in the image gives one homogeneous
equation resulting in a set of $N$ homogeneous equations, where $N$ is
the number of image points. This can be written in matrix form as: 
$$
{\bf A} x=0 
$$ where ${\bf A}$ is the $N \times 20$ estimation matrix
and $x$ is the 20 parameter vector of the unknown intermediate parameters. There is
a unique non-trivial solution if
the estimation matrix is of $rank=19$. But as in the projective case,
the system of homogeneous equations is degenerate and in this case the
rank of the estimation matrix is 16 not 19.

In \cite{ShashuaHanna:95,Stein-Shashua-cvpr97} a model-based brightness
constraint for the calibrated, small rotation model was derived in a
different way. Expanding equation~(\ref{eq:TBC.AB.calib}):
\bea
\begin{array}{lcr}
p' & \cong& \left[ I + [w']_x \right] p + k t'
\end{array}
\eea
and taking the dot product of the left-hand-side and right-hand-side terms with $s'$
(equation~\ref{eq:TBC.s1}) results in the equation:
\bea
\begin{array}{lcr}
0 & = & s'^\top p + s'^\top [w']_x p + k s'^\top t'
\end{array}
\eea
 This can be simplified to the form: 
\bea\label{eq:TBC.mot1.calib}
\begin{array}{lcr}
 k s'^\top t' + v'^\top w' + I'_t  & = & 0
\end{array}
\eea
by noting that $s'^\top p=I'_t$ and $s'^\top [w']_x p = v'^\top w'$ if
we define: 
$$ v'= p \times s' = \left(\begin{array}{c} -I_y + y(I'_t -
xI_x - yI_y)\\ I_x -x(I'_t - xI_x - yI_y)\\ xI_y - yI_x
\end{array}\right).
$$

In a similar manner, for the second motion:
\bea\label{eq:TBC.mot2.calib}
\begin{array}{lcr}
 k s''^\top t'' + v''^\top w'' + I''_t  & = & 0
\end{array}
\eea
where
$$ 
v''=p\times s'' = \left(\begin{array}{c}
-I_y + y(I''_t - xI_x - yI_y)\\
I_x  -x(I''_t - xI_x - yI_y)\\
 xI_y - yI_x
\end{array}\right). 
$$ 
Multiplying equations~(\ref{eq:TBC.mot1.calib}) by $ s''^\top t''$ and
equations~(\ref{eq:TBC.mot2.calib}) by $ s'^\top t'$ and subtracting we
obtain the following {\em 24-parameter} model-based brightness
constraint as a reduction of the tensor brightness constraint:
\be
I''_t{s'}^\top t' - I'_t{s''}^\top t'' + {s'}^\top[t'{w''}^\top]v'' -
{s''}^\top[t''{w'}^\top]v' = 0\label{eq:TBC.calib.t2}.
\ee
The resulting 24-parameter model (equation~\ref{eq:TBC.calib.t2}) is not as
compact as the 20-parameter model  we obtained
by direct substitution of the small rotation approximation into the
tensor brightness constraint equation. We find that the $24 \times N$
estimation matrix obtained using the 24-parameter model is again of rank=16.
Due to the high degeneracy of the linear equations neither the
24-parameter model nor the 20-parameter model are convenient to work
with. 

%Starting from the general projective motion model, we investigate two
%further reduced models. The next model is defined for small-angle rotations
%i.e., $A=[I+ [w']_x; t']$ and $B=[I+ [w'']_x; t'']$ where $w',w''$ are the
%angular velocity vectors and $[\cdot]_x$ is the skew-symmetric matrix
%of vector products. We obtain the following model-based brightness
%constraint as a reduction of the tensor brightness constraint:

%\be
%I''_t{s'}^\top t' - I'_t{s''}^\top t'' + {s'}^\top[t'{w''}^\top]v'' -
%{s''}^\top[t''{w'}^\top]v' = 0, \label{eq:t2}
%\ee
%where
%$$ v'=p\times s' = \left(\begin{array}{c}
%-I_y + y(I'_t - xI_x - yI_y)\\
%I_x  -x(I'_t - xI_x - yI_y)\\
%xI_y - yI_x
%\end{array}\right)
%$$
%and 
%$$ v''=p\times s'' = \left(\begin{array}{c}
%-I_y + y(I''_t - xI_x - yI_y)\\
%I_x  -x(I''_t - xI_x - yI_y)\\
% xI_y - yI_x
%\end{array}\right)
%$$ 
%The 24-parameter model contains 20 linearly independent parameters
%(instead of 27), yet still is not sufficiently constrained to produce a
%unique solution for the reduced tensor parameters (not proven here).

\subsection{The 15-Parameter, Small Motion Model}

A unique solution is obtained when we reduce the motion model further
to include infinitesimal motion using the model introduced by
Longuet-Higgins \& Prazdny  \cite{Longuet-Higgins-Prazdny:80}. The LH\&P
model assumes in addition to small rotation that $\frac{t'_z}{Z} \ll
1$.  The motion field equations for the first camera motion are then:
\begin{eqnarray}\label{eq:TBC.u.v.lh}
u'=\frac{1}{z}(t'_{1}-x t'_{3}) - w'_{3}y + w'_{2}(1+x^{2}) 
    - w'_{1} x y\\ 
v'=\frac{1}{z}(t'_{2}-y t'_{3}) + w'_{3}x - w'_{1}(1+y^{2}) 
    + w'_{2} x y \nonumber
\end{eqnarray}
By substituting equation~(\ref{eq:TBC.u.v.lh}) into the optical flow
constraint equation:
\begin{equation}
u' I_x + v' I_y + I'_t = 0
\end{equation}
and
rearranging the  terms we obtain:
\begin{equation}\label{eq:horn.weldon}
k s^\top t' + v^\top w' + I'_t = 0
\end{equation}
where $k=\frac{1}{z}$ denotes the inverse depth at each pixel location
and where $s,v$ are defined below:
\be
s=\left(\begin{array}{c}
I_x\\
I_y\\
-xI_x - yI_y
\end{array}\right)
\label{eq:s}
\ee
and
\be 
v = p \times s = \left(\begin{array}{c}
-I_y - y(xI_x + yI_y)\\
I_x  +x(xI_x + yI_y)\\
xI_y - yI_x
\end{array}\right)
\ee
Equation~(\ref{eq:horn.weldon}) which was first derived in
\cite{Negahdaripour:87,HornWeldon:88} can also be
derived directly from equation~(\ref{eq:TBC.mot1.calib}) by making the LH\&P assumptions:  $\frac{t'_z}{Z} \ll 1$, $\frac{w'_y x}{f} \ll 1$
and $\frac{w'_x y}{f} \ll 1$.

{\bf Note} that $v$ is a vector that can be computed from image
brightness, and should not be confused with $v'$ which is the $y$
component of the image flow between Image~1 and Image~2.
 

Similarly for the second motion:
\begin{equation}\label{eq:horn.weldon2}
k s^\top t'' + v^\top w'' + I''_t=0
\end{equation}
By eliminating $k$ from equations~(\ref{eq:horn.weldon}) and~(\ref{eq:horn.weldon2})
we obtain the {\em 15-parameter model-based brightness constraint}:
\begin{equation}\label{eq:the.key.eq}
\fbox{$I''_{t}s^{T}t' - I'_{t}s^{T}t'' + s^{T}[t'w''^{T} - t''w'^{T}]v = 0$}
\end{equation}
Equation~(\ref{eq:the.key.eq}) is a key equation. There is one such
equation for every point in the image. The unknowns are the motion
parameters $t'$, $t''$, $w'$ and  $w''$. The values $I'_t$, $I''_t$,
$s$ and $v$, can all be computed from image derivatives and the
coordinates of the point in Image~1. 
Given the solution to the ego-motion parameters (to be described
later), one can recover the dense depth map  from 
equations~(\ref{eq:horn.weldon}) and~(\ref{eq:horn.weldon2}). 

%A unique solution is obtained when we reduce the motion model further
%to include infinitesimal translational motion with the model introduced
%by Longuett-Higgins \& Prazdny \cite{Higgins-Prazdny80}:

%\begin{eqnarray}\label{eq:u.v.lh}
%u'=\frac{1}{z}(t'_{1}-x t'_{3}) - w'_{3}y + w'_{2}(1+x^{2}) 
%    - w'_{1} x y\\ 
%v'=\frac{1}{z}(t'_{2}-yt'_{3}) + w'_{3}x - w'_{1}(1+y^{2}) 
%    + w'_{2} x y \nonumber
%\end{eqnarray}
%from which we obtain the 15-parameter model-based brightness constraint:
%\begin{equation}
%I''_{t}s^{T}t' - I'_{t}s^{T}t'' + s^{T}[t'w''^{T} - t''w'^{T}]v = 0
%\label{eq:the.key.eq}
%\end{equation}
%where $s,v$ are defined below:
%\be
%s=\left(\begin{array}{c}
%I_x\\
%I_y\\
%-xI_x - yI_y
%\end{array}\right)
%\label{eq:s}
%\ee
%and
%\be 
%v=\left(\begin{array}{c}
%-I_y - y(xI_x + yI_y)\\
%I_x  +x(xI_x + yI_y)\\
%xI_y - yI_x
%\end{array}\right)
%\ee
%Given the solution to the ego-motion parameters (to be described
%later), one can recover the dense depth map  from the following
%equations \cite{HornWeldon:88}:

%\begin{equation}\label{eq:horn.weldon}
%k s^\top t' + v^\top w' + I'_t = 0
%\end{equation}
%\begin{equation}\label{eq:horn.weldon2}
%k s^\top t'' + v^\top w'' + I''_t=0
%\end{equation}
%where $K=\frac{1}{z}$ denotes the inverse depth at each pixel location.
%Finally, the pure translational small-motion model takes the simplest
%form:
%\begin{equation}
%I''_{t}s^{T}t' - I'_{t}s^{T}t''  = 0
%\end{equation}

Before showing how to solve equation~(\ref{eq:the.key.eq}) for both
rotation and translation we will start with the simple case of pure
translation. This case
arises in practice when we can either assure pure translation or we
have previously rectified the images (perhaps by registering a common
plane in the images as in
~\cite{Kumar-Anandan94-icpr,Irani-et-al94,Sawhney94-cvpr}). It is also
an important case from a theoretical point of view. We will show that
there are problems in recovering the motion parameters in the case of
collinear motion.  If collinear motion causes a problem for the pure
translation case, it will also create a problem for the more general
translation and rotation case.

\subsection{Pure Translation}\label{sec:TBC.pure.trans}


The pure translational small-motion model takes the simplest
form:
\begin{equation}
I''_{t}s^{T}t' - I'_{t}s^{T}t''  = 0
\end{equation}
We have one such equation for each image point and we can write it out in
the matrix form:
$$
{\bf A}t=0
$$
where:
$$
  t=\left(
  \begin{array}{cccccc}   
    t_{x}' & t_{y}' & t_{z}' & t_{x}'' & t_{y}'' & t_{z}''
  \end{array} 
\right)^{T}
$$
and ${\bf A}$ is an $N \times 6$ matrix with the $i'th$ row (corresponding to the
$i'th$ pixel) given by:
$$
  \left(
    \begin{array}{cccccc}
      I''_{t}s_{i1}  & I''_{t}s_{i2}  & I''_{t}s_{i3} 
      & I'_{t}s_{i1}  & I'_{t}s_{i2}  & I'_{t}s_{i3}
    \end{array} 
  \right)
$$
To avoid the trivial solution $t=0$ we add the
constraint $\|t\|=1$.  The least squares problem now maps to the problem of
finding  $\|t\|=1$ that minimizes: 
$$
t^{T}{\bf A}^{T}{\bf A}t=0
$$
The solution is the eigenvector of ${\bf A}^{T}{\bf A}$ corresponding
to the smallest 
eigenvalue. 

\subsubsection{The Singularity of Collinear Motion}\label{sec:TBC.singular.case}

This method fails when the two motions are in the same (or opposite)
directions. In the pure translation case
\be
k s^\top t' + I'_t=0 \label{eq:TBC.horn.weldon.pt}
\ee
and
\be
k s^\top t'' + I''_t=0   \label{eq:TBC.horn.weldon2.pt}
\ee
If the second translation vector $t''$ is proportional to the
translation vector $t'$ then equation~(\ref{eq:TBC.horn.weldon2.pt})
is simply a scaled version of equation~(\ref{eq:TBC.horn.weldon.pt})
adding no new information. The solution will therefore be ill
conditioned when $t' \approx \beta t''$ for a scalar $\beta$.  This is a
drawback in many applications (e.g. 3D reconstruction from a monocular
image sequence).  We present in Section~\ref{sec:TBC.future.work} some
early research on ways to overcome this problem. The initial results
look promising.

\subsection{Solving for Translation and Rotation}\label{sec:TBC.bilinear.solv}

%\subsubsection{show that A is rank 13 not 14}

In equation~(\ref{eq:the.key.eq}) the ego-motion parameters  are embedded in
a 15 parameter model: 6  translation parameters and the 9 outer-product
terms $[t'w''^{T} - t''w'^{T}]$. Each pixel provides one linear
equation, thus N pixels provide a system ${\bf A}c=0$ where ${\bf A}$ is an
$N\times 15$ measurement matrix and $c$ is the 15-vector of unknown
parameters. The least squares solution is the eigenvector of ${\bf
A}^\top {\bf A}$
corresponding to the smallest eigenvalue. Without noise the eigenvalue
will be zero and the solution is in the null space of  ${\bf A}^\top
{\bf A}$. The
solution is unique if the null space is one dimensional and the 
rank of ${\bf A}$ is 14. We will prove however, that the null space
is at least two dimensional. 

\begin{proposition}%{ The rank of {\bf A} is 13 not 14.}
Let ${\bf A}$ be the $N\times 15$ measurement matrix associated with the homogeneous
equation~(\ref{eq:the.key.eq}),  where $N \gg 15$. Then the null space
of the matrix ${\bf A}^\top {\bf A}$ is of dimension greater or equal to 2 and
$Rank({\bf A}) \leq 13$. 
\end{proposition}
{\bf Proof:}

\noindent
Since
$$
s^\top v = s^\top(p\times s) = 0
$$
the vector 
$$
c_{0} = (0,0,0,0,0,0,1,0,0,0,1,0,0,0,1)^{T}
$$
is in the null space of ${\bf A}$. 

\noindent
The vector $c_0$ corresponds to:
$$
\begin{array}{c}
t'=0 \\
t''=0 \\
\left[t'w''^{T} - t''w'^{T}\right] =   
  \left[
    \begin{array}{ccc}
      1 & 0 & 0\\
      0 & 1 & 0\\
      0 & 0 & 1
    \end{array} 
  \right]
\end{array}
$$
Therefore $c_0$ is clearly not an admissible solution,
and therefore the null space of ${\bf A}$ includes two vectors: the vector $c_0$
and the true solution. It follows that the rank of ${\bf A}$ is at most 13. 
\eop

The method of solution is based on the following observation. Let
$b_0$ be another vector of the null-space (for example, let $c_0$ and
$b_0$ be the two eigenvectors of ${\bf A}^\top {\bf A}$ corresponding
to the two smallest eigenvalues), then the desired solution vector $b$
is a linear combination of the two:
\be
b = b_{0} + \alpha c_{0}
\ee
We note that $c_0$ is in the null space of ${\bf A}$ regardless of noise in
the data and will therefore correspond to an eigenvalue of zero, up to
numerical precision. Theoretically, $b_0$ is also in the null space of
${\bf A}$ but only in the case of noise free data and exact motion
model. Therefore $b_0$ is the eigenvector corresponding to the second
smallest eigenvalue.  We next show how to determine $\alpha$.

\subsubsection{Solution Using the Rank=2 Constraint.}

In order to find $\alpha$ given $c_{0}$ and
$b_{0}$ we enforce the constraint that:
$$Rank[t'w''^{T} - t''w'^{T}]=2$$.
Clearly the choice of $\alpha$ will have no affect the first 6 elements of
the vector $b$. Let us arrange the last 9 elements of $b$, $b_{0}$ and
$c_{0}$ into the corresponding $3 \times 3$ matrices ${\bf B}$, ${\bf B_{0}}$ and
${\bf C_{0}}$. We are now looking for an $\alpha$ such that: 
\begin{equation}\label{eq:rank=2}
Rank({\bf B_{0}}-\alpha {\bf C_{0}}) = 2.
\end{equation}
In our case ${\bf C_{0}}$ is the identity matrix so~(\ref{eq:rank=2})
becomes:
$$
Rank({\bf B_{0}}-\alpha {\bf I}) = 2.
$$
and the solution
for $\alpha$ is given simply by the eigenvalues of ${\bf B_{0}}$. Since ${\bf
B_{0}}$ is a $3 \times 3$ matrix this results in up to three discrete
solutions. We now prove that only one of those can be a valid
solution.

\begin{theorem}{The 15-parameter model has a unique
solution.}\label{thm:LH.uniqueness}
\noindent
Let:
\be
{\bf B} = \left[ t' w''^{T} - t''w'^{T} \right]. 
\ee
Let:
\be
\tilde{\bf B} = \left[ t' w''^{T} - t''w'^{T} \right] + \alpha {\bf I}. 
\ee
Then the equation:
\be\label{eq:LH.uniqueness}
\tilde{\bf B} = \left[ t' \tilde{w}''^{T} - t''\tilde{w}'^{T} \right]
\ee
has a solution for $\tilde{w}'',\tilde{w}'$ given $t',t''$ iff $\alpha = 0$.
\end{theorem}
{\bf Proof:\ } 

The {\em if} part of the proof is trivial. For $\alpha=0$,
$\tilde{w}''=w''$ and $\tilde{w}'=w'$ is a solution for (\ref{eq:LH.uniqueness}).
 To prove the {\em only if} part, suppose there exist
$\tilde{w}'',\tilde{w}'$ that is a solution for (\ref{eq:LH.uniqueness}). 
Then:
\be\label{eq:LH.uniqueness.2}
\left[t'w''^{T} - t''w'^{T} \right] + \alpha {\bf I} = 
\left[ t' \tilde{w}''^{T} - t''\tilde{w}'^{T} \right].
\ee 
Let $q$ be a vector such that $q^{T}t'=0$,
$q^{T}t''=0$ and $q \not= \vec{0}$. Multiply left and right side of (\ref{eq:LH.uniqueness.2}) by
$q^{T}$
to get: 
\be\label{eq:LH.uniqueness.3}
q^{T}\left[t'w''^{T}\right] - q^{T}\left[t''w'^{T} \right] + q^{T}\alpha {\bf I} = 
q^{T}\left[ t' \tilde{w}''^{T}\right]  - q^{T}\left[t''\tilde{w}'^{T} \right].
\ee
Since $q^{T}t'=0$ and $q^{T}t''=0$ most of terms in
(\ref{eq:LH.uniqueness.3}) become zero and we are left with:
\be
q^{T}\alpha {\bf I} = \vec{0}. 
\ee
Since $q \not= \vec{0}$ then $\alpha = 0$.
\eop
 

\subsubsection{The Algorithm for Finding Motion Parameters}\label{sec:TBC.algorithm}
Based on the previous arguments the algorithm for finding the motion
parameters is as follows: 

\begin{enumerate}
\item
Compute the  $N \times 15$ matrix ${\bf A}$ from  equation~(\ref{eq:the.key.eq}). 
\item
Find the eigenvector corresponding to the second smallest eigenvalue of  the
matrix ${\bf A}^{T}{\bf A}$. This is  $b_{0}$. (The vector $c_{0}$
corresponds to the smallest.)
\item
The first 6 elements of $b_0$ are the translations, $t'$,$t''$. 
\item
Arrange the last 9 elements of  $b_0$ into a $3 \times 3$ matrix
${\bf B_0}$.
\item
Construct three possible rank 2 matrices ${\bf B}_{i, i=1..3}$ from
${\bf B}_{0}$:
\begin{equation}
{\bf B}_i = {\bf B}_{0} - \alpha_i {\bf I}
\end{equation}
where $\alpha_{i, i=1..3}$ are the three eigenvalues of  ${\bf B}_{0}$.
\item
For each $\alpha_{i, i=1..3}$
solve:
$$
[t'w''^{T} - t''w'^{T}] = {\bf B},
$$ 
\newline
for $w'$ and $w''$, given  $t'$ and $t''$ from step (3). This a set of
9 linear equations in the six unknowns ( $w'$,$w''$) and is solved
using least squares.
\item
From Theorem~\ref{thm:LH.uniqueness} only one of the three solutions
is correct. We select the solution that best fits the data.
For each of the three solutions use the 12 parameters $t'$, $t''$,
$w'$ and $w''$ to form the 15 intermediate parameter vector $b_i$ and then
compute the error: 
\be
E_i = b_i^T {\bf A}^T {\bf A}  b_i.
\ee 
Of the three solutions select the solution which gives the
smallest error. 
\end{enumerate}

\subsubsection{Solve as a Nonlinear Optimization}\label{sec:TBC.nonlinear.optim}

In Section~\ref{sec:TBC.algorithm} we use linear methods to solve
for the motion parameters. In order to do so we have treated the 15
intermediate parameters as linearly independent parameters while they
are in fact bilinear combinations of 12 independent parameters. It is
often possible to improve the motion estimates by using the linear
solution only as a starting point for non-linear optimization techniques
\cite{Faugeras-Papa-iccv98}.  In our experiments we did not find that
a nonlinear optimization stage improved the results. This is possibly
because the iterative framework
(Section~\ref{sec:iterative.refinement}) is in itself a form of
nonlinear optimization.

The nonlinear optimization procedure is brought here for completeness. 
Our general problem is to find motion parameters 
$t'$, $t''$, $w'$ and 
$w''$ that
minimize the cost function:
\be\label{eq:TBC.cost.function}
E(t',t'',w',w'') = b^\top {\bf A}^\top {\bf A}b
\ee
where {\bf A} is the estimation matrix derived from
equation~(\ref{eq:the.key.eq}) and $b(t',t'',w',w'')$ is the fifteen
element vector of bilinear functions of the motion parameters.


One method would be to find motion parameters $t'$, $t''$, $w'$ and 
$w''$ that
minimize equation~(\ref{eq:TBC.cost.function}). In order to avoid the
trivial solution we need to add a constraint such as: $|t'|^2 +
|t''|^2 = 1$. This becomes a 12 parameter nonlinear constrained
optimization problem.

The process can be simplified by noting that equation~(\ref{eq:the.key.eq})
can be rewritten in the form:
\be
\left(I''_t-w''^\top v''\right){s'}^\top t' - 
\left(I'_t-w'^\top v'\right){s''}^\top t'' =0.\nonumber
\ee
given a set of rotation values $w'$ and $w''$, one can compute the
least squares estimate of $t'$ and $t''$ using linear methods and also
compute the least squares error (equation~\ref{eq:TBC.cost.function}). One
can then use nonlinear search techniques to find $w'$ and $w''$ which
minimize this least squares error. This becomes a six parameter
unconstrained optimization problem.


%\section{Solving the Ego-motion Parameters}
%\label{sec:bilinear.solv}

%The ego-motion parameters from eqn.~\ref{eq:the.key.eq} are embedded in
%a 15 parameter model: 6  translation parameters and the 9 outer-product
%terms $[t'w''^{T} - t''w'^{T}]$. Each pixel provides one linear
%equation, thus N pixels provide a system $Ac=0$ where $A$ is an
%$N\times 15$ measurement matrix and $c$ is the 15-vector of unknown
%parameters. The least squares solution is the eigenvector of $A^\top A$
%corresponding to the smallest eigenvalue. The solution is unique if the
%rank of $A$ is 14. However, since $s^\top v = s^\top(p\times s) = 0$
%the vector 
%\begin{equation}
%c_{0} = (0,0,0,0,0,0,1,0,0,0,1,0,0,0,1)^{T}\nonumber
%\end{equation}
%is in the null space of $A$. Since $c_0$ is not an admissible solution,
%the rank of $A$ is 13, and the null space of $A$ is two-dimensional.
%Let $b_0$ be another vector of the null-space (say $c_0$ and $b_0$ are
%the two eigenvectors of $A^\top A$ corresponding to the two smallest
%eigenvalues), then  the
%desired solution vector $b$ is a linear combination of the two:
%\begin{equation}
%b = b_{0} + \alpha c_{0}
%\end{equation}
%In order to find $\alpha$ given $c_{0}$ and $b_{0}$ we  enforce the
%constraint that the matrix $[t'w''^{T} 
%- t''w'^{T}]$ is of  rank 2.

%Clearly the choice of $\alpha$ will have no affect the first 6 elements of
%the vector $b$. Let us arrange the last 9 elements of $b$, $b_{0}$ and
%$c_{0}$ into the corresponding $3 \times 3$ matrices $B$, $B_{0}$ and
%$C_{0}$. We are now looking for an $\alpha$ such that: 
%\begin{equation}
%Rank(B_{0}-\alpha C_{0}) = 2
%\end{equation}
%Since $C_{0}$ in our case is the identity matrix, the solution for $\alpha$
%is given by the eigenvalues of $B_{0}$. We choose the one with the smallest
%absolute value. 
%% The vector $b_{0}$ is the eigenvector corresponding to the
%% second smallest eigenvalue of the matrix $A^{T}A$. (The vector $c_{0}$
%% corresponds to the smallest.)

%Based on the previous arguments the algorithm for finding the motion
%parameters is as follows: 

%\begin{enumerate}
%\item
%Compute the  $N \times 15$ matrix $A$ from  equation (\ref{eq:the.key.eq}). 
%\item
%Find the eigenvector corresponding to the second smallest eigenvalue of  the
%matrix $A^{T}A$. This is  $b_{0}$. (The vector $c_{0}$
%corresponds to the smallest.)
%\item
%The first 6 elements of $b_{0}$ are the translations, $t'$,$t''$. 
%\item
%Arrange the last 9 elements of  $b_{0}$ into a $3x3$ matrix $B_{0}$.
%\item
%Construct a rank 2 matrix $B$ from  $B_{0}$:
%\begin{equation}
%B = B_{0} - \alpha I
%\end{equation}
%\item
%Solve:\hspace{0.25in}
%$
%[t'w''^{T} - t''w'^{T}] = B
%$, 
%\newline
%for $w'$ and $w''$, given  $t'$ and $t''$ from step (3).
%\end{enumerate}
 
%%\subsection{The singular case: }
%%\label{sec:singular.case}

%%This method fails when the two motions are in the same (or opposite)
%%directions. This is obvious in the pure translation case. If the second
%%translation vector $t''$ is proportional to the translation vector $t'$
%%then equation (\ref{eq:horn.weldon2}) is simply a scaled version of
%%equation (\ref{eq:horn.weldon}) adding no new information. Under current
%%research is a way to overcome this problem (section \ref{sec:future.work})
%%and initial results look promising.

\section{Implementation Details}\label{sec:details}

%In this section we will provide the details needed to
%build a program that actually works. 

\subsection{Computing the Depth, Smoothing and Interpolation}\label{sec:TBC.compute.depth}

After recovering the camera motion (Section~\ref{sec:TBC.algorithm}) we use equations~(\ref{eq:horn.weldon}) and~(\ref{eq:horn.weldon2}) to compute depth at every
point. 
Information is combined from both image pairs by minimizing the least
squares error: 
$$
E = \begin{array}{c}
\\
\min \arg\\
k
\end{array}  \sum_{j=1}^{2} \left(
ks^{T}t^{j} + v^{T}w^{j} + I^{j}_{t} \right)^{2}. 
$$ 
Here $j=1$ and $j=2$ denote values from the first and second image pairs
respectively. 

There are points in the image
where the brightness gradients are close to zero, (and therefore
$s^{T}t^{j} \simeq 0$)  and the
estimation of $k$ will be ill-conditioned at those points. In order to
overcome this problem we use a local region of support around the
point of interest. It is simplest form we assume the depth is
constant in the region and minimize:
\be\label{eq:TBC.depth}
E = \begin{array}{c}
\\
\min \arg\\
k
\end{array}
 \sum_{x,y \in R} \sum_{j} \beta (x,y) |s^{T}t^{j}|^{p} \left(
ks^{T}t^{j} + v^{T}w^{j} + I^{j}_{t} \right)^{2} \nonumber
\ee
where the windowing function $ \beta (x,y)$ allows one to increase the weights of the
points closer to the center of the region.
 The $|s^{T}t^{j}|^{p}$
term reduces the weight of points which have a small gradient or where the
gradient is perpendicular to the camera motion since these cases are
highly affected by noise. We used $p=1$. 

During the iteration process we typically used a region ${R}$ of
$7 \times 7$ to $11 \times 11$.  For `prettier' results the last
iteration, we typically reduced the region ${R}$ to $1 \times 1$
but added a very weak global smoothness term and performed multi-grid
membrane interpolation\cite{Numerical-recipes}. The smoothness term stabilizes regions where
there is no image gradient so very small regions of support
can be used.



\subsection{Iterative Refinement }\label{sec:iterative.refinement}

The constant brightness constraint is a linearized form of the Sum
Square Difference (SSD) criteria. The linear solution can be thought
of as a single iteration of Newton's method applied to the problem. Iterative
refinement is performed as follows:
\begin{enumerate}
\item
Calculate motion (using equation~\ref{eq:the.key.eq}).
\item
Compute depth (using equation~\ref{eq:TBC.depth}).
\item
Using the depth and motion, warp Images~2 and~3 towards Image~1 (Section~\ref{sec:TBC.image.warping}). 
\item
Compute new time derivatives $I'_t$ and $I''_t$.
\item
Compute a new motion and depth estimate.
\end{enumerate}

In the ideal case, as the final result, the warped images should
appear nearly identical to Image~1.  One must be careful and not 
simply compute the incremental changes in the translation $\delta t'$
and $\delta t''$ . As the images are warped closer together and the
translation estimate approaches zero, the system of equations~(\ref{eq:the.key.eq}) will become ill conditioned. Furthermore, since
the equations are homogeneous we must enforce a constraint such as
$|t'|^2 + |t''|^2 = 1$ to avoid the trivial solution. We would not
wish to apply such a constraint to $\delta t'$ and $\delta t''$. Thus,
one must compute the full translation model (previous iteration plus
the incremental change to the translation).  This problem does not
arise for the rotations: $w'$ and $w''$. It is in fact convenient to
warp the images using the best rotation estimate and then compute only
the incremental improvement in the rotations: $\delta w'$ and $\delta
w''$. For example we might chose to warp the images using a more exact
rotation model than the linearized small rotation model.

The method for computing the full translation model and incremental
rotation is now described. 
Let $\Psi_{1}$, $\Psi_{2}$ and $\Psi_{3}$ be the three images. Assume we
have $\hat{k}$, $\hat{t}^{j}$ , $\hat{w}^{j}$, from the previous
iteration. (Here again we use $j$ to stand for either $'$ or $''$.)  The translation components of image motion $\hat{u}_t^{j}$
and $\hat{v}_t^{j}$, and the rotational components  $\hat{u}_r^{j}$
and $\hat{v}_r^{j}$  can be
computed using $\hat{k}$, $\hat{t}^{j}$ and $\hat{w}^{j}$
in equation~(\ref{eq:TBC.u.v.lh}).
These are then used  to warp
images $\Psi_{2}$ to $\hat{\Psi}_{2}$ and $\Psi_{3}$ to $\hat{\Psi}_{3}$.
After warping, the images satisfy the brightness constraint equation:

\begin{eqnarray}\label{eq:modified.CBE}
I_{x}du' + I_{y}dv' + \hat{I}'_{t} =0\\
I_{x}du'' + I_{y}dv'' + \hat{I}''_{t} =0 \nonumber
\end{eqnarray}
where the temporal derivatives at each pixel are given by:
\begin{eqnarray}
\hat{I}'_{t} = \hat{\Psi}_{2} - \Psi_{1}\\
\hat{I}''_{t} = \hat{\Psi}_{3} - \Psi_{1}\nonumber
\end{eqnarray}
and $du^{j}$, $dv^{j}$ are the (still unknown) differences between computed
image motions and the real image motions:
\begin{eqnarray}\label{eq:du.hat}
du^{j} = u^{j} - \hat{u}_t^{j} - \hat{u}_r^{j} \\
dv^{j} = v^{j} - \hat{v}_t^{j} - \hat{v}_r^{j} \nonumber
\end{eqnarray}
Let:
\begin{equation}\label{eq:alpha}
\alpha^{j} = I_{x}\hat{u}_t^{j} + I_{y}\hat{v}_t^{j} 
\end{equation}
which can also be written as:
\begin{equation}
\alpha^{j} = \hat{k}s^{T}\hat{t}^{j} \nonumber
%\alpha^{j} = \hat{k}s^{T}\hat{t}^{j} + v^{T}\hat{w}^{j}  \nonumber
\end{equation}
Substituting equations~(\ref{eq:alpha}) and~(\ref{eq:du.hat}) in equation
(\ref{eq:modified.CBE}) we get:
\begin{eqnarray}\label{eq:modified.CBE.2}
I_{x} (u' - \hat{u}_r^{'}) + I_{y} (v' - \hat{v}_r^{'}) + (\hat{I}'_{t} - \alpha') =0 \\
I_{x} (u'' - \hat{u}_r^{''}) + I_{y} (v'' - \hat{v}_r^{''}) + (\hat{I}''_{t} - \alpha'') =0 \nonumber
\end{eqnarray}
Substituting equation~(\ref{eq:TBC.u.v.lh}) in equation
(\ref{eq:modified.CBE.2}) we get modified versions of the
equations~(\ref{eq:horn.weldon}) and~(\ref{eq:horn.weldon2})
\begin{eqnarray}
ks^{T}t' + v^{T} \delta w' + (I'_{t} - \alpha')= 0 \\
ks^{T}t'' + v^{T} \delta w'' + (I''_{t} - \alpha'')= 0 \nonumber
\end{eqnarray}
We start our first iteration with  $\hat{k}$, $\hat{t}^{j}$ , $\hat{w}^{j}$
all zero and therefore $\alpha=0$ as well. 

\subsection{Image Warping}\label{sec:TBC.image.warping}

Given an estimate of the camera motion and the depth at every point, we
can warp Image~2 towards Image~1. Image warping in general is
described in \cite{Wolberg90}.  Let $I'_{old}$ be the original Image~2
and let $I'_{new}$ be the warped image we are trying to create. We
must first define functions $x_{old}(x,y,...)$ and $y_{old}(x,y,...)$
which, given image coordinates in the new image $(x,y)$ and possibly
some extra parameters (symbolized by $,...$), return the coordinates of
the corresponding point in the old image ($x_{old},y_{old}$).

In our case the {\em extra parameters} are the camera motions $t'$ and
$w'$, and the depth map at every point $k(x,y)$. The functions
$x_{old}$ and $y_{old}$ depend on the motion model. In our
implementation we used the LH\&P small motion model
(equations~(\ref{eq:TBC.u.v.lh})). Therefore:
\be
\begin{array}{lcl}
x_{old} & =& x+u' \\
& =& x + \frac{1}{z}(t'_{1}-x t'_{3}) - w'_{3}y + w'_{2}(1+x^{2}) 
    - w'_{1} x y
\end{array}
\ee
and
\be
\begin{array}{lcl}

y_{old}  & =& y+v' \\
& =& y + \frac{1}{z}(t'_{2}-y t'_{3}) + w'_{3}x - w'_{1}(1+y^{2}) 
    + w'_{2} x y.
\end{array}
\ee


After defining $x_{old}$ and $y_{old}$ we can compute the value for
every pixel in $I'_{new}$ according to the formula:
\be
I'_{new}(x,y) = \left\{ 
\begin{array}{lr}
\lefteqn{I'_{old}(x_{old}(x,y,...),y_{old}(x,y,...)),} &  \\
    & if (x_{old},y_{old}) \in I'_{old} \\
0, &   otherwise.
\end{array}
\right.
\ee
There are two points to note:
\begin{itemize}
\item
 The values of $x_{old}$ and $y_{old}$ are in general non-integer, so
we use bilinear interpolation to compute the appropriate pixel
values. Bilinear interpolation works better than nearest
neighbor. The use of the more complex bicubic interpolation resulted in no
noticeable improvement over bilinear interpolation.
\item
If the coordinates ($x_{old}$,$y_{old}$) are outside the
coordinates of the image then we use the value $0$. These points
should be marked as invalid and not used in subsequent 
motion and depth computation. 
\end{itemize}

\subsection{Coarse-to-Fine Processing}\label{sec:TBC.coarse.to.fine}

In order to deal with image motions larger than 1 pixel we use
a Gaussian pyramid for coarse to fine processing
\cite{Bergen:92,Burt:83}.
Each of the three images is filtered by an approximation to a Gaussian
filter and subsampled to create an image of half the size (in each
dimension). This operation is performed recursively to create
a pyramid with four or five levels depending on the size of the
original image.  The following $5 \times 5$ kernel was used as the filter:
$$
\frac{1}{256} \times \left[
\begin{array}{ccccc}
1 & 4 & 6 & 4 & 1 \\
4 & 16 & 24 & 16 & 4\\
6 & 24 & 36 & 24 & 6 \\
4 & 16 & 24 & 16 & 4 \\
1 & 4 & 6 & 4 & 1 
\end{array}
\right].
$$
This filter is separable and can be implemented by convolving the
rows and columns of the image, twice each, with the filter
$f=(0.25,~0.5,~0.25)$. 

Starting from the coarsest level of we perform a few iterations to
compute motion and depth according to the scheme described in Section~\ref{sec:iterative.refinement}. The brightness derivatives
($I_x,I_y,I'_t,I''_t$) are computed using the subsampled images. After
computing the motion and depth, we create a finer depth map
from the coarse depth map using bilinear interpolation. We then use the
interpolated depth map, and the motion estimates, as the
starting values for the iterations at the finer level.  

\subsection{Field of View}\label{sec:TBC.FOV}

In this section we will first demonstrate how in the general problem
of determining camera motion from motion fields, when the field of
view is narrow, there exists an ambiguity between rotation and
translation. We will also note the errors in the estimation of the
rotation if we use the wrong focal length value. Then we show how this
ambiguity is observed in the shape of the cost function
(equation~\ref{eq:TBC.cost.function}) from three views. Finally we
describe how to achieve stable and accurate motion estimates in the
presence of rotation even without very wide fields of view.

\subsubsection{Rotation-Translation Ambiguity from Motion Fields}

It is well known that for a medium to narrow field of view, the motion field due
to camera rotation around the $Y$ axis, can be indistinguishable from the
motion field due to some translation along the $X$ axis and an
appropriate depth surface. A similar ambiguity exists
between rotations around the X axis and translations along the $Y$ axis.

In Figure~\ref{fig:TBC.rot.trans.ambig.1} we see that for
medium-narrow fields of view ($30^o$) the flow due to rotation around
the $Y$ axis is parallel to the $X$ axis. The length of the flow vectors
due to rotation is slightly shorter, closer to the $Y$ axis. A similar
flow field can be produced by camera translation along the $X$ axis with
a surface that is curved so that points in the center of the image are
more distant. There are motion flows that unambiguously indicate a
translational component to the motion: motions fields that include a
focus-of-expansion or show parallax effects. Nevertheless even these
fields are ambiguous since a rotational component could be added to
the motion, compensated for by a change in the translation and depth
to produce the same motion field. In some cases the resulting depth map
will include negative depths indicating some obvious error. This is
often formulated as the {\em depth positive constraint}\cite{Fermuller-Aloimonos-97} which can be used to
limit the range of motion ambiguity.

\begin{figure}
  \centerline{
    \epsfxsize=3.2in \epsffile[018 180 594 612]{flow-rot.ps}}

  \caption[Motion flow field due to rotation around the $Y$ axis]{\em Motion flow field due to rotation around the $Y$ axis for
wide, medium and narrow fields of view. The arrows for the medium and
wide have been offset in the Y direction for clarity. For the narrow
field of view ($f=800_{pixels}$, FOV=$30^o$ ) the flow field is large and
almost parallel to the $X$ axis. For wider fields of view ($f=400$
and  $f=200$) the $X$ component of the flow is smaller but there is a
noticeable $Y$ component. }
  \label{fig:TBC.rot.trans.ambig.1}
\end{figure}

With wider fields of view (i.e. shorter focal lengths) the ambiguity
disappears. The motion fields due to rotation cannot be modeled by
translational motion. Figure~\ref{fig:TBC.rot.trans.ambig.1} also
shows the motion fields due to rotation with fields of view of $53^o$
and $90^o$. The distribution of the $Y$ component of the
flow has the unique characteristics of rotational flow. For wide fields of view
the $Y$ component points towards the $X$ axis on the right of the image
and away from the X axis on the left of the image (for this particular
sign of the rotatation). 

The motion field for a wide field of view is unambiguous if we know
the focal length, or equivalently, the field of view (i.e. the
camera's internal parameters are known). If our estimate of the focal
length is incorrect then the motion and structure estimates which best
fit the flow field will be incorrect. If the focal length estimate is
smaller than the true focal length then the estimated angle of rotation
will be smaller than the true rotation angle (in absolute
values). Figure~\ref{fig:TBC.rot.trans.ambig.2}a shows the $X$ and $Y$
components of the motion field of a rotation of $-4.5^o$ around the $Y$
axis using a focal length of $400_{pixels}$. Figure~\ref{fig:TBC.rot.trans.ambig.2}b shows the same measurements for a
focal length of $200_{pixels}$ and a rotation of $-2.5^o$. The $Y$ components are
similar. The $X$ components are very different but these can be
accounted for by adjusting the translation and depth estimates.

\begin{figure}
  \centerline{
    \epsfxsize=3.2in \epsffile[018 180 594 612]{flow-rot-Y.ps}}
  \centerline{ (a)}
  \centerline{
    \epsfxsize=3.2in \epsffile[018 180 594 612]{flow-rot-X.ps}}
  \centerline{ (b)}
  \caption[The $X$ and $Y$ components of the motion field due a rotation
around the $Y$ axis]{\em 
    The $X$ and $Y$ components of the motion field due a rotation
around the $Y$ axis of $-4.5^o$ with a focal length of $400_{pixels}$
for a rotation of $-2.5^o$ with a focal length of $200_{pixels}$. The
arrows for the $200_{pixels}$ camera have been plotted offset right
(a) and up (b). Note that the $Y$ flow components are very similar. }
  \label{fig:TBC.rot.trans.ambig.2}
\end{figure}

\subsubsection{Rotation-Translation Ambiguity in the 15-parameter Model}

The rotation translation ambiguity can also be observed while using
the 15-parameter model. Figure~\ref{fig:TBC.rot.trans.ambig.3}
show the shape of the cost function (equation~\ref{eq:TBC.cost.function}) near the global minimum. These are
simulation results for focal lengths of $50$, $100$ and $200$. The image size
was $320 \times  240$. It is not possible to plot the full 12
dimensional surface so the figures show only the shape of the function
as we vary two motion parameters at a time. We show the effect of
varying $w'_x$ and $t'_y$ where we expect to see ambiguity and  the
effect of varying $w'_x$ and $t'_x$ where no ambiguity is expected. 
In Figure~\ref{fig:TBC.rot.trans.ambig.3}f, the long diagonal
valley shape of the cost function indicates the ambiguity of the
motion estimate for narrower fields of view. 

\begin{figure}
  \centerline{
    \hfill
    \epsfxsize=1.0in \epsffile[018 180 594 612]{costFunc-txwx50.ps}
    \hfill
    \epsfxsize=1.0in \epsffile[018 180 594 612]{costFunc-txwx100.ps}
    \hfill
    \epsfxsize=1.0in \epsffile[018 180 594 612]{costFunc-txwx200.ps}
    \hfill
  }
  \centerline{\hfill (a)  \hfill\hfill \hfill (b) \hfill \hfill  \hfill(c) \hfill}
  \centerline{
    \hfill
    \epsfxsize=1.0in \epsffile[018 180 594 612]{costFunc-txwy50.ps}
    \hfill
    \epsfxsize=1.0in \epsffile[018 180 594 612]{costFunc-txwy100.ps}
    \hfill
    \epsfxsize=1.0in \epsffile[018 180 594 612]{costFunc-txwy200.ps}
    \hfill
  }
  \centerline{\hfill (d)  \hfill\hfill  \hfill(e) \hfill  \hfill \hfill(f) \hfill}
  \caption[The shape of the cost function for various
fields of view]{\em 
The shape of the cost function around the minimum for various
fields of view.  These simulation results were computed after the
program converge to the final motion and depth estimates. For
visualization only two of the 12 parameters were changed at one
time. Note that for a narrower field of view (f) the cost function has
the shape of a narrow diagonal valley. This indicates an ambiguity
between translation along the $X$ axis and rotations around the $Y$
axis. }
  \label{fig:TBC.rot.trans.ambig.3}
\end{figure}

\subsubsection{The Solution }\label{sec:TBC.lie.to.system}

As we have seen, in order to get reasonable rotation and translation
estimates a wide field of view is required. The theoretical results of
\cite{Michaels-thesis92} and simulation results indicate that very wide
fields of view are required ($120^o$ or greater). These are hard to
achieve with standard lenses. With narrower fields of view the
iterative scheme described in Section~\ref{sec:iterative.refinement}
fails to converge. In previously published work
\cite{Stein-Shashua-cvpr97} we managed to stabilize
the results for small rotation angles by unknowingly biasing the
results towards small rotation. This was done by setting the focal
length parameter to $50$ instead of around $600$. For a $640 \times
480$ size image this implied a field of view of $150^o$ instead of the
correct $55^o$. As a result the rotation estimates were considerably
smaller that the true values. This also resulted in errors in the
translation and depth estimates.

It is important to note that while the magnitude of the angle
estimates in \cite{Stein-Shashua-cvpr97} were too small, the estimated
axis of rotation was correct. In fact the results show a linear
relationship between true angle and estimated angle.  That the
rotation values were always smaller than the true values is consistent
with the simulation results shown in
Figure~\ref{fig:TBC.rot.trans.ambig.2}. This leads us to the
following simple modification to the iterative scheme (Section~\ref{sec:iterative.refinement}):
\begin{enumerate} 
\item
For the estimation of motion use a focal length parameter $f <<
f_{true}$. 
\item
During the image warping stage use the correct focal length to rotate
Images~2 and Image~3 towards Image~1. 
\end{enumerate}

This modification required no changes to the structure of the
algorithm since the coarse-to-fine and iterative mechanisms were
already in place. The results shown in Section~\ref{sec:TBC.exp.cvpr97} show that the modified algorithm
gives accurate rotation estimates even for medium fields of view
($55^o$). 

\subsubsection{How does this affect the translation estimates?}

How does using the wrong focal length affect the translation
estimates ? Looking at equation~(\ref{eq:the.key.eq}) (and also 
(\ref{eq:horn.weldon}), (\ref{eq:horn.weldon2}) and
(\ref{eq:TBC.depth})) we see that the translations $t'$ and $t''$
always appear in a dot product with $s$ where: 
$$
s=\left(\begin{array}{c}
I_x\\
I_y\\
-xI_x - yI_y
\end{array}\right).
$$ 
The first and second terms of $s$ are not affected by the focal
length estimate. Therefore there is no effect on the translation
estimate in the $X$ and $Y$ direction, $t_x$ and $t_y$. In the third
term, $x$ and $y$ are supposed to be the normalized image
coordinates. If we use the wrong focal length then the estimated
$Z$ translation ($t_z$) will be scaled accordingly. 
$$
\frac{t_Z}{ t_{Z true}} =  \frac{f}{f_{true}}. 
$$
Using this simple relationship we can recover the true translation
direction from the estimated one even if we use a focal length value
which is much smaller than the true value. 

\subsection{Final Touches - Recompute the Depth Keeping the Motion Constant}\label{sec:TBC.final.touches}

Initially, during the first iterations and at the coarsest levels of
the pyramid, the motion estimates are not accurate and cause errors in
the depth estimates. Sometimes, particularly near the borders of the
image, these errors are large enough to affect the new derivatives,
$I'_t$ and $I''_t$ in such a way that at the next iteration even when
the motion estimates are more accurate, the depth values are still badly in
error. This tends to be a local problem and does not affect the global
motion estimates but it does create local 'holes' in the depth map. 

To fix this problem, after computing the motion at
the finest level, we go back down the pyramid and recompute the depth
while keeping the motion estimate constant. This is in effect becomes a
three camera stereo computation. 
Figure~\ref{fig:TBC.fix.holes} shows an example of the depth
estimates before and after the second pass. 

\begin{figure}
        \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir seq90iss.050803.W3p0.depth.0.ps}
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir seq90iss.050803.W3.depth.0.ps}
              \hfill }
        \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
  \caption[Local depth error fixed by a second
pass down the pyramid performing `depth-only' computation]{\em An
example of a local depth error (a) fixed by a second
pass down the pyramid performing `depth-only' computation (b). (see text)}
  \label{fig:TBC.fix.holes}
\end{figure}



%\section{Implementation details}\label{sec:details}

%\subsection{Iterative refinement  and coarse to fine processing}

%The constant brightness constraint is a linearized form of the Sum Square
%Difference (SSD) criteria. The linear solution can be thought of as a
%single iteration of Newton's method applied to the problem. Iterative
%refinement is performed as follows: First one calculates motion and depth
%using eqns.~\ref{eq:the.key.eq},\ref{eq:horn.weldon},\ref{eq:horn.weldon2}. 
%Then, using the depth and motion, images $1$ and
%$2$ are warped towards image $0$. A correction to the depth and
%motion is computed using the warped images. In the ideal case, as the final
%result the warped images should appear nearly identical to image $0$.

%One must be careful not to compute the incremental  model from the
%previous iteration because as the iterations proceed the system of
%equations of the incremental model will become badly conditioned.
%Thus, one must compute the full model (previous iteration plus the
%incremental change to the model). This is described below.


%Let $\Psi_{0}$, $\Psi_{1}$ and $\Psi_{2}$ be the three image. Assume we
%have $\hat{k}$, $\hat{t}^{j}$ , $\hat{w}^{j}$, from the previous
%iteration.  The image motions $\hat{u}^{j}$ and $\hat{v}^{j}$ can be
%computed using $\hat{K}$, $\hat{t}^{j}$ and $\hat{w}^{j}$ in 
%equation (\ref{eq:u.v.lh}).
%These are then used  to warp
%images $\Psi_{1}$ to $\hat{\Psi}_{1}$ and $\Psi_{2}$ to $\hat{\Psi}_{2}$.
%After warping, the images satisfy the brightness constraint equation:

%\begin{eqnarray}\label{eq:modified.CBE}
%I_{x}du' + I_{y}dv' + \hat{I}'_{t} =0\\
%I_{x}du'' + I_{y}dv'' + \hat{I}''_{t} =0 \nonumber
%\end{eqnarray}
%where the temporal derivatives at each pixel are given by:
%\begin{eqnarray}
%\hat{I}'_{t} = \hat{\Psi}_{1} - \Psi_{0}\\
%\hat{I}''_{t} = \hat{\Psi}_{2} - \Psi_{0}\nonumber
%\end{eqnarray}
%and $du^{j}$, $dv^{j}$ are the (still unknown) differences between computed
%image motions and the real image motions:
%\begin{eqnarray}\label{eq:du.hat}
%du^{j} = u^{j} - \hat{u}^{j} \\
%dv^{j} = v^{j} - \hat{v}^{j} \nonumber
%\end{eqnarray}
%Let:
%\begin{equation}\label{eq:alpha}
%\alpha^{j} = I_{x}\hat{u}^{j} + I_{y}\hat{v}^{j} 
%\end{equation}
%which can also be written as:
%\begin{equation}
%\alpha^{j} = \hat{k}s^{T}\hat{t}^{j} + v^{T}\hat{w}^{j}  \nonumber
%\end{equation}
%Substituting equations (\ref{eq:alpha}) and (\ref{eq:du.hat}) in equation
%(\ref{eq:modified.CBE}) we get:
%\begin{eqnarray}\label{eq:modified.CBE.2}
%I_{x} u + I_{y} v + (\hat{I}'_{t} - \alpha') =0 \\
%I_{x} u + I_{y} v + (\hat{I}''_{t} - \alpha'') =0 \nonumber
%\end{eqnarray}
%Substituting equation (\ref{eq:u.v.lh}) in equation
%(\ref{eq:modified.CBE.2}) we get modified versions of the equations
%  (\ref{eq:horn.weldon}) and  (\ref{eq:horn.weldon2})
%\begin{eqnarray}
%ks^{T}t' + v^{T}w' + (I'_{t} - \alpha')= 0 \\
%ks^{T}t'' + v^{T}w'' + (I''_{t} - \alpha'')= 0 \nonumber
%\end{eqnarray}
%We start our first iteration with  $\hat{k}$, $\hat{t}^{j}$ , $\hat{w}^{j}$
%all zero and therefore $\alpha=0$ as well. 

%In order to deal with image motions larger than 1 pixel we use
%a Gaussian pyramid for coarse to fine processing
%\cite{Bergen:92}\cite{Burt:83}.





%\subsection{Computing the depth, smoothing and interpolation.}

%After recovering the camera motions we use equations (\ref{eq:horn.weldon}) and
%(\ref{eq:horn.weldon2}) to compute depth at every
%point where $s^{T}t'$ or $s^{T}t''$ are non zero. To combine the
%information from both images and to interpolate over areas where image
%gradients are small we used Local Weighted
%Regression. 
%%This method was chosen because it is simple to implement but
%%could very well be replaced by other methods such as RBF's, B-splines or
%%thin-plate interpolation.

%Equation (\ref{eq:lwr}) shows the cost function used to
%compute the depth at a given point: 
%\small
%\begin{equation}
%\begin{array}{c}
 .%\\
%\min \arg\\
%K
%\end{array}
% \sum_{x,y \in R} \sum_{j} \beta (x,y) |s^{T}t^{j}|^{p} \left(
%ks^{T}t^{j} + v^{T}w^{j} + I^{j}_{t} \right)^{2} \label{eq:lwr}
%\end{equation}
%\normalsize
%The sum is over a region $R$ and over the two motions $j=1,\ldots,2$.
%The windowing function $ \beta (x,y)$ allows one to increase the weight of the
%closer points. We used a function created by convolving two box filters
%together. It is a crude approximation to a Gaussian. The $|S^{T}t^{j}|^{p}$
%term reduces the weight of points which have a small gradient or where the
%gradient is perpendicular to that camera motion since these cases are
%highly affected by noise. We used $p=2$. The size of the region $R$ depends
%on the amount of smoothing and interpolation required. During the iteration
%process we used a region of $5 \times 5$ or $9 \times 9$. In order to get
%'prettier' results, after the last iteration, we reduced the region
%$R$ to $3 \times 3$ but added a multigrid membrane interpolation scheme.



\section{Experiments and Results}
\subsection{Overview}

In this section we perform experiments which test various aspects of
direct estimation of motion and structure.
Experiments~\ref{sec:TBC.exp.PT.Heading},
and~\ref{sec:TBC.optical.flow.compare} deal with the pure translation
case. In experiment~\ref{sec:TBC.exp.PT.Heading}, we test the accuracy
of heading estimation for non-collinear motion.
% Experiment~\ref{sec:TBC.collinear.exp} tests the accuracy heading
% estimates obtained from collinear motion.
Experiment~\ref{sec:TBC.optical.flow.compare} tests a standard optical
flow program on some of the same input images and compares the results
with the `direct methods'.

We then move on to case where the motion involves both translation and
rotation. Experiment~\ref{sec:TBC.exp.cvpr97} revisits the image
sequence used in Stein \& Shashua \cite{Stein-Shashua-cvpr97}. We show
that the improvements described in Section \ref{sec:TBC.lie.to.system}
greatly improve the rotation estimates.  In
Experiment~\ref{sec:Head.With.Rot.exp} we repeat the heading
estimation experiments but this time the motion includes some
rotation. We also show how not taking into account even small
rotations ($\sim 0.5^o$) leads to large errors in heading estimates.
In Experiment~\ref{sec:Euclidean.exp} we test the accuracy of
Euclidean reconstruction using the simple image of a cube.  We measure
whether the right-angles between the cube faces are correctly
recovered.

In all the preceding experiments the camera was mounted on a motion
stage so that ground truth motion could be known accurately. In
Experiment~\ref{sec:Hand.Held.exp} we repeat the Euclidean
reconstruction experiments but this time with a hand-held camera. We
show that if we neglect the rotation the shape recovery is noisy. We
also show that if we neglect radial lens distortion then the
reconstruction is qualitatively good but the Euclidean measurements,
such as angles, are less accurate. This agrees with the observations
regarding lens distortion and the trilinear tensor
in~\cite{Stein-cvpr97}.  Experiments~\ref{sec:church.exp}
and~\ref{sec:BneiBrit.exp} show results with simple outdoor scenes.

%\section{Pure Translation}

\subsection{Accuracy of Heading Estimates for Pure Translation}\label{sec:TBC.exp.PT.Heading}


In this experiment we measure the accuracy of the estimation of
the camera motion direction (heading direction) in the case of pure
translation. We also gauge the effect of the non-linear lens
distortion.

\subsubsection{Experimental Procedure:}

%The camera was mounted on a motion stage with three degrees of
%computer controlled motion (Figure~\ref{fig:TBC.MotionStageFull}):
%horizontal translation, vertical translation and rotation around the
%camera's $Y$ axis. The axes of motion are highlighted in
%Figure~\ref{fig:TBC.MotionStageFull}a.  

%\begin{figure}
%  \centerline{
%    \epsfysize=3.7in \epsffile{\TbcImagesDir stage1mod.ps}}
%  \centerline{ (a)}
%  \centerline{
%    \epsfysize=3.7in \epsffile{\TbcImagesDir stage.2.ps}}
%  \centerline{ (b)}
%  \protect\caption[Two views of the optical bench and the camera mounted
%on the motion stage]{\em Two views of the optical bench and the camera mounted
%on the motion stage. The three axes of motion are highlighted in
%(a). The
%objects on the table are those that appear in Figure~\protect\ref{fig:TBC.collin.motion.input}. }
%  \label{fig:TBC.MotionStageFull}
%\end{figure}

The camera (Pulnix TM9701) was mounted on a motion stage with three degrees
of computer controlled motion: horizontal translation, vertical
translation and rotation around the camera's $Y$ axis
(Fig.~\ref{fig:TBC.Collinear.Schem}).  The camera lens was (4.9mm lens
/ $82^o$ FOV).

The camera $Y$ axis was aligned with vertical axis of the stage.
Initially the camera was positioned so that the optical axis was
aligned with the horizontal axis of the translation stage. The
accuracy of the alignment was $\pm 1^o$.  The camera was rotated
$\alpha^o$ and an image captured. Then the camera was translated
vertically $12.5mm$ and an image captured. Finally the camera was
translated horizontally $10mm$ and a third image captured.  Thus we
captured an image triplet where the camera motion was pure translation
with one of the motions vertical and the second at a heading of
$\alpha^o$. This was repeated for angles $\alpha=0^o,10^o,~...~,
90^o$. 

\begin{figure}
  \centerline{
    \epsfysize=1.9in \epsffile{\TbcImagesDir heading1.ps}}
 \caption[Schematic diagram of the pure translation experiment]
         {\em Schematic diagram of the pure translation
experiment. Camera heading is set by the angle $\alpha$. Motion~2
and Motion~3 are collinear.} 
\label{fig:TBC.Collinear.Schem}
\end{figure}
  


\subsubsection{Results:}

Motion and depth estimation was performed
with and without lens distortion correction. 
Figure~\ref{fig:TBC.MotionEstimate.PT4.9mm} shows the heading
estimates. The heading estimates are within a few degrees. Notice that
when the true heading was inside the FOV the errors were less than
$1^o$. When the true heading was outside the FOV the errors increased
to $1-2^o$. This is a case of a wide angle lens. The errors were
considerably larger when lens distortion was not taken into account
(solid line in Figure~\ref{fig:TBC.MotionEstimate.PT4.9mm}b).

\begin{figure}
  \centerline{
    \epsfxsize=1.6in \epsffile[018 180 594 612]{\PlotsDir head-estim.pt4.9mm.ps}    
    \hfill
    \epsfxsize=1.6in \epsffile[018 180 594 612]{\PlotsDir head-estim-err.pt4.9mm.ps}}
  \centerline{ \hfill(a)\hfill \hfill (b) \hfill}
    \centerline{\begin{tabular}{|ll|}\hline
       Legend: & \\
       -----  & No radial distortion correction\\
       - - -  & Images corrected for radial distortion\\\hline
    \end{tabular}}
  \caption[Heading estimate. Pure
   translation was assumed. $4.9_{mm}$ lens]{\em (a) Heading estimate from a translating camera. Pure
   translation was assumed. The $4.9_{mm}$ lens gives a $82^o$ FOV along the X
   axis. (b) Difference between heading estimate and motion stage
   reading (``True Heading''). }
  \label{fig:TBC.MotionEstimate.PT4.9mm}
\end{figure}

\subsection{Comparison with Optical Flow Techniques}\label{sec:TBC.optical.flow.compare}

After computing the motion and depth from an image triplet we can
compute the motion flow using the motion equations
(equation~\ref{eq:TBC.u.v.lh}). In this experiment we compare the optical
flow estimates obtained in this manner with the flow estimates
computed directly from two of the images using an `industry standard'
optical flow program. The program is based on code by Bergen \&
Hingorani of the Sarnoff Corp \cite{Bergen-Hingorani90,Anandan-eccv92}. We use two
sets of images which have been chosen to be particularly difficult for
optical flow programs. 

\subsubsection{Experimental Procedure}

Figure~\ref{fig:lab213.bump.ps} shows three images of a real scene:
two paper cylinders, a plaster bust and a black metal bar in front of
a background of vertical stripes.  The two motions are pure
translation, vertical and horizontal, parallel to the image plane.
\mbox{Figure~\ref{fig:lab213.bump.ps}a} shows the recovered depth map for the
sequence. The results are qualitatively correct. In
Figure~\ref{fig:lab213.bump.ps}b the texture mapping was removed for
clarity. The difficulties in the scene are explained in the
introduction (Section~\ref{sec:introduction}) and expanded upon in
the results section below. 

From the depth map and motion the motion flow was computed using
equations (equation~\ref{eq:TBC.u.v.lh}). The flow field was also computed
using the optical flow program and the results compared.

This experiment was repeated for a second set of images shown in
Figure~\ref{fig:lab263.input}a,b,c. The main difficulty here is
that most of the gradients are smooth a low contrast. The strongest
edges appear at the occluding contour of the bust and these edges do
not represent real features.


\subsubsection{Results}

Figure~\ref{fig:TBC.LK.compare.212}a shows the estimated flow
computed from the recovered depth and motion. The camera motion was
horizontal parallel to the image plane. The flow vectors have a zero
$Y$ component. Figure~\ref{fig:TBC.LK.compare.212}b shows the
magnitude of the flow vectors. \mbox{Figure~\ref{fig:TBC.LK.compare.212}c}
shows the flow estimated from two images using an optical flow
program. Near the vertical bars on the left the flow vectors have a
strong $Y$ component. This is due to the aperture problem. Even a
large window (aperture) will see image gradients in only one direction
giving no constraint on the $Y$ component of the flow. 

On the cylinder to the lower right of the bust the magnitude of the
flow vectors is too small. This is due to the fact that the edges were
horizontal and parallel to the horizontal epipolar lines. This means
that the depth estimates for theses points are unreliable
(Figure~\ref{fig:TBC.LK.compare.212}d). This problem does not occur in
Figure~\ref{fig:TBC.LK.compare.212}b because two motions one
horizontal and one vertical were used to compute the original depth map.

An enlarged detail of one of the more difficult regions is shown
Figure~\ref{fig:TBC.LK.compare.212detail}. A critical error occurs
near the intersections (in the image) of the diagonal line and the
vertical bars.

%The flow has a small $Y$ component. This error occurs because the
%program 'tracks' the intersection point as if it were a real 'feature'
%point, but the lines do not intersect in space. Although small in
%magnitude this is the more significant error because the program will
%also give a high confidence to this value.



\begin{figure}
          \centerline{ 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/flow.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/flowMag.ps}
               }
          \centerline{   \hfill (a) \hfill \hfill (b) \hfill}
          \centerline{ 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/lkFlow.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/lkFlowMag.ps}
               }
          \centerline{   \hfill (c) \hfill \hfill (d) \hfill}
\caption[Optical flow computed for
Figure~\protect\ref{fig:lab213.input.ps}]
{\em Optical flow computed for
Figure~\protect\ref{fig:lab213.input.ps}. (a) Correct optical flow
computed using depth map and recovered motion. The camera motion was
horizontal parallel to the image plane. The flow vectors have a zero
$Y$ component. (c) Optical flow computed using code from Bergen \&
Hingorani. Near the vertical bars on the left the flow vectors have a
strong $Y$ component. This is due to the aperture problem. Even a
large window (aperture) will see image gradients in only one direction
giving no constraint on the $Y$ component of the flow.  }
\label{fig:TBC.LK.compare.212}
\end{figure}


\begin{figure}
          \centerline{ 
               %\hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/flowDetail.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/lkFlowDetail.ps}}
          \centerline{   \hfill (a) \hfill \hfill (b) \hfill}
\caption[Detail of optical flow computed for
Figure~\protect\ref{fig:lab213.input.ps}]{\em (a) Detail from the top
left corner of Figure~\protect\ref{fig:TBC.LK.compare.212}a. (b) Detail from the top
left corner of Figure~\protect\ref{fig:TBC.LK.compare.212}c. Along the 
vertical bars there is a strong Y component but this is to be expected 
due to the aperture problem. Lower down, near the
intersections (in the image) of the diagonal line and the vertical
bars the flow has a small $Y$ component. This error occurs because the
program 'tracks' the intersection point as if it were a real 'feature'
point, but the lines do not intersect in space. Although smaller in
magnitude this is the more significant error because the program will
also give a high confidence to this value. }
\label{fig:TBC.LK.compare.212detail}
\end{figure}

The results for the second set of images are shown in
Figure~\ref{fig:TBC.LK.compare.262}. Errors in flow direction can be
seen near the edges of the bust. The most significant errors are
errors in magnitude of the flow in the background areas. These can
been seen clearly by comparing the mesh plots.

\begin{figure}
          \centerline{ 
               %\hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/flow263.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/flowMag263.ps}
              \hfill }
          \centerline{ 
              \hfill (a) $\uparrow$ \hfill(c) $\downarrow$ \hfill 
              \hfill (b) $\uparrow$ \hfill(d) $\downarrow$ \hfill 
              } 
          \centerline{ 
               %\hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/lkFlow263.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{lab4-images/flow/lkFlowMag263.ps}
              \hfill }
  \caption{\em Optical flow computed for
Fig.~\protect\ref{fig:lab263.input}. (a) Correct optical flow
computed using depth map and recovered motion. The camera motion was
horizontal parallel to the image plane. The flow vectors have a zero
$Y$ component. (c) Optical flow computed using code from Bergen \&
Hingorani. Errors in flow direction can be seen near the edges of the
bust. The most significant errors are errors in magnitude of the flow
in the background areas. These can been seen clearly by comparing the
mesh plots. }
\label{fig:TBC.LK.compare.262}
\end{figure}

This experiment was intended to highlight the problems with optical
flow techniques and how they are overcome by using direct
methods. While this was not a simple 'straw man' and the optical flow
program was one of the best available, a full system would based on
optical flow would have more components intended to overcome the
problems described. Together with flow vectors the optical program
would compute a confidence ellipse for each vector. These estimates
would then be combined using robust methods to reject outliers during
the motion computation. After the camera motion was estimated the flow
would be recomputed using the epipolar constraint and the new flow
would be used to compute the depth. Probably three views would be
used with different camera motions. Building such a worthy competitor
and fine tuning it to these image sets is of course beyond the scope
of this paper but this description does give an idea of the
complexity of such a system. 

%\section{Translation and Rotation}
\subsection{Translation and Rotation
(Reanalysis of Data from Stein \&
Shashua, CVPR97)}\label{sec:TBC.exp.cvpr97}

In this experiment we test the modified algorithm for recovering both
translation and rotation (Section~\ref{sec:TBC.lie.to.system}) on the
same data set used by Stein \& Shashua \cite{Stein-Shashua-cvpr97}.

\subsubsection{Experimental Procedure:}

The images used by \cite{Stein-Shashua-cvpr97} were taken with an
$8.5mm$ lens in the following way: The camera was translated first
vertically ($10mm$) and then horizontally ($5mm$) to the right.  At
this third position the camera was rotated to various angles ranging
from $-4.0^o$ to $1.0^o$. The depth in the scene ranged from $170mm$
to $400mm$. 

Figs.~\ref{fig:lab263.input}a,b,c show the three of the input
images in the sequence.  The flow due to a rotation of $-1.8^o$ in
Fig.~\ref{fig:lab263.input}b is much greater than the flow due to
the translation. The direction of
the flow due to rotation by a negative angle was in the same direction
as the image flow induced by the translation.

\begin{figure}
  \centerline{ 
    \epsfxsize=1.6in  \epsffile{lab4-images/lab262.s.ps}
    \hfill 
    \epsfxsize=1.6in \epsffile{lab4-images/lab270.s.ps} }
  \centerline{ \hfill (a) \hfill \hfill (b) \hfill  }
  \centerline{ 
    \epsfxsize=1.6in \epsffile{lab4-images/lab260.s.ps} 
    \hfill 
    \epsfxsize=1.6in \epsffile{lab4-images/lab262260270s.depth.0.ps}}
  \centerline{ \hfill (c) \hfill \hfill (d) \hfill  }
  \centerline{ 
    \epsfxsize=1.6in \epsffile{lab4-images/lab262260270s.smooth.ps} 
    \hfill 
    \epsfxsize=1.6in \epsffile{lab4-images/lab262260270s.notex2.ps}} 
  \centerline{ \hfill (e) \hfill \hfill (f) \hfill  }
\caption{\em The three input images (a,b,d) used for 3D
  reconstruction. Motion $(a \rightarrow b)$ was a horizontal
  translation with a $-1.8^o$ rotation around the $Y$ axis. Motion $(a
  \rightarrow c)$ was a vertical translation. This is a challenging
  set of images because the the texture is smooth with low contrast.
  The strongest edge 'features' appear along the occluding contour of
  the head but this edge does not in fact correspond to a real feature
  in the world. (d) The estimated depth map. The estimated rotation
  was $-1.62^o$. There are of course errors around the boundary of the
  image where there is no overlap between the images. (e), (f) 3D
  rendering of depth map in (d)}
\label{fig:lab263.input}
\end{figure}

\subsubsection{Results:}

Fig.~\ref{fig:lab263.input}d shows the recovered depth map using
images in Figs.~\ref{fig:lab263.input}a,b,c.
\mbox{Figs.~\ref{fig:lab263.input}e,f} show the 3D rendering of the
depth map. 

Fig.~\ref{fig:est-rotation.ps} shows the recovered rotation
estimates for true rotations ranging from $-4.0^o$ to
$0.8^o$. Outside this range the old version of the algorithm did not
converge. Using the new algorithm the rotation estimates are within
$5\%$ of the correct value. There is no scaling error as was observed
using the old method. 

The new algorithm also has a wider range of convergence. The original
images were not available for larger rotations but a similar setup
yielded good results for rotation for the range of $-5.0^o$ to
$5.0^o$. Since rotations can create large amounts of image flow, key
to the successful convergence are strong signals with low spatial
frequency such as the white head on a darker background. This leads to
convergence at the coarse level which is then used as the starting
point for motion estimates at finer levels. Strong repetitive high
frequency patterns such as checker board patterns are more
difficult. 

\begin{figure}
  \centerline{
    \epsfxsize=1.6in \epsffile[018 180 594 612]{lab4-images/est-rotation.ps}
    \hfill
    \epsfxsize=1.6in \epsffile[018 180 594 612]{lab4-images/est-rotation-new.ps}}
  \centerline{ \hfill (a) \hfill \hfill(b) \hfill }
%  \centerline{
%    \epsfysize=3.0in \epsffile[018 180 594 612]{lab4-images/est-rotation.ps}}
%  \centerline{ (a)}
%  \centerline{
%    \epsfysize=3.0in \epsffile[018 180 594 612]{lab4-images/est-rotation-new.ps}}
%  \centerline{ (b)}
\caption[Estimated rotation as a function of true rotation]{\em 
  Estimated rotation as a function of true rotation. (a)~Old results:
  Although there is the correct linear relationship there also appears
  to be a significant scale error (note Y axis). (b)~New results - the
  rotation estimates are correct. }
  \label{fig:est-rotation.ps}
\end{figure}

\subsection{Motion Heading Estimation with Rotation}\label{sec:Head.With.Rot.exp}

The aim of this experiment is to test the motion estimation algorithm
over a wide variety of motion directions when the motion included some
rotation. For example, this experiment 
tests whether the algorithm works when the
direction of motion is inside and outside the FOV. We also find the
error created by neglecting even small amounts of rotation. The two
translations are non-collinear. 

\subsubsection{Experimental Procedure:}

The camera setup is similar to that of experiment~1
(Section~\ref{sec:TBC.exp.PT.Heading}). The camera is rotated to a
particular heading direction $(0^o, 10^o, ... , 90^o)$. The camera was
translated vertically and then horizontally. After the horizontal
translation the camera was rotated $\pm 0.5^o $ and $\pm 1.0^o
$. Therefore, for each image heading, we have an image triplet which
is translation only and then triplets where the second motion also
includes some small rotations.  The rotation axis of the stage passed
within $5mm$ of the camera center of projection thus rotations of
$1.0^o$ produce translations of less that $1mm$.
Figure~\ref{fig:TBC.Rotation1Input} shows input images where the
horizontal translation heading was $60^o$ to the camera optical axis.

\begin{figure}
  \centerline{ 
    \epsfxsize=1.6in \epsffile{\TbcImagesDir seq60ss.005.ps}
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcImagesDir seq60ss.006.ps}}
  \centerline{ \hfill (a) \hfill \hfill (b) \hfill }
  \centerline{    
    \epsfxsize=1.6in \epsffile{\TbcImagesDir seq60ss.014.ps} 
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcImagesDir seq60ss.003.ps}}
  \centerline{ \hfill (c) \hfill \hfill (d) \hfill }
  \centerline{    
    \epsfxsize=1.6in \epsffile{\TbcImagesDir seq60ss.051403.depth.0.ps}
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcImagesDir seq60ss.051303pt.depth.0.ps}}
  \centerline{ \hfill (e) \hfill \hfill (f) \hfill }
  \caption[Experiment to test motion
  heading estimation in the presence of rotation]{\em Four of the
    input images used in the experiment to test motion heading
    estimation in the presence of rotation. Motion (a) to (d) is
    vertical motion.  Motion (a) to (b) is horizontal motion forward
    and to the left ($60^o$ from the camera optical axis).  Motion (a)
    to (c) has the same translation as (a) to (b) but also with a
    rotation of $1.0^o$. In this case the flow due to rotation of
    $1.0$ cancels out the flow due to translation in the area of the
    face and more than cancels out the flow in regions of greater
    depth. (e) Depth map estimated after computing both translation
    and rotation. The second camera motion included $-1.0^o$ rotation.
    (f)~Depth map estimated assuming pure translation when the second
    camera motion included $0.5^o$ rotation. Notice the depth errors
    in the background.(Lens $8.5mm$ lens / $52^o$
    FOV)}\label{fig:TBC.Rotation1Input}
\end{figure}

%\begin{figure}
%          \centerline{ 
%               \hfill 
%                \epsfysize=2.2in \epsffile{\TbcImagesDir seq60ss.051403.depth.0.ps}
%              \hfill 
%                \epsfysize=2.2in \epsffile{\TbcImagesDir seq60ss.051303pt.depth.0.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (a) \hfill \hfill (b) \hfill 
%              } 
%\caption[Depth map estimated after computing both translation and
%rotation]{\em (a) Depth map estimated after computing both translation
%and rotation. The second camera motion included $-1.0^o$
%rotation. (b)~Depth map estimated assuming pure translation when the
%second camera motion included $0.5^o$ rotation. Notice the depth
%errors in the background. }
%\label{fig:TBC.RotationDepthMap}
%\end{figure}

\subsubsection{Results:}

Depth and motion estimates were computed either assuming pure
translation or allowing possible rotation. The results for various
motion and motion-assumption combinations are shown in
Figure~\ref{fig:TBC.MotionEstimate.8.5mm},
Figure~\ref{fig:TBC.MotionEstimateErr.8.5mm} and
Figure~\ref{fig:TBC.RotationEstimate.8.5mm}. We notice that for even a
small rotation ($0.5^o$) the heading estimate is off by over $10^o$ if
we do not take the rotation into account. When we estimate both
rotation and translation the motion estimates are good (RMS error
$2.1^o$) although not as good as those obtained from a purely
translating camera under the pure translation assumption (RMS error
$1.44^o$).

Figure~\ref{fig:TBC.Rotation1Input}e shows the depth map estimated
from the images in Figure~\ref{fig:TBC.Rotation1Input}a),
Figure~\ref{fig:TBC.Rotation1Input}c and Figure~\ref{fig:TBC.Rotation1Input}d
allowing for both translation and rotation. The depth map is
qualitatively correct even though the second camera motion included
$-1.0^o$ rotation. Figure~\ref{fig:TBC.Rotation1Input}f shows the
depth map estimated assuming pure translation when in fact the second camera
motion included $0.5^o$ rotation. 


\begin{figure}
  \centerline{
    \epsfxsize=3.2in \epsffile[018 180 594 612]{\PlotsDir head-estim.8.5mm.ps}}
%    \centerline{\begin{tabular}{|lll|}\hline
%       Legend: & &\\
%       ----+-  & Motion: pure translation & Assumed: pure translation\\
%       ----*-  & Motion: $0.5^o$ rotation & Assumed: pure translation\\
%       - - -  &  Motion: $+1.0^o$ rotation & Assumed: possible rotation\\
%       \_.\_.\_  & Motion: pure translation & Assumed: possible rotation\\
%       .....*. &  Motion: $-1.0^o$ rotation & Assumed: possible rotation\\\hline
%    \end{tabular}}
\caption[Heading estimates from a moving camera]{\em Heading estimates
from a moving camera.  (PT - PT) Both true motion motion model was
pure translation. ($0.5^o$ - PT) True motion included $0.5^o$ and
motion model was pure translation. Note that even a small rotation, such as $0.5^o$, is enough to
create large heading estimation errors if rotation is not taken into
account. (PT- Rot, $+1.0^o$ - Rot, $-1.0^o$- Rot) Motion model
included rotation Actual rotations were $0^o, 1.0^o,
-1.0^o$.  ($8.5_{mm}$ lens / $52^o$ FOV) }
\label{fig:TBC.MotionEstimate.8.5mm}
\end{figure}

\begin{figure}
  \centerline{
    \epsfxsize=3.2in \epsffile[018 180 594 612]{\PlotsDir head-estim-err.8.5mm.ps}}
%    \centerline{\begin{tabular}{|lll|}\hline
%       Legend: & &\\
%       ----+-  & Motion: pure translation & Assumed: pure translation\\
%       - - -  &  Motion: $-1.0^o$ rotation & Assumed: possible rotation\\
%       \_.\_.\_  & Motion: pure translation & Assumed: possible rotation\\
%       .....*. &  Motion: $-1.0^o$ rotation & Assumed: possible rotation\\\hline
%    \end{tabular}}

\caption[Heading estimate error: the difference between heading
estimate and motion stage reading]{\em Heading estimate error: the
difference between heading estimate and motion stage reading. Even in
the presence of camera rotation ($+1.0^o$ - Rot, $-1.0^o$- Rot) the
heading estimation errors are small (RMS error $2.1^o$) although not
as good as obtained from purely translating camera under the pure
translation assumption (PT-PT) where the RMS error was $1.44^o$. 
%Note that
%the estimates obtained from the motions that included rotation appear
%offset by about $1^o$ to $2^o$. This offset is a correct result since
%the heading estimates are given in the second camera
%coordinates. 
($8.5mm$ lens / $52^o$ FOV)}
\label{fig:TBC.MotionEstimateErr.8.5mm}
\end{figure}

\begin{figure}
  \centerline{
    \epsfxsize=3.2in \epsffile[018 180 594 612]{\PlotsDir rotation-estim.8.5mm.ps}}
\caption[Rotation estimates from a camera translating  with some added
rotation] {\em Rotation estimates obtained from a camera translating
in different translation directions with some added rotation. Actual
rotations were $0^o, 1.0^o, -1.0^o$. (camera with $8.5_{mm}$ lens /
$52^o$ FOV)}
\label{fig:TBC.RotationEstimate.8.5mm}
\end{figure}

%\begin{figure}
%          \centerline{ 
%               \hfill 
%                \epsfysize=2.2in \epsffile{\TbcImagesDir seq60ss.051403.depth.0.ps}
%              \hfill 
%                \epsfysize=2.2in \epsffile{\TbcImagesDir seq60ss.051303pt.depth.0.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (a) \hfill \hfill (b) \hfill 
%              } 
%\caption[Depth map estimated after computing both translation and
%rotation]{\em (a) Depth map estimated after computing both translation
%and rotation. The second camera motion included $-1.0^o$
%rotation. (b)~Depth map estimated assuming pure translation when the
%second camera motion included $0.5^o$ rotation. Notice the depth
%errors in the background. }
%\label{fig:TBC.RotationDepthMap}
%\end{figure}

\subsection{Euclidean Structure Estimation  in the Presence of Rotation.  }\label{sec:Euclidean.exp}

In this section we evaluate the Euclidean structure estimation when
the motion includes rotation. We use a simple scene, a cube, in front
of a planar background. We perform the 3D reconstruction of the
scene. We then measure whether the right-angles between the cube faces
are correctly recovered.

\subsubsection{Experimental Procedure}

Figure~\ref{fig:TBC.CubeInput} shows the simple scene of a textured
cube in front of flat, low-contrast, background. The camera ($8.5mm$
lens / $52^o$ FOV) was translated vertically and then horizontally
parallel to the image plane. The camera was then rotated to $+2.0^o$
and $-1.0^o$ from its original heading.

Since we were interested in the 3D structure and not jut motion
estimates, after estimating both structure and motion we went back
down the pyramid and recomputed the structure while keeping the motion
estimates constant (Section~\ref{sec:TBC.final.touches}). The depth
estimates before and after the second pass are shown in
Figure~\ref{fig:TBC.CubeInput}d and Figure~\ref{fig:TBC.CubeInput}f
respectively. In Figure~\ref{fig:TBC.CubeInput}d there is a ``hole''
on the right of the depth map. This is due to the appearance of a
white patch near the border of the image. As with occlusions this
sudden appearance violates the constant brightness assumption. The
situation is fixed by the second pass. In general points near the edge
of the depth map are unreliable.

\begin{figure}
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir seq90iss.005.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir seq90iss.003.ps} 
              \hfill }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir  seq90iss.008.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir seq90iss.050803.W3p0.depth.0.ps}
              \hfill }
                    \centerline{ 
              \hfill (c) \hfill \hfill (d) \hfill 
              } 
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir  seq90iss.015.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir seq90iss.050803.W3.depth.0.ps}
              \hfill }
                    \centerline{ 
              \hfill (e) \hfill \hfill (f) \hfill 
              } 
\caption[A simple scene used to test 3D Euclidean
reconstruction: a cube]{\em A simple scene used to test 3D Euclidean
reconstruction: a cube constructed from cork blocks in front of a flat
background. Motion (a) to (b) was vertical. Motion (a) to (c) was
sideways translation with $-1.0^o$ rotation. Motion (a) to (e) was
same sideways translation but with $+2.0^o$ rotation. (d) initial
depth map from images (a), (b) and (c). Note how the disappearance of
a white blob in the background near the right edge of the image
creates an error in one region of the depth map. (f) After second pass
of depth estimation the ``hole'' is fixed. (see text)}
\label{fig:TBC.CubeInput}
\end{figure}

\subsubsection{Results:}

The rotation estimates were $1.68^o$ and $-0.67^o$ for true rotations
of $+2.0^o$ and $-1.0^o$ respectively. The 3D rendering of the cube is
shown in Figures~\ref{fig:TBC.CubeRes08}
and~\ref{fig:TBC.CubeRes15}. The wire-frame rendering of the overhead
views (Figures~\ref{fig:TBC.CubeRes08}b and~\ref{fig:TBC.CubeRes15}b) show the recovered angle between the
cube faces. The estimated angles were $95^o$ and $86^o$
respectively. Figure~\ref{fig:TBC.CubeRes08}d shows a side view of the
cube. The estimated angle between the top and side face of the cube is
$90^o$.

\begin{figure}
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir overhead1.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir wire1.ps}
              \hfill }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir wireside1.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir side2.ps}
              \hfill }
                    \centerline{ 
              \hfill (c) \hfill \hfill (d) \hfill 
              } 
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir side1.ps}
              \hfill }
                    \centerline{ 
              \hfill (e) \hfill 
              } 
\caption[3D rendering of the recovered scene]{\em 
3D rendering of the recovered scene. The camera motion included
$-1.0^o$ rotation. (a) Overhead view. The estimated camera position is
marked by the sphere in bottom right-hand corner. (b) An enlarged view
of the corner of the cube. The angle between the faces of the cube is
estimated $95^o$ instead of $90^o$. (c) Side view of the segmented
cube from up and to the side. (d) enlarged wire-frame side view of the
cube. The estimated angle between the top face of the cube and the
right face is $90^o$. Another view of the scene showing the
foreground/background segmentation. The depth estimates of the back
panel are noisy because they are further away from the camera and
because the texture does not have strong gradients. }
\label{fig:TBC.CubeRes08}
\end{figure}

\begin{figure}
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir side015.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcImagesDir overhead015.ps}
              \hfill }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
\caption[3D rendering of the recovered cube]{\em 3D rendering of the
recovered cube. The camera motion included $+2.0^o$ rotation.  The
cube has been segmented from the background using depth based
foreground/background segmentation. (a) View from up and to the side.
(b) Wire-frame rendering of the overhead view. From the overhead view
the estimated angle between the cube faces is $86^o$ instead of
$90^o$} 
\label{fig:TBC.CubeRes15}
\end{figure}

%\section{Experiments with a Hand-held Camera}
\subsection{Euclidean Structure
Estimation Using a Hand-held Camera}\label{sec:Hand.Held.exp}

In previous experiments the camera was mounted on a motion stage so
that ground truth motion values were available. In this experiment the
camera was hand-held to show that the system can deal with the more
natural situation. Simple objects were used to create the scene so that
it would be easier to evaluate the shape reconstruction. 

\subsubsection{Experimental Procedure}

A progressive scan camera (Pulnix TMC9701) with a wide angle lens
($4.9mm$ lens / $82^o$ FOV) was hooked up directly to the frame grabber
of an SGI Indy workstation. The cameras help in the hand with no
tripod. An image sequence of 30 frames was captured while the camera
was moved up and down and side to side the combination of which
produced a circular motion.  An effort was made to keep the camera
rotation down to a minimum.

From the image sequence three images were selected which seemed to
give two distinct motion directions.
Figures~\ref{fig:CubeHandheldInput}a,b,c show the three input
images of the cube.

\begin{figure}
          \centerline{ 
               \hfill 
                \epsfxsize=1.0in \epsffile{\TbcHandheldImagesDir cube2.015.ps}
              \hfill 
                \epsfxsize=1.0in \epsffile{\TbcHandheldImagesDir cube2.008.ps}
              \hfill 
                \epsfxsize=1.0in \epsffile{\TbcHandheldImagesDir  cube2.021.ps}
              \hfill }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill  \hfill (c) \hfill 
              } 
          \centerline{ 
               \hfill 
                \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir side1K.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir side1wireK.ps}
              \hfill}
                    \centerline{ 
              \hfill (d) \hfill \hfill (e) \hfill 
              } 
\caption{\em A simple scene used to test 3D Euclidean
  reconstruction. Camera was hand-held so exact motion between images
  (a), (b) and (c) is not known. (d) Depth reconstruction allowing for
  both rotation and translation with radial distortion correction. (e)
  Assumes pure translation. (f) Allows for both rotation and
  translation but no radial distortion correction. The cube and it's
  surrounding are recovered well while textureless areas near the
  boundaries have large errors. (d),(e) 3D renderings of the recovered
  shape. The motion estimation allowed for rotation and translation
  and the input images were preprocessed to correct for lens
  distortion. The estimated angle between cube faces is $96^o$ instead
  of $90^o$. }
\label{fig:CubeHandheldInput}
\end{figure}


\subsubsection{Results}

The motion and depth was estimated for three different cases. First
the motion and depth was estimated allowing for both rotation and
translation and taking into account the radial distortion. In the
second case motion was estimated allowing for both rotation and
translation but neglecting radial distortion. The resulting depth map
is qualitatively the same but the measured angles between the cubes
faces are larger and further from the true value of $90^o$
(Fig.~\ref{fig:CubeHandheldTop}c,d). In the third case assumed pure
translation the resulting depth map again appears to be qualitatively
the same but the 3D renderings (Figure~\ref{fig:CubeHandheldTop}e,f)
clearly show that the results are more noisy and further from true
Euclidean reconstruction.

Errors in depth reconstruction when we neglect the rotation are to be
expected according to \cite{Dutta-Snyder90}. They show that for small
total rotation angles ($|w|$) the error in depth due to an error in
rotation ($|\delta w|$) is given by:
$$
\frac{\delta Z}{Z} \approx \left( \frac{f}{| u |}\right)  \left(|
\delta w |\right)
$$ 
where $f$ is the focal length, $|u|$ is the magnitude of the image
motion. In our case by assuming pure translation we get a rotation
error of $1.45^o$ or $0.025 rad$
(Table~\ref{table:TBC.CubeHandheldRot}). The average image motions
were on the order of $15 pixels$ and the focal length $f \approx
180$. So we expect:
$$
\frac{\delta Z}{Z} \approx 0.3. 
$$
\cite{Dutta-Snyder90} refers to depth from image motion. 
The errors we observe are not so large but are definitely noticeable. 
We also note that even if we neglect the rotation none of the depth
estimates were negative so a {\em depth is positive} constraint, as
suggested by \cite{Fermuller-Aloimonos-97}, would not help in this
example. 


\begin{figure}
  \centerline{ 
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir top2K.ps}
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir top1wireK.ps}
    \hfill }
  \centerline{ 
    \hfill (a) \hfill \hfill (b) \hfill 
    } 
  \centerline{ 
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir top1noK.ps}
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir top1wirenoK.ps}
    \hfill }
  \centerline{ 
    \hfill (c) \hfill \hfill (d) \hfill 
    } 
  \centerline{ 
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir side1PT.ps}
    \hfill 
    \epsfxsize=1.6in \epsffile{\TbcHandheldImagesDir tope1wirePT.ps}
    \hfill }
  \centerline{ 
    \hfill (e) \hfill \hfill (f) \hfill 
    } 
  
  \protect\caption{\em Top view of the cube scene in
    Figure~\protect\ref{fig:CubeHandheldInput}. (a),(b) Motion and
    depth was estimated allowing for rotation. The images were
    preprocessed for radial distortion correction.  The estimated
    angle between cube faces was $97^o$. (c),(d) No radial distortion
    correction. The estimated angle between cube faces was $101^o$.
    (e),(f) Motion was estimated assuming pure translation.  There are
    many errors in the depth map and the angles estimate is $105^o$. }
\label{fig:CubeHandheldTop}
\end{figure}

\begin{table}\centering 
\protect\caption[Rotation estimates (in degrees) from the images in
Figure~\protect\ref{fig:CubeHandheldInput}]{\em Rotation estimates (in degrees) from the images in
Figure~\protect\ref{fig:CubeHandheldInput}. The camera was
hand-held so ground truth values of the motion are not known.}\label{table:TBC.CubeHandheldRot}
\begin{tabular}{||c|c|c|c||}\hline \hline
         &  $w_x$ &  $w_y$ &  $w_z$ \\
Motion 1 &   -0.7809  &   0.2262 &   -0.5520\\
Motion 2 &  -0.0538  &   0.2684  &   1.4246\\
\hline \hline
\end{tabular}
\end{table}


%\begin{figure}
%          \centerline{ 
%               \hfill 
%                \epsfxsize=3.0in \epsffile{\TbcHandheldImagesDir bust2.020.ps}
%              \hfill 
%                \epsfxsize=3.0in \epsffile{\TbcHandheldImagesDir bust2.026.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (a) \hfill \hfill (b) \hfill 
%              } 
%          \centerline{ 
%               \hfill 
%                \epsfxsize=3.0in \epsffile{\TbcHandheldImagesDir  bust2.014.ps}
%              \hfill 
%                \epsfxsize=3.0in \epsffile{\TbcHandheldImagesDir bust2.202614.depth.0.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (c) \hfill \hfill (d) \hfill 
%              } 
%\caption[Three input views from a hand-held camera and the
%resulting depth map]{\em Three input views from a hand-held camera and the
%resulting depth map. The motion model included both  rotation and
%translation. } 
%\label{fig:TBC.BustHandheldInput}
%\end{figure}

%\begin{figure}
%          \centerline{ 
%               \hfill 
%                \epsfxsize=2.4in \epsffile{\TbcHandheldImagesDir bust2side1.ps}
%              \hfill 
%                \epsfxsize=2.4in \epsffile{\TbcHandheldImagesDir bust2sidewire1.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (a) \hfill \hfill (b) \hfill 
%              } 
%          \centerline{ 
%               \hfill 
%                \epsfxsize=2.4in \epsffile{\TbcHandheldImagesDir bust2top1.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (c) \hfill
%              } 
%\protect\caption{\em 3D renderings of the depth map in
%Figure~\protect\ref{fig:TBC.BustHandheldInput}d. } 
%\label{fig:TBC.BustHandheld}
%\end{figure}

\subsection{Outdoor Scenes - Church Wall}\label{sec:church.exp}

In this experiment we tested the algorithm on outdoor scenes. These
experiments are important because they show that the photometric
constraints can be used in uncontrolled lighting conditions. 
Figures~\ref{fig:TBC.JamesInput}a,b,c show three images of a small
part of a church wall. The camera was Sont Hi8 camcorder with the lens
open as wide as possible (FOV $40^o$). The camera was mounted on a
tripod. Moving the camera horizontally ($a \rightarrow b$) was quite
smooth but moving the camera up and down included about $0.5^o$
rotation. 
Figure~\ref{fig:TBC.JamesInput}d shows the recovered depth map which
can be seen to be qualitatively correct. The Euclidean 3D renderings
of the depth map are shown in
Figure~\ref{fig:TBC.JamesRes}. 

%Figure~\ref{fig:TBC.JamesRes}c shows a
%view from the far left. This part of the wall is closest to the camera
%and so has much less noise than the section of wall far from the
%camera on the right of figure Figure~\ref{fig:TBC.JamesRes}d.  Note
%the correct right-angle between the ground and the wall in the right
%of figure Figure~\ref{fig:TBC.JamesRes}d.

\begin{figure}
          \centerline{ 
                \epsfxsize=1.6in \epsffile{\TbcJamesImagesDir james1.010.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcJamesImagesDir james1.000.ps}}
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
          \centerline{ 
                \epsfxsize=1.6in \epsffile{\TbcJamesImagesDir james1.040.ps}
              \hfill
                \epsfxsize=1.6in \epsffile{\TbcJamesImagesDir james1.104000.depth.0.ps}}
                    \centerline{ 
              \hfill (c) \hfill \hfill (d) \hfill 
              } 
\caption[Detail of church wall]{\em Detail of church wall: Three input
views and resulting depth map. (Sony Hi8 camcorder, FOV $40^o$) } 
\label{fig:TBC.JamesInput}
\end{figure}


\begin{figure}
          \centerline{ 
                \epsfxsize=3.2in \epsffile{\TbcJamesImagesDir james1.1.ps}}
          \centerline{ 
              \hfill  (a) \hfill}
          \centerline{ 
                \epsfxsize=3.2in \epsffile{\TbcJamesImagesDir james1.2.ps} }
          \centerline{ 
              \hfill  (b) \hfill}
%         \centerline{ 
%               \epsfxsize=2.9in \epsffile{\TbcJamesImagesDir james1.3.ps} 
%             \hfill 
%               \epsfxsize=2.9in \epsffile{\TbcJamesImagesDir james1.4.ps} }
%                   \centerline{ 
%             \hfill (c) \hfill \hfill (d) \hfill 
%             } 
\protect\caption{\em 3D Euclidean renderings of the depth map in
Figure~\protect\ref{fig:TBC.JamesInput}d. } 
\label{fig:TBC.JamesRes}
\end{figure}

\subsection{Outdoor Scene 2 -
Side Entrance}\label{sec:BneiBrit.exp}

Figure~\ref{fig:TBC.BneiInput} shows three views of the side entrance
to a building. The camera was a progressive scan camera (Pulnix
TMC9701) with a wide angles lens (4.9mm / FOV $82^o$). The camera
output was recoded on a Sony Hi8 camcorder. The camera was mounted on
a tripod and an effort was made to reduce the amount of rotation, but
rotation was $0.32^o$.  Figure~\ref{fig:TBC.BneiInput}d shows the
recovered depth map. The 3d Euclidean rendering of the depth map is
shown in Figure~\ref{fig:TBC.BneiRes}.

\begin{figure}
          \centerline{ 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.030.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.010.ps}}
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
          \centerline{ 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.050.ps}
              \hfill
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.305010.depth.0.ps}}
                    \centerline{ 
              \hfill (c) \hfill \hfill (d) \hfill 
              } 
\caption[Three input views of a side entrance to a building. ]{\em
Three input views of a side entrance to a building and resulting depth
map. Camera was mounted on a tripod but motion included up to $0.5^o$
rotation. }
\label{fig:TBC.BneiInput}
\end{figure}

\begin{figure}
          \centerline{ 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.1.ps}
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.2.ps} }
                    \centerline{ 
              \hfill (a) \hfill \hfill (b) \hfill 
              } 
          \centerline{ 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.3.ps} 
              \hfill 
                \epsfxsize=1.6in \epsffile{\TbcBneiImagesDir bnei1.wire1.ps}}
                    \centerline{ 
              \hfill (c) \hfill \hfill (d) \hfill 
              } 
\protect\caption[3D renderings of the depth map in
Figure~\protect\ref{fig:TBC.BneiInput}d]{\em 3D renderings of the
depth map in Figure~\protect\ref{fig:TBC.BneiInput}d. The spheres in
(a) and (b) indicate the estimated camera location. Note in overhead
view (d) true Euclidean structure is recovered with correct $90^o$
angles.  }
\label{fig:TBC.BneiRes}
\end{figure}


%\section{Experiments with real images}\label{sec:real.images}


%\subsection{Experimental procedure}

%The images were taken with a Pulnix TM9701 $\frac{1}{2}inch$ CCD video
%camera and an $8.5mm$ lens (focal length of 630 $pixels$). Image
%capture was performed using an SGI Indy Workstation at $640 \times 480
%pixel$ resolution.

%In the two image sequences the vertical motion was
%$10mm$ and the horizontal motion $5mm$. In 
%Fig.~ (\ref{fig:lab263.bump.ps}) 
%we created a scene with few strong feature
%points and edges. We also added rotation around the Y axis in
%conjunction with the horizontal motion. The depths in
%(\ref{fig:lab263.bump.ps}) range from $40cm$ to $17cm$. For pure
%translation image motions were on average $15 pixels$. With $4.0
%degrees$ of rotation around the $Y$ axis the horizontal pixels motion
%averaged $50 pixels$. 

%In Fig. (\ref{fig:lab213.bump.ps}) we created a scene which,
%although it has strong edge features, none of the point features
%correspond to true 3D features and the horizontal lines correspond to
%curved lines in 3D. The depths in (\ref{fig:lab263.bump.ps}) range
%from $62cm$ to $31cm$ and image motions were on average $9
%pixels$. The diagonal bar is balanced on the head a bit closer to the
%camera on image left. The front of the right hand cylinder was at the
%same depth as the tip of the nose.

%\subsection{Results}

%The images were processed using 5 levels of coarse to fine processing
%with 4 iterations at each level.  Varying the number of iterations
%from 1 through 4 had no qualitative impact but using a single
%iteration caused a small change in the resulting motion estimates. A
%$9 \times 9$ region was used for the local constant depth fit at all
%levels.

%Fig. (\ref{fig:lab263.bump.ps}e) shows the
%recovered depth map for the first sequence. 
%Fig. (\ref{fig:lab263.bump.ps}f) shows a 3D rendering of the surface. Image
%(\ref{fig:lab263.bump.ps}a) was used as the texture map. The Z values
%have been scaled by 18000. The rendering uses an orthographic
%projection which creates a slight distortion in the images. 

%In order to get smoother and more visually pleasing results we
%performed some smoothing using a membrane interpolation model used for
%the final stage (figures
%\ref{fig:lab263.smooth.ps}a,
%\ref{fig:lab263.smooth.ps}b). The recovered translation estimates were
%qualitatively correct even when rotation around the vertical ($Y$) axis 
%was added to the horizontal motion. 

%\begin{figure}
%          \centerline{ 
%               \hfill 
%                \epsfxsize=1.1in \epsffile{lab4-images/lab262.ps}
%              \hfill 
%                \epsfxsize=1.1in \epsffile{lab4-images/lab263.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (a) \hfill \hfill (b) \hfill 
%              } 
%          \centerline{ 
%               \hfill 
%                \epsfxsize=1.1in \epsffile{lab4-images/lab260.ps}
%              \hfill 
%                \epsfxsize=1.1in \epsffile{lab4-images/lab272.ps}
%              \hfill }
%                    \centerline{ 
%              \hfill (c) \hfill \hfill (d) \hfill 
%              } 
%  %\centerline{ (a)}
%  \centerline{ \epsfxsize=2.3in \epsffile{lab4-images/lab263.depth.ps} }
%  \centerline{ (e)}
%  \centerline{ \epsfxsize=3.0in \epsffile{lab4-images/lab263.bump.ps}}
%  \centerline{ (f)}
%  \caption{\em 
%    The three input images (a,b,c) with pure translation. (d) another
%    image taken at the same location as (b) with $-4.0degrees$ rotation.
%    The estimated depth map (e) and 3D rendering of the surface (f).  Uses
%    $9 \times 9$ region and a local constant depth model. }
%  \label{fig:lab263.bump.ps}
%\end{figure}

%\begin{figure}
%  \centerline{ \epsfxsize=2.8in \epsffile{lab4-images/lab263.smooth.ps}}
%  \centerline{ (a)}
%  \centerline{ \epsfxsize=2.8in \epsffile{lab4-images/lab263.no-tex.ps}}
%  \centerline{ (b)}
%  \centerline{
%    \epsfxsize=1.6in \epsffile[018 180 594 612]{lab4-images/est-rotation.ps}}
%  \centerline{ (c)}
%  \caption{\em 
%    (a)(b)3D rendering of the estimated surface (scaled in Z) with
%    interpolation for smoothing. (c) Estimated rotation as a function of
%    true rotation. Although there is the correct linear relationship there
%    also appears to be a significant scale error. }
%  \label{fig:lab263.smooth.ps}
%\end{figure}

%Fig. (\ref{fig:lab263.smooth.ps}c) shows the recovered rotation
%estimates for true rotations ranging from $-4.0degrees$ to
%$0.8degrees$. Outside this range the program did not
%converge. Negative rotation angles induced image motion that was in
%the same direction as the image motion induced by the
%translation. There seems to be a bad scaling error in the rotation
%estimates.  Fig. (\ref{fig:lab213.bump.ps}d) shows the recovered depth
%map for the second sequence. The results are qualitatively correct. In
%Fig. (\ref{fig:lab213.bump.ps}f) the texture mapping was removed for
%clarity.



%\section{Discussion and future work}
\section{Summary and Future Work}

We have presented a general relationship between the spatio-temporal
derivatives of three frames and the ego-motion parameters of the two
motions. This relationship was derived first for the general case
where the ego-motion model comprises of the 27 coefficients of the
trilinear tensor, and then for the small-motion model of
Longuet-Higgins \& Prazdny. 

On practical grounds, these relationships are model-based brightness
constraints which provide a linear constraint per pixel in the image
--- thereby providing a direct structure and motion estimation that
cuts through the aperture problem and without prior detection of
feature points. The linear constraints of the 27-parameter projective
model are degenerate and cannot lead directly to a unique solution for
the tensor coefficients although a unique solution can be found
using quadratic admissibility constraints. This added complexity led us to
implement in practice  the simpler small-motion model of
Longuet-Higgins \& Prazdny.

The implementation details of these model-based brightness constraints
are important and include four critical elements: (i) embedding of
the computations within a coarse-to-fine (Gaussian pyramid) framework,
(ii) Newton iterations over the brightness constraint equation,
(iii) post-processing of smoothness and surface interpolation for
obtaining visually pleasing results, and  (iv) obtaining both stability and
accuracy by assuming a very short focal length  for the motion
estimation, but using the correct focal length for the image warping. 

The algorithm was tested on a set of challenging real image
situations which contain very little ``good features'' (local regions
with significant variability of gradient direction) which are necessary for
optical flow and discrete point matching algorithms. Yet, we obtained a
faithful dense depth map of the scene due to the fact that the
algorithm is not hindered by the presence of aperture effects.

\subsection{Future Work - Collinear Motion}\label{sec:TBC.future.work}

The method, as described so far, fails when the two motions are in the
same (or opposite) directions. This is a drawback in many applications
(e.g. 3D reconstruction from a monocular image sequence).  We present
here some early research on ways to overcome this problem. The initial
results look promising (\cite{Stein-thesis}).  We will
investigate the pure translation case because if it fails for pure
translation, it will also fail if there is some rotation.

For pure translation:
\be
I''_{t}s^\top t' - I'_{t}s^\top t'' =0
\ee
and 
\bea
k s^\top t'' + I''_t=0 \label{eq:TBC.horn.weldon2.3}\\
k s^\top t' + I'_t=0 \label{eq:TBC.horn.weldon.3} 
\eea
If the second translation vector $t''$ is proportional to the
translation vector $t'$ then equation~(\ref{eq:TBC.horn.weldon2.3}) is
simply a scaled version of equation~(\ref{eq:TBC.horn.weldon.3})
adding no new information and the solution is ill conditioned. Another
way of interpreting equations~(\ref{eq:TBC.horn.weldon.3}) and
~(\ref{eq:TBC.horn.weldon2.3}) is that a motion $t'$ in one direction
will create a change in the image $I'_t$ which exactly the opposite of
the change $I''_t$ induced by a motion $t''$ of equal magnitude and
opposite direction to $t'$.

The LH\&P model which assumes ($\frac{t_z}{Z}<<1$) and was used to
derive equations~(\ref{eq:TBC.horn.weldon.3})
and~(\ref{eq:TBC.horn.weldon2.3}) is of course not exact unless the
motion in the $Z$ direction is zero, but in any case the contribution
to the optical flow due to translation in the $Z$ direction is small
unless the field of view is very wide. Therefore, for small motion,
the LH\&P model will be quite accurate and the equations will be
ill-conditioned.
 
This situation is closely related to the case of reconstruction from
line correspondences. There is a critical combination of a line
configuration and camera motions which can be stated in the following
theorem:

\begin{theorem}{ Critical configuration of lines and motion}\label{thm:critical.config}

Let ${\cal S}$ be a set of lines in 3D which have a common
intersecting line $L$ (i.e., $S\wedge L=0$ for all $S\in{\cal S}$). Let $s$,$s'$,$s''$ be the projections of the line  $S\in{\cal
S}$ in the three views: Image~1, Image~2 and Image~3. Let $o$, $o'$ and
$o''$ be the corresponding camera centers. If  $o$, $o'$ and
$o''$ all lie on a line $S_i \in{\cal S}$ the corresponding lines
$s$,$s'$,$s''$ do not provide a unique solution to the camera
motion. 
\end{theorem}

We currently have no proof, but simulation experiments confirm the
theorem.  This line configuration is the Linear Line Complex (LLC)
which is examined in~\cite{Stein-Shashua-eccv98}.
In~\cite{Stein-Shashua-eccv98} it is shown that for general motion,
the tensor has a linear degeneracy in the case of an LLC, but a unique
solution is obtained by taking into account nonlinear constraints.
Here, for the specific case where the camera motion is along a line
belonging to the LLC, there is a whole family of solutions. The
epipole could lie anywhere on the projection of the line $L$ onto the
image. Thus, there is one degree of uncertainty in addition to the
scale factor ambiguity.

If we take the common to be the line at infinity on the $XY$ plane we
get the particular case where all the lines lie on planes which are
parallel to the image plane, and the translations are collinear and
also parallel to the image plane.  This is a critical condition
according to Theorem~\ref{thm:critical.config} and the direction of
translation in the $XY$ plane cannot be recovered.  In such a
configuration, lines in Image~1 will be parallel to corresponding
lines in Image~2 and Image~3 (assuming no rotation).

How does this relate to our case? By using the first order
approximation of the optical flow constraint equation to get the
point-line-line correspondences, we have created, as an artifact, a
situation where all corresponding lines are parallel. In reality this
is not the case, unless, all the lines came from planes parallel to
the image plane.


%For example, imagine a 3D line $S$ that is not parallel to the image
%plane. Suppose the camera moves in a straight line (collinear motion)
%parallel to the image plane. The orientation of the image of the line
%the $s$
%will change as the camera moves. As a second example imagine part of
%an image
%where the iso-brightness contours are not straight but circular. In
%general, as the camera moves, the image gradient at a point $(x,y)$ in
%the image, again changes in orientation. These changes in orientation of
%the iso-brightness contour were not taken into
%account. 

\subsubsection{A Possible Solution}

The Optical Flow Constraint Equation
\cite{Horn-Schunk81} was based on the first
order Taylor expansion of $I(x,y,t)=I(x',y',t+\delta t)$ where $t$ is
used to denote time (not translation). Keeping the second order terms
of the Taylor expansion results in:
\bea
\lefteqn{I(x+\delta x,y+\delta y,t+\delta t) = } \nonumber \\
&&I(x,y,t) + \delta x I_x + \delta
y I_y + \delta t I_t + \nonumber \\
&&  \frac{1}{2!}\left(\delta x^2 E_{xx} + \delta y^2 E_{yy} + \delta t^2
E_{tt} + 2 \delta x \delta y  E_{xy} + \right. \nonumber \\
&&\hspace{0.2in} \left. 2 \delta x \delta t E_{xt} + 2
\delta y \delta t E_{yt}\right) + e
\eea
Replacing $\delta x = u' \delta t$ and  $\delta y = v' \delta t$ yields:
\bea
\lefteqn{I(x+u \delta t,y+v \delta t,t+\delta t)  =  }\nonumber \\
&&I(x,y,t) + u' \delta t
I_x + v' \delta t I_y +  \delta t I_t + \nonumber \\
&& \frac{1}{2!}\left((u \delta t)^2 E_{xx} + (v' \delta t)^2 E_{yy} + \delta t^2
E_{tt} + \right. \nonumber \\
&& \hspace{0.2in} \left. 
2 u' v' \delta t^2  E_{xy} + 2 u' \delta t^2 E_{xt} + 2
v' \delta t^2 E_{yt}\right) \nonumber \\
&&+ e
\eea
We notice that two of the second order terms are also linear with $u$
and $v$. If we keep those two terms, (at this point, apart from
convenience, we have no formal justification why we should keep those
two terms and not all the 2nd order terms) and apply the constant
brightness assumption: 
$$ 
I(x,y,t)=I(x+\delta x,y+\delta y,t+\delta t)
$$ we get for the first motion:
\be
u'(I_x + I_{xt}\delta t) + v'(I_y + I_{yt}\delta t) + I'_t=0
\ee
Since $I'_x = I_x + I_{xt}\delta t$ and  $I'_y = I_y + I_{yt}\delta t$ this becomes:
\be
u'I'_x + v'I'_y  + I'_t=0
\ee
(The original work by Horn \& Schunk
used: $ u'\frac{(I'_x + I_x)}{2} + v'\frac{(I'_y + I_y)}{2} + I'_t=0
$). Note that $I'_x$ and $I'_y$ are computed at coordinates $(x,y)$ not
$(x',y')$. Similarly for the second motion:
\be
u''I''_x + v''I''_y + I''_t=0
\ee
Likewise $I''_x,I''_y$ are the brightness
gradients in Image~3 at coordinates  $(x,y)$.  
That leads us to the equation:
\be
I''_{t}s'^\top t' - I'_{t}s''^\top t'' =0
\ee
where:
\be
\begin{array}{ccccc}
s' = \left(\begin{array}{c}
  I'_x\\
  I'_y\\
  - x I'_x - y I'_y
\end{array}
\right) &&&&
d'' = \left(\begin{array}{c}
  I''_x\\
  I''_y\\
  - x I''_x - y I''_y
\end{array}
\right) 
\end{array}
\ee
Since now $s' \neq s''$ the equations are better conditioned. This
depends on how much change there was in the images gradients. 
As we have noted, this depends on whether the surface is parallel to
the image plane or not. 

A modification is required to the coarse-to-fine framework and the
iterative refinement which are described in
Section~\ref{sec:TBC.coarse.to.fine} and
Section\ref{sec:iterative.refinement} respectively. For motion
estimation, Image~1 is warped (using forward warping) towards Image~2
and Image~3. For depth estimation Image~2 and Image~3 are warped
towards Image~1.

We have shown an initial solution to the case of collinear motion with
pure translation and have successfully implemented these ideas. It is
not clear how to transfer this to the general motion case and this
remains an open area of study.

\small
\bibliographystyle{myplain}
\bibliography{\bibDir amnonbib,\bibDir vision}

\end{document}