Fast and Accurate Algorithms for Projective Multi-Image Structure from Motion.
AbstractWe describe algorithms for computing projective structure and motion from a multi-image sequence of tracked points. The algorithms are essentially linear, work for any motion of moderate size, and give accuracies similar to those of a maximum-likelihood estimate. They give better results than the Sturm/Triggs factorization approach and are equally fast and they are much faster than bundle adjustment. Our experiments show that the (iterated) Sturm/Triggs approach often fails for linear camera motions. In addition, we study experimentally the common situation where the calibration is fixed and approximately known, comparing the projective versions of our algorithms to mixed projective/Euclidean strategies. We clarify the nature of dominant-plane compensation, showing that it can be considered a small-translation approximation rather than an approximation that the scene is planar. We show that projective algorithms accurately recover the (projected) inverse depths and homographies despite the possibility of transforming the structure and motion by a projective transformation.
Introduction
This paper extends our previous multi-image structure-from-motion (SFM) algorithms [18][26][16] [24] from
the Euclidean to the projective context. Previously, we assumed that the camera calibration was known and
fixed; the new versions of our algorithms handle sequences with varying and unknown calibrations. We also
adapt our algorithms for the common situation where the calibration is measured up to moderate errors and
is known to be fixed.
As in our previous work, we aim for an approach that is fast and accurate on sequences with small
"signal," i.e., with small (translational) image displacements. The small signal makes the reconstruction
problem difficult, so to be effective an algorithm must exploit all the information in the sequence: we
require an intrinsically multi-image approach [19]. 1 On the other hand, the small displacements simplify
the problem of finding correspondence and limit the effects of correspondence outliers on the reconstruction,
since the correspondence errors cannot be large. Motivated by these factors, as in [18][26][16], we disregard
the correspondence problem and present an algorithm that reconstructs from pre-tracked point features over
a large number of images.
See [21][20] for extensions of our approach to a "direct method" ([5][8]) that reconstructs directly from
the image intensities as well as from pre-tracked point and line data.
For small displacements, most scene points are visible in all the images. This makes it possible to use
fast factorization methods, as we do here [18][26][16][31][30]. It also means that one can use multiple images
to eliminate correspondences outliers. We have found that a good way to get reliable feature tracks is to
compute correspondences first for image pairs and then restrict to features that are tracked consistently over
all images. Eliminating features lowers the signal, so an intrinsically multi-image approach, that effectively
exploits the consistent feature tracks over the entire sequence, becomes even more crucial.
As in our past work, we focus on scenes with a large range of depths from the camera (large perspective
effects), which factorization approaches previous to this work could not handle. For shallow-depth scenes,
one can use algorithms such as [31] to complement our approach. Our assumptions of large depth range and
small image displacements imply that the camera motions are small, and we have designed our approaches
to exploit this.
We present two classes of algorithms. One works only for motions that are "truly" general, i.e., with
camera positions that do not lie on a single plane, but we expect it to be more effective for such motions than
the other approach. In [16], and below, we describe how one can determine self-consistently whether or not
1 By an intrinsically multi-frame approach (which can be either batch or recursive), we mean one that reconstructs directly
from many images, rather than first reconstructing from subsets of a few images and combining the results.
the motion is general enough for this approach to work. The second approach works for any moderate-sized
motions, but it is most effective when the camera moves roughly along a line, and we propose it mainly for
this type of motion. However, our experiments show that it also works well when the camera moves on a
plane or makes unrestricted 3D motions.
Unlike the projective factorization approach of Sturm and Triggs [30], our approach needs no initial guess
for the projective depths. It succeeds in situations where the Sturm/Triggs algorithm [30] fails and usually
gives better results. It factorizes a smaller matrix than that used by Sturm and Triggs and, thus, should be
faster on large problems.
Our discussion illustrates a number of interesting theoretical points. We demonstrate the close analogy
between recovering rotations in the Euclidean context and recovering 2D projective transforms (homogra-
phies) in the projective one. We show that compensating for the dominant plane as in [9][27][11] can be
understood as a small-translation approximation rather than an approximation that the scene is planar. We
show that algorithms can recover the (projected) inverse depths and homographies even in projective SFM.
Some of our experiments focus on the common situation where the camera calibration remains fixed over
the sequence and is approximately known. In experiments with general motion, we show that our algorithm
gives a fast computation of the projective structure which approaches the accuracy of the maximum-likelihood
least-squares estimate. A mixed Euclidean/projective strategy does somewhat better than a
completely projective one. In experiments with linear motions, we find that the Euclidean approach with
a slightly wrong calibration does better than the projective one for intermediate translation directions, but
worse for translations close to the image plane or viewing direction. Our results again compare well with
those of a projective maximum-likelihood estimate. Finally, we show that our approach also gives better
results than the Sturm/Triggs approach when the camera moves on a plane.
1.1

Summary

The next section describes notation and preliminary results that we need to derive our algorithms. Section
2.1 shows that one can accurately recover the "projected" inverse depths despite the freedom to alter the
structure by a projective transform. Section 2.2 characterizes the image displacements caused by the camera's
motion and calibration changes, and Section 2.3 characterizes the flows due to infinitesimal 2D projective
transforms, i.e., infinitesimal homographies.
Section 3 describes our general-motion algorithm in fully projective and mixed Euclidean/projective
versions. Results for this algorithm are presented in Section 3.4. Section 4 describes our algorithm specialized
for the linear-translation case. Section 4.3 gives the results of experiments with it, and Section 4.4 explains
how to extend it to handle any translational motion.
Preliminaries
Notation and Definitions. We use MATLAB notation: a semi-colon separates entries in a column vector,
a comma or space separates entries in a row vector, and a colon indicates a range of indices. We denote the
average value of a quantity X by hXi. If M is a matrix (of any size) with some singular values much bigger
than the rest, we refer to the large singular values as the leading singular values, and to the corresponding
singular vectors as the leading singular vectors.
Given a vector V, define [V] 2 as the length-2 vector consisting of the first two components of V. If V is
a 3D point, define the ideal image point corresponding to V by V j [V] 2 =V z . Note that we take the "ideal
image" to have focal length 1. Define the 3D point v 1 corresponding to the 2d point v 1 by v 1 j
\Theta

If M is a 3 \Theta 3 matrix and v 1 is a 2D image point, let M   v 1 denote the image point obtained from v 1
after multiplying by we use the notation v 1 \Theta v 2 to mean
If V is a 3D vector, we have the identity M
\Theta
z
' \Theta
z
(1)
[MV] z
Note that in general
Assume we have N quantities i a indexed by a. We define fig to denote the length-N column vector
whose a-th element is i a . We also use the notation
\Upsilon to indicate a column vector, typically a "long" vector
with length greater than 3. Similarly, if we have quantities  ab indexed by a and b; we let denote the
matrix whose (a; b)-th entry is  ab .
Let there be N p points tracked over N I images, where we label the images by 1g. We
take the zeroth image as the reference image. Define Pm j
to be the m-th 3D point
in the coordinate system of the zeroth image. Let p i
denote the m-th image point in the i-th
image and let
.
Let K i denote the calibration matrix for the i-th image, where we define these with respect to the ideal
image, e.g.,
(neglecting noise). We take the calibration matrices to have the standard upper-diagonal form
22 K i0
Let T i and R i represent the translation and rotation from the reference image 0 to the i-th image, and
We define the motion of the 3D point P under R i and T i by P
z ). We refer to the T 0i as the 3D epipoles, since the epipoles in the
reference image are given by e
2.1 Projective Inverse Depths
In our experiments, we report results for the (projected) inverse depths. We show that one can recover these
accurately in projective SFM.
In the projective context, one can recover the structure only up to a projective transform \Pi; where \Pi is
a 4 \Theta4 matrix. The structure changes under \Pi by:
where the \Delta m are scalars. Dividing by Z 0
m and Zm , we rewrite this as
~
where Pm are the noise-free ideal image coordinates.
Since we adopt the coordinate system of the zeroth image, we only need to consider \Pi that leave the
image points fixed in this image. Then, for generic scenes, the first three rows of \Pi must have the form
\Theta

, where 1 3 is the 3 \Theta 3 identity matrix and 0 3 is a length-3 column of zeros. This implies
Up to small effects of noise, one can accurately recover the inverse depths up to an additive plane and a
multiplicative scale.
Define a N p \Theta N p projection matrix QL annihilating the three vectors f1g ; fxg, fyg. Then (2) implies
one can accurately recover the projection QL
up to scale and the small effects of noise. In our
experiments, we compare the accuracy of our various algorithms in recovering this projection.
2.2 Image Displacements due to Varying Calibration
Consider two images, e.g., images 0 and i; and neglect the noise. We have
Note that
\Theta
(one can eliminate K i from the denominator on the right since
\Theta
1). Let
represent the image point one would get by a pure translation without changing the original calibration
Write the image displacements as the sum of a pure translational piece and a remaining piece that
contains the rotational effects:
where
Recall that T
z
\Theta
z
and one can rewrite d i
Rm as
Tm
Tm
Thus, the d i
Rm represent the displacements due to a pure 2D projective transform K i R
. The form
of (3) leads naturally to our algorithms described below.
2.3 Infinitesimal 2D Projective Transforms
Our algorithm requires characterizing the form of an infinitesimal 2D projective transform or homography.
Suppose that a set of image points p 0
m is given by a small 2D projective transform M of image points pm .
Write the 3 \Theta 3 matrix M as
F I
is the 3 \Theta 3 identity matrix,  is small, F is a 2 \Theta 2 matrix, and I and J are 2 \Theta 1 vectors. For the
2D image displacements,
For each of the 8 parameters in F; I ; J; we define a length-2 vector h (b)
giving the corresponding first-order
displacement due to the infinitesimal homography. Thus
xm ym xm ym y 2
Define the 8 \Theta 2 matrix hm such that the b-th column of h equals h (b)
.
2.4 Additional Definitions for the Algorithms
The homography flow vectors
\Psi (b) . Assuming N p image points, we define eight length-2N p "homog-
raphy flow" vectors
\Psi (b) corresponding to the h (b)
defined above, with
fxyg \Phi
Define the 2N p \Theta 8 matrix \Psi such that its b-th column equals
\Psi (b) .
Translational flow vectors
\Phi (b) . Define the three length-2N p translational-flow vectors
\Phi z ,
corresponding to the translational displacements in (3), by
\Phi
f0g
\Phi
Let
Image displacement matrix. We organize the observed image displacements d i
m into a (N I \Gamma 1) \Theta 2N p
matrix D; by putting all the x and then the y coordinates for a given image on a single row.
Bias Compensation Matrix C. Define a
We use
F
Note that one can compute products of C \Gamma1=2 and C with O (N) computation.
Modified translation vectors
Ca . For a 2 fx;
a g be the length-(N I \Gamma 1) vector whose
elements are the T 0i
a . Define
a
3 General-Motion Algorithm
Assumptions. As in the approach of [16], the algorithm requires that the translational motion not be too
large (e.g., with jT 0 j =Z  1=3) and that the camera positions do not lie in a plane. One can automatically
detect when the camera is moving on a plane or line and use the approach of Algorithm II below [26].
3.1 Algorithm (Proj)
Step P1: Homography compensation. Assuming that the translations are zero, we recover the homographies
separately between the reference image and each subsequent image i. We
compensate for these homographies, defining the compensated image i by p i
. Let the
image displacements d i
m and displacement matrix D now refer to the compensated image points p i
cm .
Step P2a: Homography elimination. Using Householder matrices [25][4], we compute a (2N
annihilates the subspace generated by the 8 vectors
projection (i.e., the rows of H are orthogonal).
Step P2b: Singular value decomposition. We define the modified displacement matrix
and compute its singular-value decomposition (SVD). Let
A (1;2;3) be the three leading right singular vectors
of DCH , and let A j
\Theta
A (2) ;
A (3)

Step P3a: Depth recovery. We recover the depths by solving the linear system
\Theta
\Gamma\Phi
\Phi y
\Gamma\Phi
\Phi z
\Gamma\Phi
for the Z \Gamma1
m and the 3 \Theta 3 matrix U .
Step P3b: Translation recovery. Using the recovered values for the
\Phi a , we solve the linear system
Cx
Cy
Cz
\Phi y
\Phi z
for the
Ca . The fT 0
a g are given by fT 0
a
Ca .
Step P4: Improved homography recovery. Let the T 0i
represent the current estimates of the T 0i .
For each i; we solve the linear system
\Theta d i
\Theta I i
\Theta pmA J iT pm ; (10)
for I i , F i , J i . Our estimate of the residual homography from the reference image to image i is
Step P5: Iteration (optional). We iteratively repeat until convergence the following: 1) Compensate
for the residual homographies newly computed in Step P4; 2) Repeat Steps P2b-P4.
3.2 Discussion
Step P1: Homography compensation. This is a preprocessing step and sometimes unnecessary. Its
purpose is to make the effective motion smaller. One could also apply this step in the Sturm/Triggs approach
[30], since this approach, like ours, works best for small motions.
The errors in recovering the M i scale with the size of the translational image displacements,
O
O
where
fi fi in the correction term denote the average sizes of the image noise, inverse depths, and
translations. Under our moderate-translation assumption, the
fi fi are small.
Compensating for the M i is equivalent to compensating the dominant plane 2 [9] [27][11]. Our formulation
makes it clear that one can understand this compensation as a small-translation approximation rather than
as an approximation that the scene is planar.
2 Though the initial step of our method is the same as in [9], our subsequent algorithm differs in two important respects: it
is linear, and it corrects for the first-order errors in the compensation of the dominant plane.
Steps P2a and b: Homography elimination and SVD. One can compute H and its products with
O (N p ) computation [16]. Our current implementation computes the products of H in O(N 2
LAB's overhead makes our O (N p ) implementation quite slow.
After the homography compensation in Step P1, we have
z
\Theta
z
where
We neglect the denominator in (11) and D, since it causes a second order correction O
. Multiplying
D by H eliminates the first-order corrections due to ffiM i , so up to second order we get a bilinear
expression for
Cx
Cy
Cz
z
Step 2b is based on approximating DCH as bilinear in the structure and motion.
Our use of H to annihilate the residual homographies derives from the optical-flow technique of Jepson
and Heeger [10], as generalized by [25] to apply to arbitrary, rather than regular, image-point configurations.
For sideways translations, the denominator in (11) exactly equals 1. If we can succeed in making the
residual homography ffiM i small, for example by means of the iteration in Step P5, then DCH becomes
exactly bilinear for sideways translations, no matter how big. Thus, our algorithm is capable of giving good
results for large sideways translations. 3 The experiments of [14] with the Euclidean version of our approach
show that iterating typically does make the residual rotations small, even with large translations. Rotations
are the Euclidean equivalent of homographies, so this implies that iterating will also make the ffiM i small,
and that our approach will succeed for large sideways translations. We have confirmed this on a few test
sequences, see below. However, our approach is intended mainly for small or moderate translations.
As discussed in [16], multiplying by C \Gamma1=2 reduces the bias due to singling out the reference image for
special treatment.
This algorithm requires that DCH has 3 singular values that are much larger than the rest. This is the
precise form in which we impose our general-motion assumption. If DCH does not have 3 large singular
values for a given sequence, one should use Algorithm II below rather than the approach of this section.
3 We presented an algorithm specialized for sideways motions in [26].
Step P3a: Depth recovery. (8) reflects the fact that the recovered right singular vectors
A (b) must
generate approximately the same subspace as the three translational flow vectors
\Phi x;y;z .
As discussed in Section 2.1, one can recover the Z \Gamma1
only up to an additive plane and multiplicative
scale. One can see this explicitly from the fact that H annihilates all contributions to the
\Phi x;y;z from the
components of
are length-N p vectors.
One can solve (8) with O (N p ) computation, as discussed in [16]. In our current implementation, we use
a simple O(N 3
Step P4. Improved homography recovery. From (11), we get
\Theta d i
\Theta
\Theta
where
denotes the size of the error in estimating the T 0i and
denotes the average size of T 0i . As
before,
fi fi and j denote the average sizes of the inverse depths and the noise. Define I i ; F i , J i as in (4)
by
Then (13) and (5) give the linear system (10).
The linear system (10) does not determine the residual homographies ffiM i unambiguously. This just
reflects the projective covariance of the equations: one is free to change the projective basis, which changes
the values of the ffiM i , Z \Gamma1
, and the T 0i . As discussed in Section 2.1, one can characterize the projective
transforms that leave the reference image fixed by their effects on the Z \Gamma1
m . We remove the ambiguity in
recovering the ffiM i by setting to zero the part of ffiM i that would add a plane to the inverse depths.
Remark. One could also estimate improved values of the image points as in Step 4 of the algorithm of
[16], but we have not implemented this.
3.3 Variations
The algorithm we have presented is fully projective: it can handle arbitrary changes in the linear calibrations
images. But, most sequences are taken with a single camera, so that the calibration remains
essentially constant over the sequence (with the possible exception of the focal length). We have considered
a few simple modifications of Step P1 that exploit this. One variation, the fixed algorithm, recovers
the homographies between the reference and subsequent images by a least-squares optimization under the
assumption that the calibration K is fixed. (Like the original Step P1, it assumes that the translations
are zero.) The other variations exploit the fact that the calibration error is typically small, so that one
can neglect the homographies K i R(K) \Gamma1 or approximate them by pure rotations R i . The Proj-nocomp
variation of our algorithm does not compensate for the homographies at all, simply eliminating Step P1.
The version Proj-Unrot approximates K by a pure rotation. Instead of Step P1, it computes
and compensates for the best rotations transforming the reference image to the subsequent images [16].
Finally, we have compared our projective algorithms to a Euclidean version of the same algorithm, which
we described in [16].
3.4 General-Motion Experiments
In the following synthetic experiments, the motion, structure, and image noise varied randomly for each
trial. The sequences consisted of 15 images of points with a field of view (FOV). The 3D depths
varied from 20  Z  100. In Experiments 1-4 we randomly chose each translation component (with
respect to the zeroth camera position) such that \GammaT max  T x;y;z  Tmax and added random rotations up
to a maximum of about 20 ffi . In Experiment 5, the motion consists of a scene rotation by up to about
(This corresponds to very large camera translations in camera-centered coordinates and poses a serious
challenge for our approach, since it is targeted for moderate translations. Other algorithms, e.g., the two-image
or Tomasi/Kanade approach, are more appropriate for this situation [19][31].) The noise was one-pixel
Gaussian, assuming a 512 \Theta 512 image and the specified FOV.
We simulated calibration error by multiplying the images derived as above by a matrix K; with
This corresponds to a calibration error in the focal length of 5-10%.

Table

1 shows the results for several versions of our algorithm in comparison to the maximum-likelihood
estimates. The entries give the mean error in degrees over all trials between the recovered and ground-truth
values of QL
(Section 2.1). Proj refers to our algorithm of Section 3.1, and P-U refers to
Proj-Unrot described above. MLE gives the results for the maximum-likelihood least-squares estimate,
computed by a standard Levenberg-Marquardt steepest-descent approach starting from the ground truth.
Expt #Seqs Proj P-U Euc MLE Tmax

Table

1: Experiments 1-5: general motion. Average errors in degrees in the projected inverse depths
1g.
Euc is our Euclidean algorithm [16].
Proj-Unrot gives the best of the tabulated results. The fixed results (not shown) are comparable, but
not good enough to justify their extra computational cost. For small translations, with our
algorithm should be most accurate, Proj-Unrot does only 11% worse than the MLE. Euc also does well.
For its results are only 12% worse than those of Proj-Unrot. One of the reasons for the good
performance of Proj-Unrot and Euc is that the calibration "error," i.e., the difference between K and
the identity matrix, is moderate.
As expected, our algorithm does relatively less well than the MLE in Experiment 5. However, because of
the large baselines in this experiment, the resulting translational image displacements constrain the reconstruction
so strongly that our algorithms still give good results, i.e., Proj-Unrot averages better than 7
error.
We also computed the results (not shown) for another variation of our algorithm, where instead of Step
P1 we computed all the homographies between the reference and subsequent images in a single optimization.
(One way of doing this is to use the 2D version of the Sturm/Triggs approach [30][15].) Since Step P1
computes the transform between image 0 and image i separately for each i; it overweights the reference
image noise, but by computing all the transforms simultaneously we can avoid this. However, we found that
the single-optimization approach did not improve our results.
Real-Image Sequence. We tested our algorithm on the Castle real image sequence available from CMU.

Figure

1 shows the first image of this sequence. We generated the images for this experiment by multiplying
the correctly calibrated images (for unit focal length) by 12:2 and then shifting by (0:403; 0:093) (to center
the images and scale them to unity). Thus the assumed focal length was incorrect by a factor of 12:2. Our
pure projective approach Proj gave an error for the projected inverse depths of :40 ffi , compared to :04 ffi for
the MLE. Euc gave an error of 28:2 ffi . The largeness of the Euclidean error is due to the very large error in
the assumed focal length.
Since the scene in this sequence has a very shallow depth, an affine approach such as the Tomasi/Kanade

Figure

1: Castle image
algorithm [31] is more appropriate than ours [16].
Comparison to the Sturm/Triggs Approach. The Sturm/Triggs algorithm [30] is a factorization
method that, like ours, deals well with small motions and large perspective effects. We compared our
approach to this algorithm, using an implementation of it that we created previously for other purposes.
We have optimized the code for the Sturm/Triggs algorithm to about the same extent as for ours. In the
Sturm/Triggs implementation, we followed the advice of [32]: before applying the algorithm, we first centered
and scaled each image to a unit box, and then normalized each homogenous image point to unit norm. We
initialized the algorithm by setting all the unknown projective depths equal to one [32][1][7], as is appropriate
for small motions. (Since we aim for a true multi-image technique that works for low signal-to-noise, we do
not compute the projective depths from two or three images as in the original algorithm of [30].) We have
proven in [15] that the iterative version of the Sturm/Triggs algorithm converges, and we also tested our
approach against the iterated version. The implementation of Steps P2 and P3 of our algorithm requires
about 80 lines of MATLAB code, and Steps P4 and P5 require an additional 30. The iterative version of the
Sturm/Triggs approach requires about 80 lines.
We created synthetic sequences using the ground-truth structure from two real-image sequences: the
UMASS/Martin-Marietta rocket-field sequence [3] and the UMASS PUMA sequence [12]. The points in
the rocket-field sequence range from 17 to 67 in depth and cover an effective FOV of 37 ffi , while the PUMA
points range from 13 to 32 in depth and cover an effective FOV of 33 ffi . Table 2 gives the parameters of the
experiments, and Table 3 shows the errors in recovering the projected inverse depths, the epipoles in the
zeroth image, and the homographies. As Table 2 indicates, in most experiments each image had a different,
randomly chosen calibration. Define r (a) to be a random number chosen uniformly in the interval [a; \Gammaa].
For each of the N I images, we selected the calibration matrix via
\Theta
\Theta

\Theta
\Theta
\Theta
\Theta
chose r separately for each entry of the matrix.)

Table

2 gives the size of the added noise in pixels, where we define the size of a pixel by taking the
maximum magnitude of the image point coordinates to correspond to 256 pixels. Note that this is after
applying the calibration matrix K: Since the shift in the camera center due to
\Theta

1:2;3
displaces the image
region from the origin, this noise is usually larger than it would be for an image region centered on the origin.
We define the epipole error for a sequence as the average over images
true
true
calc and e i
true are the calculated and true epipoles in the reference image. One can show, as in Section
2.1, that projective transforms that leave the reference image fixed also leave these epipoles fixed.
We obtain the homography error for a sequence by averaging the homography error over all images
with i  1. By the "homography error," we mean the error in recovering the K projective
transform leaves the reference image fixed, one can show that it changes this matrix by
true
(up to noise), where the V i are length-3 vectors. If G i is a 3 \Theta 3 homography matrix for the i-th image,
define the corresponding invariant homography matrix by
invar
where we use the backslash notation of MATLAB, i.e., e i
true is the matrix division of e i
true into G i .
With this definition, invar
is invariant (up to noise) to projective transforms that leave the reference
image fixed.
As before, we denote the recovered homographies by M i . We define the homography error for image
i by the angle in degrees between the length-9 vectors V true and V calc , where these contain the entries,
respectively, of invar
and invar
Expt. Struct Uncalib oe T oe ROT Noise S/T Proj
deg. pix. itermax itermax

Table

2: Parameters for the experiments comparing Proj and the Sturm/Triggs approach (S/T). The 'Struct'
column indicates the structure used to generate the sequences. A 1 in the 'Uncalib' column indicates that we
introduced different calibrations for each image (see text). We choose T x , T y , T z independently as Gaussian
variables. The columns labelled oe T tabulate the standard deviations of the T x;y;z : oe ROT characterizes the
typical size of the rotations. The 'noise' column indicates the standard deviation of the Gaussian noise
added. The remaining two columns indicate the maximum number of iterations allowed, respectively, for
the Sturm/Triggs approach and for Step P5 of the Proj algorithm.
Our results show that our algorithm Proj usually gives better results for the structure and motion
than the Sturm/Triggs approach. One iteration of Proj takes slightly longer than one iteration of the
Sturm/Triggs approach (recall that our current implementation is slower than necessary). However, the
preprocessing Step P1, that is, the initial homography recovery, accounts for most of the computation time.

Table

3 shows that the factorization part of our algorithm, Steps P2-P3, is about four times faster than the
Sturm/Triggs approach. Also, the computation times of Experiment 7 suggest that our approach converges
more quickly than the Sturm/Triggs method; see also Section 4.4.
We also checked the performance of our algorithm on a sequence with large sideways translations. With
the PUMA structure, oe (refer to Table 2 and the
definition of r (a) above), our algorithm gave errors for the projected depths of 1:6 sequence, the
largest translation had size 13.4, which is comparable to 13:9, the distance of the closest 3D point to the
camera's reference position.
4 Algorithm II
The algorithm presented above assumes that the translational motion is sufficiently general, i.e., that the
camera locations do not lie close to any line or plane. In this section, we present a version of our algorithm
that deals with the common case of a camera moving along a line. One can extend it to deal with more
Expt. Time (sec) Z \Gamma1 Epis Epis Homog Homog Cycles
9 0.042 0.19 0.15 8.6 8.4 8.1 2.7 6.9 1.8 5.4 1
9(med) 0.04 0.19 0.15 6.8 6.7 7.1 1.1 2.8 1.1 4.9 1
4.3 3.2 14 1
11(med) 0.04 0.19 0.15 9.1 11 11 1.4 2.8 2.3 5.2 1
19(med)

Table

3: Results for the experiments with parameters given in the previous table. 'NoComp' and `NC'
refer to Steps P2 and P3 of the Proj algorithm, without any computation or initial compensation for the
homographies. The columns labelled 'S/T' give results for the Sturm/Triggs approach. We present errors
for the projected inverse depths ('Z \Gamma1 '), the epipoles ('Epis'), and the homographies ('Homog'). See text
for an explanation of the error measures used. The 'Cycles' column indicates the number of cycles used by
Proj in Step P5. All quantities shown represent the mean of the results over 100 random trials, except for
those rows labelled by '(med)', where the quantities represent the median over 100 random trials.
general motions such as a camera moving on a plane or in 3D, as we describe briefly below.
Definitions. Let the unit vector T 0 denote the direction of the T 0i , and let  i j
Let
\Theta d
denote the (N I \Gamma 1) \Theta N p matrix whose (i; m)-th element equals (pm \GammaT 0 ) \Theta d i
be the N p \Theta N p projection matrix that annihilates the subspace generated by the eight
vectors
ae
\Theta h (b)
oe
One can show that this subspace is five-dimensional.
4.1 Algorithm Description
Step L0: Rescaling. We transform all images so that they center on the origin and have the same scale.
We choose the scale so
for the reference image. (At the end of the algorithm, we
transform the recovered unknowns back to the coordinate system of the original reference image.)
Step L1: Homography compensation. Assuming that the translations are zero, we recover the
homographies separately between the reference image and each subsequent image i.
Define the compensated image i by p i
. Let the image displacements d i
m and displacement
matrix D now refer to the compensated image points p i
cm .
Step L2a: Using Householder matrices [25][4], we compute a (N
annihilates the subspace generated by the six length-N p vectors f1g
is a projection.
Step L2b: Linear translation recovery. We solve the linear system
for T 0 .
Step L3: Iterative improvement of recovered translation. Starting from the previous estimate
we minimize the error
trace
\Theta d
\Theta d
with respect to T 0 . Take T
Step L4: Improved homography recovery. The same as in Step P4 of our general-motion approach.
Step L5: Iteration (optional). The same as in Step P5 of our general-motion approach.
4.2 Discussion
Step L0. As in [6] and [18][13], this step reduces the bias of our linear algorithm in Step L2.
Step L2. This step of our projective algorithm is exactly the same as in our Euclidean algorithm of [18];
the only difference is that the projective algorithm computes an estimate of KT rather than T: However,
in the projective case, the linear algorithm exploits more of the available information and gives a better
approximation to the result of minimizing the full error. Since HL annihilates only one more degree of
freedom than Q 5 ; the projective linear algorithm uses all but one of the available constraints, while the
Euclidean linear algorithm forgoes three of the available constraints. One can show that the additional
length-N p vector annihilated by HL is
z
z
For linear motion, (3) becomes
z
Assume we have compensated for the homographies in Step L1. At the ground-truth value for
to noise, and
independent of the denominator in (16). Thus, if the residual homographies ffiM i are small, our algorithm
works for any size and direction of translation [18]. This also holds for Step L3. If Step L1 does not suffice
to make the ffiM i small, the iteration in Step L5 may. Thus, we expect the algorithm to work even for very
large translations, as the Euclidean version of our approach has been shown to do [14].
Step L3. For e =T 0 in or near the FOV, instead of Steps L2 and L3, one could use the recent method of
Srinivasan [29][28] to find the true global minimum of (15) directly. This is important, since [14][2] showed
that the least-squares error typically has several local minima for e near the image region. We have not yet
implemented Srinivasan's approach.
Remark. One could also recover the Z \Gamma1
m as in [18][23] and improve the homography computation as in
Step P4 of the general-motion algorithm. Since these are straightforward transcriptions of previous methods,
we do not discuss them here.
4.3 Linear-Translation Experiments
We generated sequences as in Experiments 1-4 except that the camera translated in constant steps of 0:2
along a line (with random rotations as before). We created 400 sequences, choosing random directions for
the image-plane projections of "
T and systematically varying T
y from :01 to 4 in steps of :01.

Figure

2 shows the mean angular errors in recovering T 0 for three versions of our algorithms: a pure Euclidean
approach [18][23], the projective approach Algorithm II, and Algorithm II with Step L1 replaced
by initial compensation of the rotations. We derived the curve labeled MLE by first using Levenberg-Marquardt
to compute the maximum-likelihood least-squares projective reconstruction and then computing
the calibration and motions from this by least-squares minimization.
The Euclidean estimate is worse than the projective ones when T is parallel to the image plane or to " z,
but it is comparable 4 for intermediate T. The overall median of these results is slightly better than those of
4 We verified that the Euclidean approach gives bad results for T ? "z simply because of the calibration error. With no
calibration error, it does perform well for T parallel to the image plane.
(degrees)
y

Figure

2: Angular errors in recovering T 0 . x axis shows the lower limit of a bin of size 0:25 in T
each data point is an average over 25 trials from the indicated bin. Cyan circles: MLE; blue squares:
Euclidean. x and + are for our linear projective algorithm, with initial compensation for the rotations and
homographies, respectively.
the projective algorithms.
We also ran Algorithm II (with initial projective compensation) on all 55 choices of image pairs from
the Castle sequence (Figure 1). Our algorithm recovered T 0 with an average error of 1:1 ffi , compared to 0:71 ffi
for the MLE. The median errors for the two approaches were respectively :80 ffi and :40 ffi , and the maximum
errors were 6:6 ffi and 4:4 ffi . Though the assumed focal length was incorrect by a factor of 12:2 (Section 3.4),
our approach did nearly as well as the MLE.
Because of the large error in the focal length, the variation of initially compensating for a rotation instead
of a homography yielded relatively large errors of about 10 ffi in the rotation. Our Euclidean algorithm usually
recovered T 0 accurately. Apart from 13 outliers (possibly due to local minima), the Euclidean algorithm
recovered T 0 with an average error of 5:1 ffi .
Comparison to the Sturm/Triggs Approach. To compare Algorithm II to the Sturm/Triggs ap-
proach, we created synthetic sequences of 16 images and 32 points using the structure from the PUMA
sequence. Let "
T true denote the true translation direction and "z the viewing direction in the reference image.
We systematically varied the angle ' true ji
from 0 ffi to 86 ffi in increments of 2:9 ffi . For each
selected ' true , we created sequences, where we chose the projection of "
T true on the x-y plane randomly
for each sequence. For each image, we randomly chose the rotation R i with oe with a standard
deviation of 5 ffi for the rotation around each axis, and we randomly and uniformly chose the translation
magnitudes up to a maximum of 1. Recall that the PUMA depths range from 13-32. We introduced varying
Angular
error
in
recovered
epipole
(deg.)

Figure

3: The angular errors in recovering the 3D epipole T 0 for the Sturm/Triggs approach and Algorithm
II, plotted versus ' true
trials. 'Triangles' show the
result of one Sturm/Triggs iteration; 'diamonds' show the Sturm/Triggs result after convergence or after a
cut off of 300 iterations if the algorithm did not converge by this number. '+' shows the linear estimate of
Algorithm II, Step L2, and '*' shows the result after the nonlinear minimization in Step L3 of Algorithm II.
calibrations as in Experiments 6-19, and added Gaussian noise of 0:05 pixels, assuming that the maximum
magnitude of the image point coordinates corresponded to 256 pixels.

Figure

3 compares the results of Algorithm II and the Sturm/Triggs approach for the 3D epipoles.
Each data point represents the mean over the 30 trials for the indicated value of ' true =i
: For
our approach, we plot results for the error
calc
true is the true value of the 3D epipole and
calc is the value recovered. Since the Sturm/Triggs approach recovers a different 3D epipole T 0i
S=T for each
image, we plot the average error
ii
our approach, we show the errors for
the initial linear estimate Step L2 and for one iteration of the nonlinear estimate Step L3. For comparison,
we show the Sturm/Triggs estimates after one iteration and after the algorithm either converges or reaches
300 iterations.
One iteration of the Sturm/Triggs approach does much worse than our linear estimate, and the converged
Sturm/Triggs result is also much worse. Our linear estimate is almost as good as the nonlinear estimate.
Our nonlinear algorithm averaged about 2 seconds of computation and the linear algorithm took fractions
of a second, while the iterated Sturm/Triggs approach averaged about 9 seconds.
We also applied the iteration of Step L5 and allowed Algorithm II to converge. The average over all

Figure

4: The first image of the rocket-field sequence.
sequences of the number of cycles till convergence was ! 4, and the average time till convergence was 8
seconds, less than the 9 seconds for the Sturm/Triggs approach. The Sturm/Triggs approach averaged 177
iterations. However, we used a much stricter convergence criterion for the Sturm/Triggs approach, to give
it a chance to converge to an accurate result. As proven in [15], the Sturm/Triggs approach minimizes
a particular error function. We defined the algorithm to have converged when that error changed by less
iterations. We defined Algorithm II to have converged when the residual homography
recovered in Step L4 satisfied
For an additional comparison, we ran the Step L5 iteration of Algorithm II on 100 sequences similar
to those above, defining convergence by
On average, Step L5 converged in 7.5
cycles and always in less than 15.

Figure

5 shows similar results for the homography recovery, for the homography error defined in Section
3.4.
We also tested Algorithm II and the Sturm/Triggs approach on the rocket-field real-image sequence [3],
see

Figure

4. This sequence has large translations ranging up to 7:5 in magnitude (recall that the depths vary
from 17 to 67). The camera moves approximately along a line, with T i deviating from its average direction
by up to 1:5 ffi . Steps L2 and L3 of our approach recovered T 0 with errors of 5:8 ffi and 6:4 ffi , respectively, and
Steps L1 and L4 recovered the homographies with average errors of 2:8 ffi and 2:0 ffi . The initial Sturm/Triggs
estimate had average errors of 14:7 ffi and 15:5 ffi for T 0 and the homographies, respectively. After 100 iterations,
the Sturm/Triggs approach gave errors of 3:3 ffi and 7:1
We have found that the Sturm/Triggs approach gives poor results for linear motion when the motion is
forward (or backward). We created sequences with 16 images of 32 points for randomly chosen structures,
where the depths varied between 40 and 60, the translation direction was (0:1; 0:1; 1) ; the
randomly up to 1, the rotations had oe there was zero noise. We allowed the Sturm/Triggs
approach up to 3000 iterations to converge, and defined it to have converged when its error changed by less
iterations. We proved in [15] that the Sturm/Triggs algorithm does eventually converge
to a local minimum of the error. On 9 of the trials, the algorithm had large errors in recovering the
epipoles, with an average error of 44 ffi for these trials. The maximum error of our linear estimate in Step L2
was 0:15 ffi . For a similar set of sequences with rotations of up to 15 ffi , the Sturm/Triggs approach gave
large errors on 25 out of our approach again gave nearly perfect results. We also found that
the Sturm/Triggs approach often failed to converge correctly when we used the structure from the PUMA
sequence to generate sequences.
These failures are produced by a bad initial choice of the projective depths, which causes the Sturm/Triggs
approach to converge to an incorrect local minima of its error. But, for small forward or backward motions,
there is no way to get good initial estimates for the projective depths-one cannot compute them accurately
from any small number of images, as [30] proposed to do. Also, the projective error surface is flat for
epipoles away from the forward-motion direction [17], and it tends to have local minima near the image
points [14][2][17]. If the initial choice of projective depths corresponds to an epipole far from the image
region, these factors could cause a local minimum to intercept the algorithm and prevent it from converging
to the true global minimum.
4.4 Extension to General Motion
We describe how to extend Algorithm II to any motion [26].
Step L'0. Rescaling.
Step L'1: Homography compensation.
Step L'2: Define H and HL and compute the SVD of DCH . Let NS  3 be the number of large
singular values of DCH , and let
A
be the right and left singular vectors corresponding to these
singular values. Let oe (s) denote the leading singular values of DCH , and define the (2N
oe (1)
Step L'0 is the same as Step L0, and Steps L'1-L'2 are the same as Steps P1, P2 in the general-motion
algorithm.
Homography
errors

Figure

5: The errors in recovering the homographies for the Sturm/Triggs approach and Algorithm II,
plotted versus ' true
all images in
each sequence. 'Triangles' show the result of one Sturm/Triggs iteration; `diamonds' show the Sturm/Triggs
result after convergence or after a cut off of 300 iterations if the algorithm did not converge by this number.
'+' shows the initial estimate of Algorithm II, Step L1, and `*' shows the result after Step L5 of Algorithm
II has converged.
Step L'3a. For each s with 1  s  NS ; define
d
A
A
Np+m
. That is, we define the  d
m so that H T
hn
d
x
d
y
Step L'3b: Linear estimates of effective translation directions. Let the unit vector T 0(s) be the
translation direction corresponding to
A . For each s, we solve
ni
\Theta
d
in the least-squares sense for T 0(s) .
Step L'3c: Refinement of effective translation directions. For each s, we refine the estimate for
T 0(s) by minimizing
d
with respect to T 0(s) .
Step L'4a: Isolate translational flow. Recall that hm is a 8 \Theta 2 matrix such that the b-th column of
hm equals h (b)
m . For each s, we compute a length-8 vector w (s) by solving
pm \GammaT 0(s)
hmw
Compute the
\Delta from H T
Step L'4b: Depth recovery. Solve the linear system
z
oe
for the Z \Gamma1
m and the 4NS constants
Step L'4c: Full translation recovery. We recover the translations via
\Theta \Phi
x
y
z
Step L'5: Improved homography recovery. The same as in Steps L4 and P4 of the algorithms
described above.
Step L'6: Iteration (optional). The same as in Steps P5 and L5 of the algorithms described above.
4.4.1 Discussion
This algorithm is the projective version of our Euclidean algorithm in [26].
Step L'3: Translation-direction recovery. Each leading singular vector corresponds to an effective
translation direction T 0(s) . We recover the T 0(s) exactly as before in Algorithm II.
Step L'4a,b: Isolate translational flow; depth recovery. For each leading singular vector, the hmw
is the effective infinitesimal homography. As before, one cannot determine this uniquely due to projective
covariance, and we specify it by setting to zero the components of w (s) that would add a plane to the Z \Gamma1
.
The constants a are necessary since Step L'4a recovers the homographies up to an ambiguity.
If s 0 6= s; the Z \Gamma1
m corresponding to the recovered
may differ by an additive plane from the Z \Gamma1
We need the scales  (s) to fix the scale of the Z \Gamma1
between different singular vectors.
We introduce the singular value oe in (18) to emphasize the singular vectors with larger singular values,
since these have less noise sensitivity.
4.4.2 Experiments
We created 100 sequences with general motion as in Experiments 6-19, using the PUMA structure, varying
calibration, oe (refer to Table 2). We ran the extended version of
Algorithm II without the iterative refinement of translations in Step L'3c and without the iteration of Step
L'6. Without Steps L'3c and L'6, the extended algorithm is purely linear, with no iteration apart from the
SVD and linear equation solving. We compared our results to those of the Sturm/Triggs (S/T) algorithm,
which we allowed up to 300 iterations to converge. On average, S/T required 135 iterations to converge.
This noniterative version of our algorithm averaged 0:21 seconds of computation, compared to 7:0 seconds
for the iterated Sturm/Triggs (IST) approach. It gave an average error for the projected Z \Gamma1
m of 1:39 ffi
compared to 1:63 ffi for IST. Its average error in computing the epipolar directions was 8:8
5:3 ffi . Its average homography error was 0:16 ffi and that of the IST was 0:52 ffi .
A single iteration of S/T gave an average error of 29:7 ffi for the epipoles and 3:2 ffi for the homographies.
We created a second set of 100 sequences with planar motion with the same parameters as the first.
We achieved planar motions by creating the T i as before and then setting the third singular value of the
matrix of translations to zero. The average errors for the projected Z \Gamma1
were 1:62
The average errors for the epipolar directions were 4:4 (IST). The average errors for the
homographies were 0:11
On average, IST took 140 iterations to converge. The first iteration of S/T gave an average error of 33 ffi
for the epipolar directions and 3:7 ffi for the homographies. We also ran IST for 11 iterations, to check its
speed of convergence. This took on average 0:61 seconds and gave average errors of 5:1 ffi , 22:9 ffi , and 2:4 ffi for
the projected inverse depths, epipolar directions, and homographies, respectively.
4.4.3 Approximately Linear or Planar Motions
What happens when we apply Algorithm II or its extension with a too restrictive motion model-for
example, assuming planar motion when the camera centers do not lie exactly on a plane? To the extent that
we can neglect the denominator in (3)-and doing so causes small second order effects when the translations
are small-this causes no additional error in the structure recovery. "Applying our algorithm with a too
restrictive motion model" means using NS of the leading singular vectors when the dimensionality of the
translational motion is greater than NS , i.e., we neglect the smaller components of the translational motion.
This does not prevent us from recovering a good approximation to the translations projected into an NS -
dimensional subspace. Since each singular vector gives a separate estimate of the Z \Gamma1
m , the only effect on
the structure recovery of restricting to a subset of the singular vectors is that we lose information that could
have been used to improve the estimates.
We ran our extended Algorithm II with sequences created as in our above test of this
algorithm on 100 general-motion sequences. As before, we ran the algorithm without the iterative refinement
of translations in Step L'3c and without the iteration of Step L'6. We obtained an average error for the
projected inverse depths of 1:37 ffi , which is less than the result with than the IST error.
We projected the true 3D epipoles into the plane of the epipolar directions recovered by our algorithm. The
average error in recovering the 3D epipoles in this plane was 6:4 ffi . The average error of Step L'5 in recovering
the homographies was 0:20 ffi .
These results suggest that it can be advantageous to ignore the smaller of the leading singular values
and vectors when these have a good deal of noise contamination. In the sequences above, the x, y, and
z components of the translational motion have similar magnitudes, but the smallness of the FOV causes
the signal (i.e., the image displacements) from the z-translations to be much smaller than for the other
directions. One can detect this from the smallness of the third leading singular value of DCH , which for
these sequences is often not much bigger than the noise effects.
Once one has recovered the structure accurately using a subset of the leading singular vectors, one can
return to calculate the translational components corresponding to the neglected leading singular vectors.
Conclusions
We presented fast projective structure-from-motion algorithms which have comparable accuracy to the MLE
and appear superior to the Sturm/Triggs approach [30]. These algorithms work for any motion as long as
the camera displacements are not too big, with jTj =Z ! 1=3 [16]. We showed experimentally that the
Sturm/Triggs approach often fails for linear camera motions, especially forward or backward motions. We
speculated that the recent characterizations of the projective least-squares error surface in [17][2][14] may
offer part of the explanation for this. We studied the advantages of a pure projective approach, versus
a mixed Euclidean/projective strategy, for the common situation when the calibration is fixed and partly
known. We showed that algorithms can recover the (projected) inverse depths and homographies even in
projective SFM. We clarified the nature of dominant-plane compensation, showing that it can be considered
as a small-translation approximation rather than a planar-scene approximation.



--R

"Recursive Structure and Motion from Image Sequences using Shape and Depth Spaces,"
"Optimal Structure from Motion: Local Ambiguities and Global Esti- mates,"
"A data set for quantitative motion analysis,"
Matrix Computations
"Direct Multi-Resolution Estimation of Ego-Motion and Structure from Motion,"
"In Defense of the Eight-Point Algorithm,"
"Projective structure and motion from from image sequences using subspace methods,"
"Direct Methods for Recovering Motion,"
"A Unified Approach to Moving Object Detection in 2D and 3D Scenes,"
"Linear subspace methods for recovering translational direction,"
"Direct recovery of shape from multiple views: A parallax based approach,"
"Sensitivity of the Pose Refinement Problem to Accurate Estimation of Camera Parameters,"
"Removal of translation bias when using subspace methods,"
"Three Algorithms for 2-Image and  2-Image Structure from Motion,"
"Fast and Accurate Self-Calibration"
"A Multi-frame Structure from Motion Algorithm under Perspective Projection,"
"A New Structure from Motion Ambiguity,"
Recovering Heading and Structure for Constant-Direction Motion,"
"A Critique of Structure from Motion Algorithms,"
"Direct Multi-Frame Structure from Motion for Hand-Held Cameras,"
"Structure from Motion using Points, Lines, and Intensities,"
"Fast Algorithms for Projective Multi-Frame Structure from Motion,"
"Computing the Camera Heading from Multiple Frames,"
"Multiframe Structure from Motion in Perspective,"
"A Linear Solution for Multiframe Structure from Motion,"
"Structure from Linear and Planar Motions,"
"Simplifying Motion and Structure Analysis Using Planar Parallax and Image Warping,"
"Extracting Structure from Optical Flow Using the Fast Error Search Technique,"
"Fast Partial Search Solution to the 3D SFM Problem,"
"A factorization based algorithm for multi-image projective structure and motion,"
"Shape and motion from image streams under orthography: A factorization method,"
"Factorization methods for projective structure and motion,"
--TR

--CTR
Lionel Moisan , Brenger Stival, A Probabilistic Criterion to Detect Rigid Point Matches Between Two Images and Estimate the Fundamental Matrix, International Journal of Computer Vision, v.57 n.3, p.201-218, May-June 2004
John Oliensis, Exact Two-Image Structure from Motion, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.12, p.1618-1633, December 2002
P. Anandan , Michal Irani, Factorization with Uncertainty, International Journal of Computer Vision, v.49 n.2-3, p.101-116, September-October 2002
John Oliensis, The least-squares error for structure from infinitesimal motion, International Journal of Computer Vision, v.61 n.3, p.259-299, February/March 2005
Amit K.  Roy Chowdhury , R. Chellappa, Stochastic Approximation and Rate-Distortion Analysis for Robust Structure and Motion Estimation, International Journal of Computer Vision, v.55 n.1, p.27-53, October
structure from motion;linear algorithms;bundle adjustment;shape from X;sturm/triggs factorization;low level vision;dominant plane;factorization;projective geometry;projective multiframe structure from motion