2 Early works

This section covers some of the relatively early works on SfM, which introduced the basic problem formalism and the factorization-based methods. We first consider the seminal work by Longuet-Higgins (Reference Longuet-Higgins1981), who introduced the first linear method based on point correspondences, later named the eight-point algorithm, to solve the SfM problem for a pair of cameras. Specifically, for a pair of cameras, Longuet-Higgins (Reference Longuet-Higgins1981) aimed to estimate the relative camera motion, that is, relative rotation and translation, and the 3D coordinates of the scene points captured by these cameras. Here we represent the $i$ th camera by $(R_{i},\mathbf{t}_{i})$ , where $R_{i}\in \text{SO}(3)$ is the camera’s orientation matrix and $\mathbf{t}_{i}\in \mathbb{R}^{3}$ is its position. Considering that we can always fix the (intrinsic) coordinate system of one of the two cameras to be the global coordinate system (i.e. we can set $R_{i}=I$ and $\mathbf{t}_{i}=\mathbf{0}$ , since an absolute coordinate system cannot be determined from point correspondences), solving for the relative motion for a camera pair, and for the scene points based on the computed motion, then corresponds to actually solving the SfM problem for the pair. The problem setupFootnote ¹ involves a simple pinhole camera model, in which a scene point $\mathbf{P}\in \mathbb{R}^{3}$ is represented in the $i$ th image plane by $\mathbf{p}_{i}\in \mathbb{R}^{3}$ , as shown in Figure 2.1; see also Özyeşil et al. (Reference Özyeşil, Singer and Basri2015, Section 4) and Hartley and Zisserman (Reference Hartley and Zisserman2000), for example. Here, we obtain $\mathbf{p}_{i}$ by first representing $\mathbf{P}$ , in the $i$ th camera’s coordinate system, as $\mathbf{P}_{i}=R_{i}^{T}(\mathbf{P}-\mathbf{t}_{i})=(\mathbf{P}_{i}^{x},\mathbf{P}_{i}^{y},\mathbf{P}_{i}^{z})^{T}$ , and then projecting it onto the $i$ th image plane via $\mathbf{p}_{i}=(f_{i}/P_{i}^{z})\mathbf{P}_{i}$ , where $R_{i}\in \text{SO}(3)$ denotes the orientation matrix, $\mathbf{t}_{i}\in \mathbb{R}^{3}$ denotes the location of the focal point, and $f_{i}\in \mathbb{R}^{+}$ is the focal length of the $i$ th camera. For a pair of cameras $i$ and $j$ as in Figure 2.1, we can restate the coplanarity of the points $\mathbf{P},\mathbf{t}_{i}$ and $\mathbf{t}_{j}$ , that is, $[(\mathbf{P}-\mathbf{t}_{i})\times (\mathbf{P}-\mathbf{t}_{j})]^{T}(\mathbf{t}_{i}-\mathbf{t}_{j})=0$ , in terms of the observable corresponding point measurements $\mathbf{p}_{i},\mathbf{p}_{j}$ (see Figure 3.1 below for a depiction of a set of corresponding points between a pair of images) and the camera parameters as

where $[\mathbf{t}]_{\times }$ is the skew-symmetric matrix corresponding to the cross-product with $\mathbf{t}$ , and $E_{ij}=[R_{i}^{T}(\mathbf{t}_{j}-\mathbf{t}_{i})]_{\times }R_{i}^{T}R_{j}$ is the essential matrix for the cameras $i$ and $j$ . Also, let $T_{ij}:=[R_{i}^{T}(\mathbf{t}_{j}-\mathbf{t}_{i})]_{\times }$ and $R_{ij}:=R_{i}^{T}R_{j}$ , yielding $E_{ij}=T_{ij}R_{ij}$ . Longuet-Higgins (Reference Longuet-Higgins1981) made the crucial observation that, among the various equivalent restatements of the coplanarity of $\mathbf{P},\mathbf{t}_{i}$ and $\mathbf{t}_{j}$ , the useful property of (2.1), known as the epipolar constraint, is that it provides a basis for the estimation of the structured essential matrices in terms of the observable corresponding points. In other words, by fixing the undetermined scale for the entries of $E_{ij}$ (e.g. the Frobenius norm $\Vert E_{ij}\Vert _{F}=1$ , or as in Longuet-Higgins (Reference Longuet-Higgins1981) $\Vert \mathbf{t}_{j}-\mathbf{t}_{i}\Vert =1$ , which is equivalent to $\Vert E_{ij}\Vert _{F}=\sqrt{2}$ ), we can solve for $E_{ij}\in \mathbb{R}^{3\times 3}$ from eight linearly independent epipolar constraints (2.1) corresponding to eight 3D points.Footnote ² Longuet-Higgins (Reference Longuet-Higgins1981) did not provide any algorithm for this purpose, nor did he consider the effects of uncertainties (i.e. noise) in the corresponding point measurements; however, the usual approach in the literature is simply to minimize the sum of squared errors in the epipolar constraints subject to the scale constraint, which is equivalent to finding the singular vector corresponding to the smallest singular value of the resulting data matrix. After estimating $E_{ij}$ (up to a sign to be determined later), Longuet-Higgins (Reference Longuet-Higgins1981) solved for $T_{ij}$ (up to another sign, again to be determined later) by eliminating $R_{ij}$ using $E_{ij}E_{ij}^{T}=T_{ij}T_{ij}^{T}$ . He then solved for $R_{ij}$ , based on its orthogonality, and for the scene points using simple algebraic equations (see Longuet-Higgins Reference Longuet-Higgins1981 for details). Lastly, the undetermined signs of $E_{ij}$ and $T_{ij}$ are fixed by requiring the scene points to lie in front of both cameras. If the scene points are behind both cameras, the sign of $T_{ij}$ is altered, and the relative rotation $R_{ij}$ and the scene points are recalculated. If the scene points fall behind one of the two cameras, the sign of $E_{ij}$ is altered, and $T_{ij},R_{ij}$ and the scene points are recalculated. Further, Longuet-Higgins (Reference Longuet-Higgins1981) briefly discussed some of the degenerate cases of eight 3D points for which the eight-point estimation cannot be used, for example the cases of all points corresponding to the vertices of a cube, at least seven points lying on a plane, six of the eight points corresponding to the vertices of a regular hexagon, or at least four points lying on a line. We note that although Longuet-Higgins (Reference Longuet-Higgins1981) argued for the accuracy of the eight-point algorithm, it was in fact observed to be sensitive to noise, which resulted in the construction of alternative methods, such as a normalized version of the eight-point algorithm (Hartley Reference Hartley1997), for relative pairwise motion estimation. For further reading on alternative methods to obtain epipolar geometries between images, see for example Zhang (Reference Zhang1998) and Hartley and Zisserman (Reference Hartley and Zisserman2000).

Figure 2.1.

The pinhole camera model.

Another early work in the literature is that of Tomasi and Kanade (Reference Tomasi and Kanade1992), who introduced the factorization method for SfM. To model the SfM problem for objects that are relatively distant compared to their sizes, Tomasi and Kanade (Reference Tomasi and Kanade1992) assumed an orthographic camera model, in which the 3D points are measured via parallel projections onto the image plane, consequently ignoring the camera translation along the optical axis. Let the location of the $j$ th 3D point on the $j$ th camera plane be given by $\mathbf{q}_{ij}=(x_{ij},y_{ij})\in \mathbb{R}^{2}$ . In the noiseless case, these orthographic measurements of $m$ 3D points by $n$ cameras are represented in terms of a measurement matrix $W\in \mathbb{R}^{2n\times m}$ satisfying

The crucial observation made by Tomasi and Kanade (Reference Tomasi and Kanade1992) is that if the origin of the global reference system is chosen to be at the centre of the 3D structure points (implying that each row of $W$ sums to zero), then $W$ can be decomposed into

where $\mathbf{u}_{k}$ and $\mathbf{v}_{k}$ represent the orientation of the $k$ th camera (since, in the orthographic model, the third direction given by $\mathbf{u}_{k}\times \mathbf{v}_{k}$ is parallel to the viewing direction of the cameras) and $\mathbf{P}_{k}$ is the $k$ th 3D point. As a result, $W$ has rank $3$ . In the noisy case, which may result in a higher rank for $W$ , Tomasi and Kanade (Reference Tomasi and Kanade1992) considered a valid measurement matrix to be given by the best rank- $3$ approximation to the noisy $W$ in the Frobenius norm sense, which is given by the singular value decomposition (SVD). Note that replacing $R$ and $S$ with $RQ$ and $Q^{-1}S$ , for an arbitrary invertible matrix $Q\in \mathbb{R}^{3\times 3}$ , results in the same measurement matrix $W$ . As a result, in order to compute $R$ and $S$ from $W$ , Tomasi and Kanade (Reference Tomasi and Kanade1992) proposed first computing a decomposition of $W$ using SVD, and then solving for $Q$ using the orthonormality of the $n$ pairs $(\mathbf{u}_{k},\mathbf{v}_{k})$ . That is, for $W=U_{W}\unicode[STIX]{x1D6F4}_{W}V_{W}^{T}$ representing the SVD of $W$ , they first let $\hat{R}=U_{W}\unicode[STIX]{x1D6F4}^{1/2}$ and ${\hat{S}}=\unicode[STIX]{x1D6F4}^{1/2}V_{W}^{T}$ , and then computed the solution $R=\hat{R}Q$ and $S=Q^{-1}{\hat{S}}$ , where $Q$ is given to be a solution to the non-linear system of $3n$ equations

where $\hat{\mathbf{u}}_{k}$ and $\hat{\mathbf{v}}_{k}$ represent the $k$ th and $(k+n)$ th rows of $\hat{R}$ (see (2.3)). They also extended the factorization method to the case of occlusions (i.e. to the case when some scene points are not visible by all of the cameras) and provided empirical evidence demonstrating the stability of their factorization approach on various experimental scenarios.

Sturm and Triggs introduced an extension of the factorization methodology for the multiview SfM with perspective cameras. They considered the basic image projection equations

where, for the $j$ th 3D point represented (by abuse of notation) in the homogeneous coordinates by $\mathbf{P}_{j}\in \mathbb{R}^{4}$ and $C_{i}$ denoting the $i$ th camera matrix, $\mathbf{p}_{ij}\in \mathbb{R}^{3}$ is the projection of the $j$ th point onto the $i$ th camera plane (in the homogeneous coordinates) and $\unicode[STIX]{x1D706}_{ij}$ are the undetermined scales, called projective depths. Similarly to Tomasi and Kanade (Reference Tomasi and Kanade1992), Sturm and Triggs (Reference Sturm and Triggs1996) collected these projective measurements into a rank- $4$ measurement matrix $Y\in \mathbb{R}^{3n\times m}$ given by

The crucial point here is that we do not have direct access to the projective depths $\unicode[STIX]{x1D706}_{ij}$ , and if $\unicode[STIX]{x1D706}_{ij}$ were to be accurately estimated, which constitutes the main contribution of Sturm and Triggs (Reference Sturm and Triggs1996), a factorization technique similar to that of Tomasi and Kanade (Reference Tomasi and Kanade1992) could easily be applied. In order to achieve this goal, Sturm and Triggs (Reference Sturm and Triggs1996) considered the linear system of equations for pairs of cameras $(i,j)$ given by

where the fundamental matrix $F_{ij}\in \mathbb{R}^{3\times 3}$ satisfies $E_{ij}=K_{i}^{T}F_{ij}K_{j}$ , for $K_{i}$ denoting the $i$ th camera calibration matrix, $E_{ij}$ denotes the essential matrix for the pair $(i,j)$ (see (2.1)), and $\mathbf{e}_{ij}$ is the epipole on the $i$ th image plane (see Figure 2.1). Considering the homogeneous linear equations (2.7) for the minimal number of $n-1$ camera pairs, Sturm and Triggs (Reference Sturm and Triggs1996) estimated the projective depths up to global scale. After substitution of the recovered $\unicode[STIX]{x1D706}_{ij}$ in (2.6), they proceeded similarly to Tomasi and Kanade (Reference Tomasi and Kanade1992) and factorized $Y$ into its camera and structure components using SVD. However, this time there is no need to compute an undetermined multiplier $Q$ since the camera and structure parts of $Y$ do not have particular forms. Lastly, Sturm and Triggs (Reference Sturm and Triggs1996) concluded with empirical results to evaluate the performance of their algorithm. A closely related approach, with an additional section for the extension of the algorithm in Sturm and Triggs (Reference Sturm and Triggs1996) to the case of lines, instead of points, as feature measurements, is given by Triggs (Reference Triggs1996).

An excellent account of factorization-based methods, for various camera models and alternative feature measurements, can also be found in the survey by Kanade and Morris (Reference Kanade and Morris1998).

3 Camera location estimation

In this section we discuss the existing literature focusing primarily on the camera location estimation part of the SfM problem. We consider the camera location estimation methods based on corresponding point estimates between pairs of images. In the majority of the methods we discuss, the corresponding point estimates are used to make relative motion measurements between pairs of images, which can be decomposed into pairwise rotational and translational measurements. The generic approach for these methods is to estimate the camera orientations separately based on the pairwise rotational measurements, and to use these camera orientation estimates together with the pairwise translational measurements in order to solve for the camera locations. On the other hand, we also discuss some of the methods in the literature that, to some degree, deviate from the generic recipe mentioned above, and aim to jointly estimate camera orientations and locations, use local structure estimates for location estimation, or consider triplets of cameras.

In order to illustrate the generic approach mentioned above, consider the simple pinhole camera model introduced in Section 2 and depicted in Figure 2.1. As discussed in Section 2, the essential matrix estimates, computed using the epipolar constraints (2.1) in the presence of sufficiently many corresponding scene points (see Figure 3.1 for a real example of pairwise point correspondences) and knowledge of the intrinsic camera parameters, can be uniquely factorized into relative rotational and translational parts, that is, into estimates of $R_{i}^{T}R_{j}$ and $[R_{i}^{T}(\mathbf{t}_{j}-\mathbf{t}_{i})]_{\times }$ . In general, the relative rotational parts are used for the estimation of camera orientations $R_{i}$ ; see Tron et al. (Reference Tron, Zhou and Daniilidis2016) for a detailed account of the existing methods for camera orientation estimation. In order to estimate the camera locations, the estimated orientations can be used directly with the translational parts to obtain estimates of the pairwise directions $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}=(\mathbf{t}_{i}-\mathbf{t}_{j})/\Vert \mathbf{t}_{i}-\mathbf{t}_{j}\Vert$ . An alternative approach is to go back to the epipolar constraints (2.1) to rewrite them as a set of linear homogeneous equations in the camera locations $\mathbf{t}_{i}$ given by

where $\mathbf{y}_{ij}^{k}$ are functions of $R_{i},R_{j}$ and the $k$ th corresponding points $\mathbf{p}_{i}^{k},\mathbf{p}_{j}^{k}$ , and $m_{ij}$ denotes the number of corresponding points for the cameras $i$ and $j$ . These equations can then be used directly for camera location estimation, or to obtain pairwise direction estimates.

Figure 3.1.

Two images from the Notre Dame Cathedral set in Snavely et al. (Reference Snavely, Seitz and Szeliski2006), with corresponding feature points, extracted using the scale-invariant feature transform (SIFT: Lowe Reference Lowe2004).

A crucial point to note about the estimates of the pairwise directions $\unicode[STIX]{x1D6FE}_{ij}$ is that they lack any scale information pertaining to the distance $\Vert \mathbf{t}_{i}-\mathbf{t}_{j}\Vert$ between camera pairs. The difficultiesFootnote ³ arising from this homogeneity can be regarded as the most important contributor to the diversity of the camera location estimation methods in the literature.

A relatively direct approach for location estimation from pairwise directions is based on the homogeneous system of equations

which encapsulates the fact that the difference vectors $\mathbf{t}_{i}-\mathbf{t}_{j}$ are parallel to the pairwise directions $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}$ . Brand, Antone and Teller (Reference Brand, Antone and Teller2004) estimated the locations by minimizing the sum of squared errors in the system obtained by replacing $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}$ in (3.2) with their estimatesFootnote ⁴ and using constraints to fix the global scale and translation, given respectively by $\sum _{i}\Vert \mathbf{t}_{i}\Vert ^{2}=1$ and $\sum _{i}\mathbf{t}_{i}=\mathbf{0}$ , to prevent the trivial solution of $\mathbf{t}_{i}\equiv \mathbf{t}$ for some $\mathbf{t}\in \mathbb{R}^{3}$ , which is summarized as

Note that, considering the $n$ locations stacked into a vector $\mathbf{t}\in \mathbb{R}^{3n}$ , whose $i$ th $3\times 1$ block is equal to $\mathbf{t}_{i}$ , we can rewrite the cost function of (3.3) as $\mathbf{t}^{T}L\mathbf{t}$ , for some appropriate $L\in \mathbb{R}^{3n\times 3n}$ , and the constraints as $\Vert \mathbf{t}\Vert ^{2}=1$ and $V\mathbf{t}=\mathbf{0}$ (for $V=I_{3}\otimes \mathbf{1}_{n}^{T}$ ). Then the solution to the least-squares problem (3.3) is given in terms of the (normalized) eigenvector corresponding to the fourth smallest eigenvalue of $L$ , whose null space includes the three rows of $V$ . Brand et al. (Reference Brand, Antone and Teller2004) also (heuristically) identified various classes of problems, which they referred to as ‘problem pathologies’, resulting in relatively unstable solutions. Among the listed ‘pathologies’, an important one is the ‘underconstrained’ class of pathologies which possess additional degrees of freedom resulting in more than one set of location estimates (not related to each other via a global translation and/or scale transformation) that may be significantly different. As sources of these underconstrained pathologies, Brand et al. (Reference Brand, Antone and Teller2004) identified two (rather extreme) conditions, namely disconnectedness of the measurement graphFootnote ⁵ and the collinearity of all the directions at a single node, and also provide heuristic solutions to resolve the extra degrees of freedom; these two conditions are in fact only specific examples of a rich set of algebraic conditions defining underconstrained problems, as was fully characterized later by Özyeşil et al. (Reference Özyeşil, Singer and Basri2015) and Özyeşil and Singer (Reference Özyeşil and Singer2015). A similar location estimation method based on the system of equations (3.2) is given in Govindu (Reference Govindu2001). In this method, instead of directly minimizing the sum of squared errors in the system (3.2), an iteratively reweighted least-squares solver, where in the $(k+1)$ th iteration each equation in (3.2) is weighted with $1/\Vert \mathbf{t}_{i}^{k}-\mathbf{t}_{j}^{k}\Vert$ , where $\mathbf{t}_{i}^{k}$ denotes the solution of the $k$ th iteration, is used to democratize the contribution of each pair of cameras in the total cost function. Another approach, similar in spirit to these least-squares methods, is introduced in Arie-Nachimson et al. (Reference Arie-Nachimson, Kovalsky, Kemelmacher-Shlizerman, Singer and Basri2012). This time, instead of the system of equations (3.2), a least-squares method is used to minimize the sum of squared errors in the epipolar constraints (3.1).

Figure 3.2.

Estimates for a noisy, synthetic example of the camera location estimation problem for $100$ cameras, demonstrating the clustering estimates of the least-squares (LS) method of Brand et al. (Reference Brand, Antone and Teller2004) versus estimates of the relatively stable semidefinite relaxation (SDR) solver of Özyeşil et al. (Reference Özyeşil, Singer and Basri2015). The segments represent the error incurred by the SDR solution compared to the ground truth.

Even though the least-squares methods of Brand et al. (Reference Brand, Antone and Teller2004), Govindu (Reference Govindu2001) and Arie-Nachimson et al. (Reference Arie-Nachimson, Kovalsky, Kemelmacher-Shlizerman, Singer and Basri2012) are computationally very efficient and, for relatively small datasets, produce acceptable location estimates, the quality of their estimates degrades significantly for large, sparse, unordered and noisy datasets, that is, for most of the real datasets studied in the literature. Typically, they tend to produce spurious location estimates that cluster tightly around a small number of points: see Figure 3.2 for an example. This undesired tendency can be intuitively explained by considering the least-squares optimization problem (3.3), for which a clustering part of the solution results in a very small contribution to the cost function, and a few locations, usually corresponding to nodes in the measurement graph having few edges, are separated from the cluster to satisfy the constraints. In order to resolve this difficulty, several different methods have been proposed in the literature. A relatively straightforward starting point, introduced by Özyeşil et al. (Reference Özyeşil, Singer and Basri2015), is to require the locations to satisfy a maximal proximity condition by incorporating ‘repulsion constraints’ of the form $\Vert \mathbf{t}_{i}-\mathbf{t}_{j}\Vert \geq c$ , for some fixed $c>0$ (taken to be $c=1$ in Özyeşil et al. Reference Özyeşil, Singer and Basri2015), into the least-squares formulation of Brand et al. (Reference Brand, Antone and Teller2004). These non-convex constraints convert the least-squares formulation into the notoriously difficult problem

To approximate this computationally difficult problem, Özyeşil et al. (Reference Özyeşil, Singer and Basri2015) first introduce the rank- $1$ , positive semidefinite matrix variable $T\in \mathbb{R}^{3n\times 3n}$ (where $n$ denotes the number of cameras), the $ij$ th $3\times 3$ block $T_{ij}$ of which is given by $T_{ij}=\mathbf{t}_{i}\mathbf{t}_{j}^{T}$ . Using $T$ , the cost function and the non-convex constraints of (3.4) are given, respectively, by linear functions

and $\operatorname{Tr}({\mathcal{L}}_{ij}(T))\geq 1$ of $T$ , where the operator ${\mathcal{L}}_{ij}(\cdot )$ satisfies

As a result of the linearity in $T$ of the cost function and the constraints, the non-convexity of (3.4) is fully represented by the non-convex rank- $1$ constraint on $T$ (note that positive semidefiniteness is a convex constraint). As a result, Özyeşil et al. (Reference Özyeşil, Singer and Basri2015) remove this final non-convex constraint to obtain the semidefinite relaxation formulation given by

The approximate solution $\mathbf{t}_{i}$ in (3.4) is obtained to be the $i$ th $3\times 1$ block of the leading eigenvector of the solution of the SDR (3.6). The tightnessFootnote ⁶ of this relaxation up to relatively high levels of noise is demonstrated empirically by Özyeşil et al. (Reference Özyeşil, Singer and Basri2015), who also prove a stability of recovery result, which quantifies the amount of distortion in the location estimates to be of the order of the noise present in the pairwise directions, under fairly general assumptions for the noise model. To improve the quality of their location estimates, Özyeşil et al. (Reference Özyeşil, Singer and Basri2015) propose a ‘robust pairwise direction estimation’ method based on the system (3.1), where the $\mathbf{y}_{ij}^{k}$ are viewed as noisy samples from the two-dimensional subspace orthogonal to $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}$ and efficient robust subspace recovery techniques are used to estimate each pairwise direction $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}$ . Further, and perhaps more importantly, Özyeşil et al. (Reference Özyeşil, Singer and Basri2015) provide a complete characterization of well-posed location recovery from pairwise directionsFootnote ⁷ by demonstrating its equivalence to the existing results in the field of parallel rigidity theory. For instance, while camera orientation estimation only requires the connectivity of the measurement graph, connectivity alone is insufficient to make a camera location estimation well-posed; see Figure 3.3 for an example. The complete characterization of well-posed location recovery problems (for arbitrary dimensions $d\geq 2$ ), via (generically) parallel rigid measurement graphs, is summarized in Theorem 3.1.

Figure 3.3.

(a) A noiseless example of the camera location estimation problem for five locations on a connected graph (we assume noiseless direction measurements between location pairs having an edge), which is not parallel rigid (in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ ). Non-uniqueness is demonstrated by two non-congruent location solutions $\{\mathbf{t}_{1},\mathbf{t}_{2},\mathbf{t}_{3},\mathbf{t}_{4},\mathbf{t}_{5}\}$ and $\{\mathbf{t}_{1},\mathbf{t}_{2},\mathbf{t}_{3},\mathbf{t}_{4}^{\prime },\mathbf{t}_{5}^{\prime }\}$ , each of which can be obtained from the other by an independent rescaling of the solution for one of its maximally parallel rigid components. (b) Maximally parallel rigid components of the formation in (a). (c) A parallel rigid formation (in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ ) obtained from the formation in (a) by adding the extra edge $(1,4)$ linking its maximally parallel rigid components.

Özyeşil et al. (Reference Özyeşil, Singer and Basri2015) also provide an efficient randomized algorithm for testing well-posedness. Further, they introduce an efficient alternating direction augmented Lagrangian method (ADM) to solve the SDR (3.6), and a distributed algorithm to apply the SDR method to large sets of images.

An alternative two-step procedure, studied by Wilson and Snavely (Reference Wilson and Snavely2014), is composed of the detection of outliers among the pairwise direction measurements, via a procedure called 1DSfM, followed by a (non-convex) location estimation method; see Figure 3.4 for a simplified illustration of 1DSfM.

Figure 3.4.

A simplified illustration (courtesy of Wilson and Snavely Reference Wilson and Snavely2014) of outlier detection via 1DSfM. (a) Ground truth camera locations with (potentially noisy) pairwise directions available for camera pairs connected with arrows. (b) Available pairwise directions, with an outlier edge shown in red, and two directions, $\hat{i}$ and $p$ , for one-dimensional (1D) projections. (c) Projections of pairwise directions onto $\hat{i}$ and $p$ , resulting in two variants of the minimum feedback arc set problem. (d) Solutions to the 1D problems in (c), depicting the violation (satisfaction) of the 1D constraints by the outlier edge $(3,0)$ in the 1D problem induced by $\hat{i}$ ( $p$ ), since outlier edges may be consistent with only some of the solutions to the 1D subproblems.

The main idea in the outlier detection step is to project all of the pairwise directions onto a one-dimensional (random) subspace and try to solve the resulting ordering problemFootnote ⁸ as consistently as possible, from which potential outliers can be identified as having larger inconsistencies. The detection procedure ‘solves’ several of these randomized one-dimensional problems and identifies the outliers as those directions having an inconsistency score larger than a certain value. After removal of outlier direction measurements, the locations are estimated by an unconstrained maximization, using the Levenberg–Marquardt algorithm, of the non-convex cost function $\sum _{ij}\hat{\boldsymbol{\unicode[STIX]{x1D6FE}}}_{ij}^{T}(\mathbf{t}_{i}-\mathbf{t}_{j})/\Vert \mathbf{t}_{i}-\mathbf{t}_{j}\Vert$ that measures the alignment of the pairwise directions $(\mathbf{t}_{i}-\mathbf{t}_{j})/\Vert \mathbf{t}_{i}-\mathbf{t}_{j}\Vert$ with their ‘cleaned’ estimates $\hat{\boldsymbol{\unicode[STIX]{x1D6FE}}}_{ij}$ . Wilson and Snavely (Reference Wilson and Snavely2014) also provide extensive experimental results for large datasets, demonstrating the accuracy and efficiency of their method.

Instead of the systems of equations (3.1) or (3.2), on which the various least-squares methods are based, an alternative approach aiming to prevent the clustering phenomena is to consider the system of equations

where $\unicode[STIX]{x1D6FC}_{ij}:=\Vert \mathbf{t}_{i}-\mathbf{t}_{j}\Vert$ . Note that if we ignore this defining (non-convex) relation between the $\unicode[STIX]{x1D6FC}_{ij}$ and $\mathbf{t}_{i}$ , (3.8) turns out to be a linear system in $\unicode[STIX]{x1D6FC}_{ij}$ and $\mathbf{t}_{i}$ for given $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}$ . This system forms the basis of the location estimation method used by Tron and Vidal (Reference Tron and Vidal2014), where the cost function is the sum of squares of the errors in (3.8) and the additional constraints $\unicode[STIX]{x1D6FC}_{ij}\geq 1$ are used to prevent the trivially collapsing solution $\mathbf{t}_{i}\equiv \mathbf{0},\unicode[STIX]{x1D6FC}_{ij}\equiv 0$ . That is, for location estimation, Tron and Vidal (Reference Tron and Vidal2014) solve

In fact, the main theme in Tron and Vidal (Reference Tron and Vidal2014) is to provide a distributed algorithm for camera motion estimation, that is, joint estimation of orientations and locations. In this respect, they provide a consensus-based three-step algorithm, and local convergence and partial stability results (for the rotational part). The quadratic location estimation method (3.9) is used to obtain an initial set of locations (i.e. it constitutes the second step of their distributed algorithm), and it is then fused with the rotational part to obtain a joint motion estimation method. In terms of the camera location estimation problem with given pairwise directions, the crucial contribution of Tron and Vidal (Reference Tron and Vidal2014) is using (3.8) as an alternative to (3.1) or (3.2), together with the relaxed repulsion constraints $\unicode[STIX]{x1D6FC}_{ij}\geq 1$ . As reported by Özyeşil and Singer (Reference Özyeşil and Singer2015), for example, this framework provides location estimates with acceptable quality, which do not tend to cluster for large and noisy datasets, while preserving computational efficiency. A closely related optimization method based on (3.8) and the constraints $\unicode[STIX]{x1D6FC}_{ij}\geq 1$ has been introduced by Moulon, Monasse and Marlet (Reference Moulon, Monasse and Marlet2013). The cost function of this method is given by the maximum absolute deviation in the system (3.8).Footnote ⁹ Moulon et al. (Reference Moulon, Monasse and Marlet2013) also introduce robust rotation and pairwise direction estimation methods. We note, however, that such maximal deviation penalty methods (i.e.  $\ell _{\infty }$ -norm minimization methods) are highly sensitive to outliers in the pairwise direction estimates.

4 Structure estimation

This section considers some of the key contributions in the literature mainly concerning the 3D structure estimation part of the SfM problem. These works have a relatively wide spectrum, including incremental and global methods, joint structure and motion estimation methods, and methods specifically targeting very large-scale SfM. Typically, the number of 3D structure points in SfM is much greater than the number of corresponding cameras: for an unordered collection of about $10^{3}$ cameras, for example, it is reasonable to expect of the order of $10^{6}$ structure points. Therefore, the fundamental challenge is the efficient processing of the given data in a consistent way, while maintaining robustness to noise in the data.

The primary algorithms for joint refinement of the 3D structure, the camera motion and possibly also the intrinsic camera parameters (also called the camera calibration parameters) are known as bundle adjustment methods. As indicated in Triggs et al.’s excellent survey of bundle adjustment (see Triggs et al. Reference Triggs, McLauchlan, Hartley and Fitzgibbon2000, Figure 9), their origin can be traced back to the works of Gauss and Legendre on least-squares estimation in the late eighteenth and early nineteenth centuries. However, in this section we will focus on more recent developments, and more generally on structure estimation methods that mainly aim to facilitate the processing of large unordered sets of images.

Triggs et al. (Reference Triggs, McLauchlan, Hartley and Fitzgibbon2000) discussed several different aspects of bundle adjustment, including the basic projection model and problem parametrization, error modelling and the choice of cost function, the main types of numerical optimization algorithms, the sparse problem structure as induced by the network of variables (and how to exploit it to improve reconstruction accuracy), gauge invariance, and quality control to detect outliers and characterize the overall accuracy of the estimated parameters. They presented the projection model and the problem parametrization in a generalized, rather abstract framework to incorporate various scene models (e.g. relatively simple isolated 3D features made up of points, lines or planes, or more complicated models involving dynamics, photometry or complex objects linked by constraints) and camera models (e.g. perspective, affine or orthographic projections, or less frequently encountered models such as rational polynomial cameras). To define a feature prediction error (also sometimes referred to as the total reprojection error), Triggs et al. (Reference Triggs, McLauchlan, Hartley and Fitzgibbon2000) considered a simple scene model comprising static, individual (abstract) 3D features $\mathbf{X}_{j}$ imaged by a set of cameras corresponding to the position and intrinsic calibration parameters $\mathbf{C}_{i}$ . In this framework we are given noisy measurements $\bar{\mathbf{x}}_{ij}$ of the image features $\mathbf{x}_{ij}$ representing the true image of features $\mathbf{X}_{j}$ in image $i$ , that are assumed to satisfy a predictive model given by $\mathbf{x}_{ij}=\mathbf{x}(\mathbf{C}_{i},\mathbf{X}_{j})$ . The feature prediction errors for the measurements $\bar{\mathbf{x}}_{ij}$ are then defined by

To reconstruct the scene, a measure of these prediction errors is minimized, where bundle adjustment corresponds to the refinement part of this optimization, and requires an initial set of estimates for the camera and structure parameters (e.g. computed using a simpler technique). In summary, as noted by Triggs et al. (Reference Triggs, McLauchlan, Hartley and Fitzgibbon2000, p. 302), bundle adjustment ‘is essentially a matter of optimizing a complicated nonlinear cost function (the total prediction error) over a large nonlinear parameter space (the scene and camera parameters)’, the accuracy of which is highly dependent on the initial estimates of camera and structure parameters (explaining, for example, the wide variety of initial motion estimation methods), and on the choice of the parametrization of the parameter space; for a more detailed discussion see Triggs et al. (Reference Triggs, McLauchlan, Hartley and Fitzgibbon2000, Section 2.2). Moreover, Triggs et al. studied the choice of the cost function defining the measure of the errors in (4.1), reaching the main conclusion that robust, statistically based error metrics allowing the presence of outliers in the measurements are appropriate choices. They discussed three main categories of algorithms for the numerical optimization of the chosen cost function in detail, namely the second-order Newton-style methods, first-order methods (demonstrating linear asymptotic convergence) and the sequential methods incorporating a series of observations one by one (as opposed to solving for the whole system globally). For second-order methods, which have fast asymptotic convergence rates but relatively high costs per iteration, they explored several different techniques for improving efficiency, mainly exploiting the sparsity patterns of the Hessian of the cost function. They also discussed various first-order methods, including the simple gradient descent, linear and non-linear Gauss–Seidel methods and Krylov subspace methods. For the third category of sequential algorithms, they considered low-rank update techniques for updating the Hessian and various filtering techniques. Triggs et al. (Reference Triggs, McLauchlan, Hartley and Fitzgibbon2000) also provided methods for handling theoretical and algorithmic difficulties that arise from several gauge symmetries involved in the problem, and (essentially) concluded with various analytical and heuristic techniques for quality control, that are used to maintain algorithmic stability and robustness to outlier measurements, and to characterize the overall accuracy of the estimates.

A relatively recent work on large-scale bundle adjustment (for tens of thousands of images) is that of Agarwal, Snavely, Seitz and Szeliski (Reference Agarwal, Snavely, Seitz and Szeliski2010). The main challenge identified by Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010) is the low sparsity levels (i.e. relatively high density) present in the connectivity graphsFootnote ¹³ of unstructured sets of images, such as community photo collections, as opposed to the higher sparsity of the structured image sets (e.g. street-level datasets such as images captured from a moving truck); for an example see Figure 4.1. The decrease in sparsity drastically increases the computational cost of the sparse factorization methods employed in classical bundle adjustment methods. To address this difficulty, the inexact (truncated) Newton-type algorithm introduced by Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010) makes use of conjugate gradients for the Newton step computation combined with various preconditioners. In more detail, they consider a non-linear least-squaresFootnote ¹⁴ bundle adjustment problem stated (abstractly) as

where $\mathbf{x}\in \mathbb{R}^{p}$ , $\mathbf{F}:\mathbb{R}^{p}\rightarrow \mathbb{R}^{q}$ . Let $J(\mathbf{x})\in \mathbb{R}^{q\times p}$ denote the Jacobian of $\mathbf{F}$ at $\mathbf{x}$ and let

Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010) consider a Levenberg–Marquardt (LM) approach to solving (4.2), which, to compute the update $\mathbf{x}\rightarrow \mathbf{x}+\unicode[STIX]{x1D6E5}\mathbf{x}$ at each iteration (while ensuring convergence by controlling the step size $\Vert \unicode[STIX]{x1D6E5}\mathbf{x}\Vert$ ), solves the regularized least-squares problem

where $D(\mathbf{x})$ is a non-negative diagonal matrix, for example, given by

and $\unicode[STIX]{x1D707}>0$ determines the regularization strength. Note that if the regularized Hessian $H_{\unicode[STIX]{x1D707}}(\mathbf{x}):=J(\mathbf{x})^{T}J(\mathbf{x})+\unicode[STIX]{x1D707}D(\mathbf{x})^{T}D(\mathbf{x})$ is positive definite, we can solve (4.3) via the normal equations

Here, the crucial point is that $H_{\unicode[STIX]{x1D707}}(\mathbf{x})$ usually turns out to have a special block structure induced by the sparsity pattern of the bundle adjustment problem, which simplifies the solution to (4.4). For instance, assuming that the components of $\mathbf{F}(\cdot )$ only depend on a single camera or a single 3D feature, which is usually the case, then we get

where $B$ and $C$ are block diagonal (with blocks of size ${<}10$ , and hence can be efficiently inverted), and $E$ is block sparse. Therefore, by also considering the appropriate block representations $\unicode[STIX]{x1D6E5}\mathbf{x}=[\unicode[STIX]{x1D6E5}\mathbf{y},\unicode[STIX]{x1D6E5}\mathbf{z}]$ and $-g=[v,w]$ , and solving for $\unicode[STIX]{x1D6E5}\mathbf{z}$ by $\unicode[STIX]{x1D6E5}\mathbf{z}=C^{-1}(w-E^{T}\unicode[STIX]{x1D6E5}\mathbf{y})$ , we obtain

Here the Schur complement $S$ of $C$ in $H_{\unicode[STIX]{x1D707}}$ given by $S=B-EC^{-1}E^{T}$ , and known as the reduced camera matrix, has a block sparse structure with a non-zero $ij$ th block if and only if the $i$ th and the $j$ th images have at least one common feature point. As noted by Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010), computationally speaking, the reduction of (4.4) to (4.6), via the Schur complement trick given above, allows us to obtain a solution to (4.4) using classical Cholesky factorization methods (that do not exploit the sparsity pattern in  $S$ ) in $O(n^{2})$ space and $O(n^{3})$ time complexity, where $n$ denotes the number of cameras. These computational costs become prohibitive for large datasets, and hence an alternative is to use decomposition algorithms that exploit the sparsity in $S$ . However, as the datasets grow and the level and simplicity of sparsity in $S$ decreases (as argued by Agarwal et al. Reference Agarwal, Snavely, Seitz and Szeliski2010 for the case of community photo collections), even the construction and storage of $S$ becomes computationally challenging. As a result, instead of trying to solve (4.3) exactly, for example by using the Schur complement trick discussed above, they propose using an inexact Newton method based on an approximate iterative linear system solver, such as conjugate gradient (CG), possessing a termination rule given by

where $0<\unicode[STIX]{x1D702}_{k}\leq \unicode[STIX]{x1D702}_{0}<1$ , $k$ denoting the LM iteration number, is used as a forcing sequence to guarantee convergence. Further, Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010) consider several preconditioners $M$ for which the matrix $M^{-1}H_{\unicode[STIX]{x1D707}}(\mathbf{x})$ of the equivalent linear system $M^{-1}H_{\unicode[STIX]{x1D707}}(\mathbf{x})\unicode[STIX]{x1D6E5}\mathbf{x}=-M^{-1}\mathbf{g}(\mathbf{x})$ has reduced condition number, which therefore improves the CG convergence rate. Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010) conclude their paper by extensively comparing the experimental stability and efficiency of their preconditioned CG method with other preconditioners, and with the classical methods factoring the Schur complement, obtaining significant improvements in accuracy and efficiency, in both time and memory, using their proposed algorithm.

Figure 4.1.

Sparsity patterns of the adjacency matrices of the connectivity graphs for two sets of images: (a) a structured set of photos (taken from a moving truck), and (b) an unstructured collection of community photos corresponding to a search for the term ‘Dubrovnik’ in Flickr. Figure courtesy of Agarwal et al. (Reference Agarwal, Snavely, Seitz and Szeliski2010).

For large, unordered and hence computationally challenging datasets, Crandall, Owens, Snavely and Huttenlocher (Reference Crandall, Owens, Snavely and Huttenlocher2011) have introduced an alternative global approach, whose formulation is based on first computing a coarse initial solution via a hybrid discrete-continuous optimization, using discrete belief propagation (BP) within a Markov random field (MRF) formulation of pairwise camera and camera-point constraints to estimate camera parameters and a continuous LM improvement, and then refining the resulting solution with classical bundle adjustment. This method also employs various alternative sources of data such as noisy geotags and vanishing point estimates extracted from images. More concretely, the method first solves for camera orientations by minimizing a robust cost function of the error in the (classical) pairwise rotational consistency equations (i.e. the equations $R_{ij}=R_{i}^{T}R_{j}$ , for which we have access to noisy measurements of some of the relative rotations $R_{ij}$ ), using a mixture of BP (on a discretization of the 3D sphere) and LM refinement. Having solved for the orientations $R_{i}$ , the joint estimation of the camera locations and structure points is based on the minimization of (a robust cost function of) the total angles of deviation between pairwise translations $\mathbf{t}_{i}-\mathbf{t}_{j}$ and pairwise directions $\boldsymbol{\unicode[STIX]{x1D6FE}}_{ij}$ (see e.g. (3.2)) and between 3D point-to-camera displacements (given by $\mathbf{P}_{k}-\mathbf{t}_{i}$ between the $k$ th 3D point $\mathbf{P}_{k}$ and camera $i$ ) and ‘ray directions’ (given in the noiseless case by $\mathbf{r}_{ik}=R_{i}^{T}K_{i}^{-1}\mathbf{p}_{ik}$ , for the known intrinsic matrix $K_{i}$ of camera $i$ , and the projective measurement $\mathbf{p}_{ik}$ of $\mathbf{P}_{k}$ by camera $i$ ). This minimization is again performed with a hybrid BP and LM method. Even though, in principle, such non-convex methods are highly sensitive to initial points, and hence susceptible to local minima, Crandall et al. (Reference Crandall, Owens, Snavely and Huttenlocher2011) have been able to improve its sensitivity and obtain, as empirically demonstrated, high-quality reconstructions by also making use of external, additional sources of information like geotags and vanishing points. Further, the parallelizable structure of the methods in Crandall et al. (Reference Crandall, Owens, Snavely and Huttenlocher2011) allows the realization of significant improvements in computational efficiency.

A relatively popular main approach to handling challenging datasets is to employ sequential processing of the available images, that is, to incorporate the images one by one, or in small groups, to obtain a solution for the whole set based on a ‘core solution’ computed from a small subset of the available data. A systematic method of this sequential type was introduced by Snavely, Seitz and Szeliski (Reference Snavely, Seitz and Szeliski2008b ), whose main idea was to identify a small subset of the available images, called the skeletal set, by first estimating the accuracy of pairwise reconstructions of overlapping images and then extracting a subset whose reconstruction accuracy approximates that of the total set. The remaining images are then incorporated by position estimation, and a final bundle adjustment can be used to improve accuracy. The computation of the skeletal set is rather complicated,Footnote ¹⁵ but the main idea is to identify a (directed) measure of proximity between cameras, identified as nodes of the image graph, given by the (translational) uncertainty in the pairwise reconstruction (obtained by fixing the location and orientation of one of the cameras in the pair and then repeating it for the other one) as quantified by the trace of the covariance matrix given by the translational part of the inverse of the Hessian of the reduced camera system (4.6) (i.e. the Schur complement $S=B-EC^{-1}E^{T}$ for the camera pair); see Figure 4.2 for a simplified depiction. The skeletal set problem is to identify the smallest subset of cameras, subject to the constraint that the distance (length of the shortest path, as defined by the total uncertainty of the path) between any pair of cameras in the (directed) subgraph induced by this subset is at most $t$ times longer than their distance in the original digraph (induced by the set of all cameras). The factor $t$ is called the stretch factor, and controls the amount of additional uncertainty introduced between camera pairs by removing cameras from the available dataset. Also, the skeletal set is required to satisfy the constraint of being uniquely realizable, which is modelled and controlled via a structure that they call the pair graph. Unfortunately, this problem (even if we ignore the additional realizability constraint), known as computing a minimum $t$ -spanner for the digraph of camera pairs, is NP-complete for general graphs. Hence, Snavely et al. (Reference Snavely, Seitz and Szeliski2008b , Section 4) employ a two-step heuristic approximation to compute the skeletal set. The empirical results provided by Snavely et al. (Reference Snavely, Seitz and Szeliski2008b ) demonstrate a significant increase in overall SfM efficiency, and the accuracy of the reconstructed structure as compared to applying existing SfM techniques to the full set of images.

Figure 4.2.

A toy example of the construction of the skeletal graph in Snavely et al. (Reference Snavely, Seitz and Szeliski2008b ). (a) Nodes of the image graph represent images and its edges represent pairwise reconstructions (in fact the image graph is directed, and each connected pair has two directed edges; only one is shown for simplicity). Each (directed) edge $(I,J)$ is endowed with a covariance matrix $C_{IJ}$ representing the uncertainty in image $J$ relative to $I$ . (b) Relative uncertainty between the nodes $P$ and $Q$ is quantified by the trace of the covariance matrices (ellipses) chained up along the shortest path (arrows) between the nodes. (c) Removing an edge lengthens the shortest path, resulting in larger estimated covariance (larger ellipses). (d) The solid edges are those of the skeletal graph, while the dotted edges have been removed. The black (interior) nodes form the skeletal set, and are to be reconstructed first, while the grey (leaf) nodes are added to it using position estimation. The skeletal graph is constructed by trying to minimize the number of interior nodes, while bounding the increase in relative uncertainty between all pairs of nodes in the original graph in (a). Figure courtesy of Snavely et al. (Reference Snavely, Seitz and Szeliski2008b ).

Another incremental method, similar in spirit to Snavely et al. (Reference Snavely, Seitz and Szeliski2008b ), has been studied by Havlena, Torii and Pajdla (Reference Havlena, Torii and Pajdla2010). Their main idea for improving the efficiency of the overall SfM pipeline is to identify a small subset of input images by computing an approximate minimally connected dominating set, for which each of the remaining images has significant overlaps with at least one of the images in this set, and to employ task prioritization to promote easier image matching. Moreover, instead of computing actual image matchings, which is computationally prohibitive for thousands of cameras if performed on each pair of images, they employ fast image indexing techniques based on large image vocabularies to measure visual overlaps between images. To identify image similarities, they use a ‘bag-of-words approach’ (Havlena et al. Reference Havlena, Torii and Pajdla2010, Section 2.1). The computation of a minimally connected dominating set for the induced (unweighted, undirected) graph, defined as a minimum-sized connected subgraph, of which each vertex of the original graph is either a member or a neighbour, is in general NP-hard. Hence, they use a polynomial-time approximation algorithm from the literature for this task (Havlena et al. Reference Havlena, Torii and Pajdla2010, Algorithm 1). For the 3D structure, they make use of the atomic 3D models from image triplets approach of Havlena, Torii, Knopp and Pajdla (Reference Havlena, Torii, Knopp and Pajdla2009), and grow the structure by using a prioritization queue, which orders the different tasks of either creating a new atomic 3D model or incorporating one image into a given 3D model (Havlena et al. Reference Havlena, Torii and Pajdla2010, Section 3). They conclude with various experimental results, including image sets captured by omnidirectional cameras, demonstrating the efficiency and accuracy of their approach.

A notable work of an incremental nature is the photo tourism system of Snavely et al. (Reference Snavely, Seitz and Szeliski2006), which was also discussed in Section 1. Although the SfM part of their paper is primarily based on a rather classical, generic type of sequential incorporation of images, one by one, using the sparse bundle adjustment algorithm of Lourakis and Argyros (Reference Lourakis and Argyros2009), it provides a very rich set of tools for content exploration, such as 3D scene visualization, object-based photo browsing, geolocation from images, or object annotation; see Figure 4.3 for a simplified illustration of the photo tourism system.

Figure 4.3.

Photo tourism system of Snavely et al. (Reference Snavely, Seitz and Szeliski2006): (a) An unstructured collection of images, e.g. downloaded from the Internet, as input to the system. (b) SfM step composed of camera motion estimation and 3D structure reconstruction. (c) Content exploration, such as 3D scene visualization, object-based photo browsing, geolocation from images, object annotation, etc. Figure courtesy of Snavely et al. (Reference Snavely, Seitz and Szeliski2006).

For sequential SfM, Snavely et al. identify an initial pair of cameras sharing a large set of feature points, which also have a known baseline to improve the accuracy of the reconstructed initial 3D points. The selection criterion for a new camera is to have a maximal number of feature points that correspond to the already reconstructed 3D points. This iterative procedure is repeated until no remaining camera observes any reconstructed 3D point; for additional details, such as detection of outlier tracks, geo-registration, and specific algorithms used, see Snavely et al. (Reference Snavely, Seitz and Szeliski2006, Section 4). Even though the overall computational cost of the proposed SfM pipeline can be fairly high for relatively large datasets,Footnote ¹⁶ it provides an accurate alternative to SfM, and more importantly, the rich set of tools it provides for content exploration strikingly demonstrates the wide range of applicability and usefulness of SfM.

Yet another remarkable SfM approach for processing large (specifically city-scale) unordered datasets has been introduced by Agarwal, Snavely, Simon, Seitz and Szeliski (Reference Agarwal, Snavely, Simon, Seitz and Szeliski2009). The key paradigm they explore is the parallel processing of the various stages in SfM, and hence its main contributions are the identification of challenges and introduction of solutions in the parallelization of the SfM pipeline. In its multi-stage image matching method, where each stage consists of a proposal step (of possible image pairs sharing common features) and a verification step, they employ several different techniques, such as vocabulary tree-based whole-image similarity and query expansion, while simultaneously controlling the computational efficiency of the distributed framework; for details see Agarwal et al. (Reference Agarwal, Snavely, Simon, Seitz and Szeliski2009, Sections 2.2, 2.3). For the SfM part, they use the skeletal set method of Snavely et al. (Reference Snavely, Seitz and Szeliski2008b ) to reduce the computational cost significantly, and they compute the initial camera motion and 3D structure for the skeletal set via the sequential SfM method of Snavely et al. (Reference Snavely, Seitz and Szeliski2006). To improve upon the computed motion and structure, they perform bundle adjustment, where, depending on the size of the problem, they choose between a preconditioned CG method (as in Agarwal et al. Reference Agarwal, Snavely, Seitz and Szeliski2010) for large datasets, and an exact step LM algorithm for smaller datasets. Finally, Agarwal et al. (Reference Agarwal, Snavely, Simon, Seitz and Szeliski2009) demonstrate the accuracy and efficiency of their overall approach via extensive experimental results for various city-scale datasets.

Interestingly, after camera rotations have been estimated, the rest of the SfM pipeline may be reformulated as a bipartite location recovery problem. Namely, camera locations and 3D structure points form the two sets of nodes of the bipartite graph, and directions between cameras and scene points obtained from correspondences between pairs of views form the edges. This observation, coupled with fast algorithms for ShapeFit and LUD, has potentially substantial implications for real-time robotics applications, as discussed in Section 5. Hand, Lee and Voroninski (Reference Hand, Lee and Voroninski2015), in follow-up work to Hand et al. (Reference Hand, Lee and Voroninski2017), have proved a set of corruption-robust exact recovery theorems for ShapeFit in the bipartite setting. As obtained after the correspondence estimation step of the SfM pipeline, Hand et al. (Reference Hand, Lee and Voroninski2015) assume the availability of measurements (partially corrupted with adversarial noise) of camera-to-point directions $(\mathbf{t}_{i}-\mathbf{p}_{j})/\Vert \mathbf{t}_{i}-\mathbf{p}_{j}\Vert$ , where $\mathbf{t}_{i}\in \mathbb{R}^{d}$ is the $i$ th location and $\mathbf{p}_{j}\in \mathbb{R}^{d}$ the $j$ th structure point. Using the convex ShapeFit algorithm introduced in Hand et al. (Reference Hand, Lee and Voroninski2017) for the simultaneous estimation of locations and structure points, they prove a (high-dimensional) exact recovery result for a set of i.i.d. Gaussian locations and structure points with high probability if the underlying bipartite measurement graph is sampled from an Erdős–Rényi ensemble, the problem dimension $d$ is sufficiently large, and at most a constant fraction of direction measurements involving any particular location or structure point are corrupted; for a precise statement of this main result see Hand et al. (Reference Hand, Lee and Voroninski2015, Theorem 1). Although this result is in the high-dimensional regime, it is the first simultaneous exact location and structure recovery result for partially corrupted camera-to-point direction measurements. The proofs in Hand et al. (Reference Hand, Lee and Voroninski2017) crucially rely on the existence of sufficiently many triangles in random graphs, while any bipartite graph cannot contain a triangle. The results in Hand et al. (Reference Hand, Lee and Voroninski2015) establish the necessary properties for exact recovery in high dimensions from corrupted relative-direction observations, similar to Hand et al. (Reference Hand, Lee and Voroninski2017, Theorem 1), by relying instead on the existence of randomly distributed four-cycles in random bipartite graphs. Hand et al. (Reference Hand, Lee and Voroninski2015) also provide empirical evidence on simulated data supporting their theoretical findings, and Goldstein et al. (Reference Goldstein, Hand, Lee, Voroninski and Soatto2016) evaluate bipartite ShapeFit on real datasets, providing evidence of robust recovery in the bipartite setting.

5 Simultaneous localization and mapping (SLAM)

Simultaneous localization and mapping (SLAM) consists in estimating a depth map of the 3D environment, as well as localizing the agent which is moving throughout this environment. SLAM is a fundamental problem in robotics, providing crucial information for autonomous navigation of cars, drones and consumer robots. Pre-built maps of the environment can be used to provide much useful information for navigation, provided that localization within those maps can be performed accurately. The correct assertion that a vehicle or robotic entity is in a previously observed location is called a loop closure. SLAM is effectively a special case of the SfM problem, in that it is specific to robotics and augmented/virtual reality applications, yet retains many of the same difficulties.

A multitude of combinations of sensing techniques have been used for SLAM, including LIDAR, inertial sensors, GPS and vision. LIDAR is a form of active sensing, which consists of measuring laser time of flight to infer depth in real time. While LIDAR to a large extent trivializes the depth estimation problem and also eases the task of loop closure, it is an expensive sensor with low resolution, and suffers performance degradation in sub-optimal weather conditions such as rain and fog. Nonetheless, LIDAR is the centrepiece of many of today’s approaches to autonomous navigation, by allowing pre-mapping of urban environments and accurate loop closures in such maps, which provides information about the static part of the environment in real time. Localization can also be aided by differential GPS, which has very high accuracy but is also expensive. The holy grail of SLAM research is a purely vision-based approach, since vision sensors are cheap and ubiquitous, and provide very rich and high-resolution information about the environment. However, due to the aperture problem and various other ambiguities of the image formation process, using vision alone for SLAM is a difficult task and remains at the forefront of robotics research, in particular due to its high computational complexity. In this section we will cover state-of-the-art vision-based and visual-inertial SLAM approaches as well as some of the particular difficulties faced in robotics applications.

6 Other topics

This section briefly discusses some of the additional topics in the SfM literature, not considered in previous sections. These topics include methods for feature extraction and matching, alternative core measurement types and camera models, methods aiming to handle symmetries and ambiguities in images, and widely utilized sources of data and SfM software.