Fast and Accurate Camera Covariance Computation for Large 3D Reconstruction
Abstract
Estimating uncertainty of camera parameters computed in Structure from Motion (SfM) is an important tool for evaluating the quality of the reconstruction and guiding the reconstruction process. Yet, the quality of the estimated parameters of large reconstructions has been rarely evaluated due to the computational challenges. We present a new algorithm which employs the sparsity of the uncertainty propagation and speeds the computation up about ten times w.r.t. previous approaches. Our computation is accurate and does not use any approximations. We can compute uncertainties of thousands of cameras in tens of seconds on a standard PC. We also demonstrate that our approach can be effectively used for reconstructions of any size by applying it to smaller subreconstructions.
Keywords:
uncertainty, covariance propagation, Structure from Motion, 3D reconstruction1 Introduction
Threedimensional reconstruction has a wide range of applications (e.g. virtual reality, robot navigation or selfdriving cars), and therefore is an output of many algorithms such as Structure from Motion (SfM), Simultaneous location and mapping (SLAM) or Multiview Stereo (MVS). Recent work in SfM and SLAM has demonstrated that the geometry of threedimensional scene can be obtained from a large number of images [1],[14],[16]. Efficient nonlinear refinement [2] of camera and point parameters has been developed to produce optimal reconstructions.
The uncertainty of detected points in images can be efficiently propagated in case of SLAM [16],[28] into the uncertainty o threedimensional scene parameters thanks to fixing the first camera pose and scale. In SfM framework, however, we are often allowing for gauge freedom [18], and therefore practical computation of the uncertainty [9] is mostly missing in the state of the art pipelines [23],[30],[32].
In SfM, reconstructions are in general obtained up to an unknown similarity transformation, i.e., a rotation, translation, and scale. The backward uncertainty propagation [13] (the propagation from detected feature points to parameters the of the reconstruction) requires the “inversion” of a Fischer information matrix, which is rank deficient [9],[13] in this case. Naturally, we want to compute the uncertainty of the inner geometry [9] and ignore the infinite uncertainty of the free choice of the similarity transformation. This can be done by the MoorePenrose (MP) inversion of the Fisher information matrix [9],[13],[18]. However, the MP inversion is a computationally challenging process. It has cubic time and quadratic memory complexity in the number of columns of the information matrix, i.e., the number of parameters.
Fast and numerically stable uncertainty propagation has numerous applications [26]. We could use it for selecting the next best view [10] from a large collection of images [1],[14], for detecting wrongly added cameras to existing partial reconstructions, for improving fitting to the control points [21], and for filtering the mostly unconstrained cameras in the reconstruction to speed up the bundle adjustment [2] by reducing the size of the reconstruction. It would also help to improve the accuracy of the iterative closest point (ICP) algorithm [5], by using the precision of the camera poses, and to provide the uncertainty of the points in 3D [27].
2 Contribution
We present the first algorithm for uncertainty propagation from input feature points to camera parameters that works without any approximation of the natural form of the covariance matrix on thousands of cameras. It is about ten times faster than the state of the art algorithms [19],[26]. Our approach builds on top of GaussMarkov estimation with constraints by Rao [29]. The novelty is in a new method for nullspace computation in SfM. We introdice a fast sparse method, which is independent on a chosen parametrization of rotations. Further, we combine the fixation of gauge freedom by nullspace, from Förstner and Wrobel [9] and methods applied in SLAM, i.e., the block matrix inversion [6] and Woodbury matrix identity [12].
Our main contribution is a clear formulation of the nullspace construction, which is based on the similarity transformation between parameters of the reconstruction. Using the nullspace and the normal equation from [9], we correctly apply the block matrix inversion, which has been done only approximately before [26]. This brings an improvement in accuracy as well as in speed. We also demonstrate that our approach can be effectively used for reconstructions of any size by applying it to smaller subreconstructions. We show empirically that our approach is valid and practical.
Our algorithm is faster, more accurate and more stable than any previous method [19],[26],[27]. The output of our work is publicly available as source code which can be used as an external library in nonlinear optimization pipelines, like Ceres Solver [2] and reconstruction pipelines like [23],[30],[32]. The code, datasets, and detailed experiments will be available online https://michalpolic.github.io/usfm.github.io.
3 Related work
The uncertainty propagation is a well known process [9],[13],[18],[26]. Our goal is to propagate the uncertainties of input measurements, i.e. feature points in images, into the parameters of the reconstruction, e.g. poses of cameras and positions of points in 3D, by using the projection function [13]. For the purpose of uncertainty propagation, a nonlinear projection function is in practice often replaced by its first order approximation using its Jacobian matrix [8],[13]. For the propagation using higher order approximations of the projection function, as described in Förstner and Wrobel [9], higher order estimates of uncertainties of feature points are required. Unfortunately, these are difficult to estimate [9, 25] reliably.
In the case of SfM, the uncertainty propagation is called the backward propagation of nonlinear function in overparameterized case [13] because of the projection function, which does not fully constrain the reconstruction parameters [22], i.e., the reconstruction can be shifted, rotated and scaled without any change of image projections.
We are primarily interested in estimating inner geometry , e.g. angles and ratios of distance, and its inner precision [9]. Inner precision is invariant to changes of gauge, i.e. to similarity transformations of the cameras and the scene [18]. A natural choice of the fixation of gauge, which leads to the inner uncertainty of inner geometry, is to fix seven degrees of freedom caused by the invariance of the projection function to the similarity transformation of space [9],[13],[18]. One way to do this is to use the MoorePenrose (MP) inversion [24] of the Fisher information matrix [9].
Recently, several works on speeding up the MP inversion of the information matrix for SfM frameworks have appeared. Lhuillier and Perriollat [19] used the block matrix inversion of the Fisher information matrix. They performed MP inversion of the Schur complement matrix [34] of the block related to point parameters and then projected the results to the space orthogonal to the similarity transformation constraints. This approach allowed working with much larger scenes because the square Schur complement matrix has the dimension equal to the number of camera parameters, which is at least six times the number of cameras, compared to the mere dimension of the square Fisher information matrix, which is just about three times the number of points.
However, it is not clear if the decomposition of Fisher information matrix holds for MP inversion without fulfilling the rank additivity condition [33] and it was shown in [26] that approach [19] is not always accurate enough. Polic et al. [26] evaluated the state of the art solutions against more accurate results computed in high precision arithmetics, i.e. using 100 digits instead of 15 significant digits of double precision. They compared the influence of several fixations of the gauge on the output uncertainties and found that fixing three points that are far from each other together with a clever approximation of the inversion leads to a good approximation of the uncertainties.
The most related work is [29], which contains uncertainty formulation for GaussMarkov model with constraints. We combine this result with our new approach for nullspace computation to fixing gauge freedom.
Finally, let us mention work on fast uncertainty propagation in SLAM. The difference between SfM and SLAM is that in SLAM we know, and fix, the first camera pose and the scale of the scene which makes the information matrix full rank. Thus one can use a fast Cholesky decomposition to invert a Schur complement matrix as well as other techniques for fast covariance computation [16, 17]. Polok, Ila et al. [15],[28] claim addressing uncertainty computation in SfM but actually assume full rank Fisher information matrix and hence do not deal with gauge freedom. In contrary, we solved here the full SfM problem which requires dealings with gauge freedom.
4 Problem formulation
In this section, we describe basic notions in uncertainty propagation in SfM and provide the problem formulation.
The set of parameters of threedimensional scene is composed from cameras and points in 3D. The th camera is a vector , which consist of internal parameters (i.e. focal length and radial distortion ) and external parameters (i.e. rotation and camera center ). Estimated parameters are labelled with the hat .
We consider that the parameters were estimated by a reconstruction pipeline using a vector of observations . Each observation is a 2D point in the image detected up to some uncertainty that is described by its covariance matrix . It characterizes the Gaussian distribution assumed for the detection error and can be computed from the structure tensor [7] of the local neighbourhood of in the image . The vector is a projection of point into the image plane described by camera parameters . All pairs of indices are in the index set that determines which point is seen by which camera
(1)  
(2) 
Next, we define function and vector as a composition of all projection functions and related detection errors
(3) 
This function is used in the nonlinear least squares optimization (Bundle Adjustment [2])
(4) 
which minimises the sum of squared differences between the measured feature points and the projections of the reconstructed 3D points. We assume the as a block diagonal matrix composed of blocks. The optimal estimate , minimising the Mahalanobis norm, is
(5) 
To find the formula for uncertainty propagation, the nonlinear projection functions can be linearized by the first order term of its Taylor expansion
(6)  
(7) 
which leads to the estimated correction of the parameters
(8) 
Partial derivatives of the objective function must vanishing in the optimum
(9) 
which defines the normal equation system
(10)  
(11) 
The normal equation system has seven degrees of freedom and therefore requires to fix seven parameters, called the gauge [18], namely a scale, a translation and a rotation. Any choice of fixing these parameters leads to a valid solution.
The natural choice of covariance, which is unique, has the zero uncertainty in the scale, the translation, and rotation of all cameras and scene points. It can be obtained by the MP inversion of Fisher information matrix or by GaussMarkov Model with constraints [9]. If we assume a constraints , which fix the scene scale, translation and rotation, we can write their derivatives, i.e. the nullspace , as
(12) 
Using Lagrange multipliers , we are minimising the function
(13) 
that has partial derivative with respect equal to zero in the optimum (as in Eqn. 9)
(14) 
This constraints lead to the extended normal equations
(15) 
and allow us to compute the inversion instead of MP inversion
(16) 
5 Solution method
We next describe how to compute the nullspace and decompose the original Eqn. 16 by a block matrix inversion. The proposed method assumes that the Jacobian of the projection function is provided numerically and provides the nullspace independently of the representation of the camera rotation.
5.1 The nullspace of the Jacobian
The scene can be transformed by a similarity transformation
(17) 
depending on seven parameters for translation, rotation, and scale without any change of the projection function . If we assume a difference similarity transformation, we obtain the total derivative
(18) 
Since it needs to hold for any , the matrix
(19) 
is the nullspace of . Next, consider an order of parameters such that 3D point parameters follow the camera parameters
(20) 
The cameras have parameters ordered as and the projection function equals
(21) 
where projects vectors from to by (i) first dividing by the third coordinate, and (ii) then applying image distortion with parameters . Note that function can be chosen quite freely, e.g. adding a tangential distortion or encountering a rolling shutter projection model [3]. Using Eqn. 17, we are getting for
(22)  
(23)  
(24) 
Note that for any parameters , the projection remains unchanged. It can be checked by expanding the equation above. Eqn. 24 is linear in and . The differences of and are as follows
(25)  
(26) 
The Jacobian and the nullspace can be written as
(27) 
where is the observation, i.e. the pair . The columns of are related to transformation parameters . The rows are related to parameters . The derivatives of differences of scene parameters and with respect to the transformation parameters are exactly the blocks of the nullspace
(28) 
where is the skew symmetric matrix such that for all .
Eqn. 24 is not linear in rotation . To deal with any rotation representation, we can compute the values of for all using Eqn. 18. The columns, which contain blocks , are orthogonal to the rest of the nullspace and to the Jacobian .
The system of equations can be rewritten as
(29) 
where is composed as a blockdiagonal matrix from the red submatrices (see Fig. 1) of . The matrix is composed from red submatrices as
(30) 
The matrix is composed of the green submatrices (see Fig. 1) of multiplied by the minus green submatrices of . The solution to this system is
(31) 
where is computed by a sparse multiplication, see Fig. 1. The inversion of is the inversion of a sparse matrix with blocks on the diagonal.
5.2 Uncertainty propagation to camera parameters
The propagation of uncertainty is based on Eqn. 16. The inversion of extended Fisher information matrix is first conditioned for better numerical accuracy as follows
(32)  
(33)  
(34) 
by diagonal matrices , which condition the columns of matrices , . Secondly, we permute the columns of to have point parameters followed by the camera parameters
(35) 
where is an appropriate permutation matrix. The matrix is a full rank matrix which can be decomposed and inverted using a block matrix inversion
(36) 
where is the symmetric Schur complement matrix of point parameters block
(37) 
Matrix is a sparse symmetric block diagonal matrix with blocks on the diagonal, see Fig. 2. The covariances for camera parameters are computed using the inversion of with the size for our model of cameras (i.e., )
(38) 
where is the left top submatrix of and is the corresponding subblock of scale matrix .
6 Uncertainty for subreconstructions
The algorithm based on GaussMarkov estimate with constraints, which is described in Section 5, works in principle properly for thousands of cameras. However, largescale reconstructions with thousands cameras would require a large space, e.g. 131GB for Rome dataset [20], to store the matrix for our camera model , and its inversion might be inaccurate due to rounding errors.
Fortunately, it is possible to evaluate the uncertainty of a camera from only a partial subreconstruction comprising cameras and points in the vicinity of . Using subreconstructions, we can approximate the uncertainty computed from a complete reconstruction. The error of our approximation decreases with increasing size of a subreconstruction. If we add a camera to a reconstruction, we add at least four observations which influence the Fisher information matrix as
(39) 
where the matrix is the Fisher information matrix of the added observations. We can propagate this update using equations in Section 5 to the Schur complement matrix
(40) 
which has full rank. Using Woodbury matrix identity
(41) 
we can see that the positive definite covariance matrices are subtracted after adding some observations, i.e. the uncertainty decreases.
We show empirically that the error decreases with increasing the size of the reconstruction (see Fig. 3). We have found that for 100–150 neighbouring cameras, the error is usually small enough to be used in practice. Each evaluation of the subreconstruction produces an upper bound on the uncertainty for cameras involved in the subreconstruction. The accuracy of the upper bound depends on a particular decomposition of the complete reconstruction into subreconstructions. To get reliable results, it is useful to decompose the reconstruction several times and choose the covariance matrix with the smallest trace.
The theoretical proof of the quality of this approximation and selection of the optimal decomposition is an open question for future research.
7 Experimental evaluation
We use synthetic as well as real datasets (Table 1) to test and compare the algorithms (Table 2) with respect the accuracy (Fig. 3) and speed (Fig. 5). The evaluations on subreconstructions are shown in Figs. 5, 5(a), 5(b). All experiments were performed on a single computer with one 2.6GHz Intel Core i76700HQ with 32GB RAM running a 64bit Windows 10 operating system.
#  Dataset  

1  Cube  6  15  60 
2  Toy  10  60  200 
3  Flat  30  100  1033 
4  Daliborka  64  200  5205 
5  Marianska  118  80 873  248 511 
6  Dolnoslaskie  360  529 829  226 0026 
7  Tower of London  530  65 768  508 579 
8  Notre Dame  715  127 431  748 003 
9  Seychelles  1400  407 193  2 098 201 
Compared algorithms
are listed in Table 2.
The standard way of computing the covariance matrix is by using the MP inversion of the information matrix using the Singular Value Decomposition (SVD) with the last seven singular values set to zeros and inverting the rest of them as in [26]. There are many implementations of this procedure that differ in numerical stability and speed. We compared three of them. Alg. 1 uses high precision number representation in Maple (runs 22 hours on Daliborka dataset), Alg. 2 denotes the implementation in Ceres [2], which uses Eigen library [11] internally (runs 25.9 minutes on Daliborka dataset) and Alg. 3 is our Matlab implementation, which internally calls LAPACK library [4] (runs 0.45 seconds on Daliborka dataset). Further, we compared Lhuilier [19] and Polic [26] approaches, which approximate the uncertainty propagation, with our algorithm denoted as Nullspace bounding uncertainty propagation (NBUP).
#  Algorithm 

1.  MP inversion of using Maple (Kanatani [18]) (Ground Truth) 
2.  MP inversion of using Ceres (Kanatani [18]) 
3.  MP inversion of using Matlab (Kanatani [18]) 
4.  MP inversion of Schur complement matrix with correction term (Lhuillier [19]) 
5.  TE inversion of Schur complement matrix with three points fixed (Polic [27]) 
6.  Nullspace bounding uncertainty propagation (NBUP) 
The accuracy
of all algorithms is compared against the Ground Truth (GT) in Fig. 3. The evaluation is performed on the first four datasets which have reasonably small number of 3D points. The computation of GT for the fourth dataset took about 22 hours and larger datasets were uncomputable because of time and memory requirements. We decomposed information matrix using SVD, set exactly the last seven singular values to zero and inverted the rest of them. We also used 100 significant digits instead of 15 digits used by a double number representation. The GT computation follows approach from [26].
The covariance matrices for our camera model (comprising rotation, camera center, focal length and radial distortion) contain a large range of values. Some parameters, e.g. rotations represented by the Euler vector, are in units while other parameters, as the focal length, are in thousands of units. Moreover, the rotation is in all tested examples better constrained than the focal length. This fact leads to approximately mean absolute value in rotation part of the covariance matrix and approximately mean value for the focal length variance. Standard deviations for datasets 14 and are about for rotations and for focal lengths. To obtain comparable standard deviations for different parameters, we can divide the mean values of rotations by and focal length by . We used the same approach for the comparison of the measured errors
(41) 
The error shows the differences between GT covariance matrices and the computed ones . The matrix
(42) 
has dimension and normalises the error to percentages of the absolute magnitude of the original units. Symbol stands for elementwise division of matrices (i.e. equals for ).
Fig. 3 shows the comparison of the mean of the errors for all cameras in the datasets. We see that our new method, NBUP, delivers the most accurate results on all datasets.
Speed
of the algorithms is shown in Fig. 5. Note that the MP inversion (i.e. Alg. 13) cannot be evaluated on medium and larger datasets 59 because of memory requirements for storing dense matrix . We see that our new method NBUP is faster than all other methods. Considerable speedup is obtained on datasets 79 where our NBUP method is about 8 times faster.


Uncertainty approximation
on subreconstructions was tested on datasets 59. We decomposed reconstructions several times using a different number of cameras inside smaller subreconstructions, and measured relative and absolute errors of approximated covariances for cameras parameters. Fig. 6 shows the decrease of error for larger subreconstructions. There were 25 subreconstructions for each with the set of neighbouring cameras randomly selected using the view graph. Note that Fig. 5(a) shows the mean of relative errors given by Eqn. 41. Fig. 5(b) shows that the absolute covariance error decreases significantly with increasing the number of cameras in a subreconstruction.
Fig. 5 shows the error of the simplest approximation of covariances used in practice. For every camera, one hundred of its neighbours using viewgraph were used to get a subreconstruction for evaluating the uncertainties. It produces upper bound estimates for the covariances for each camera from which we selected the smallest one, i.e. the covariance matrix with the smallest trace, and evaluate the mean of the relative error .
8 Conclusions
Current methods for evaluating of the uncertainty [19],[26] in SfM rely 1) either on imposing the gauge constraints by using a few parameters as observations, which does not lead to the natural form of the covariance matrix, or 2) on the MoorePenrose inversion [2], which cannot be used in case of medium and largescale datasets because of cubic time and quadratic memory complexity.
We proposed a new method for the nullspace computation in SfM and combined it with Gauss Markov estimate with constraints [29] to obtain a fullrank matrix [9] allowing robust inversion. This allowed us to use efficient methods from SLAM such as block matrix inversion or Woodbury matrix identity. Our approach is the first one which allows a computation of natural form of the covariance matrix on scenes with more than thousand of cameras, e.g. 1400 cameras, with affordable computation time, e.g. 60 seconds, on a standard PC. Further, we show that using subreconstruction of roughly 100300 cameras provides reliable estimates of the uncertainties for arbitrarily large scenes.
9 Acknowledgement
This work was supported by the European Regional Development Fund under the project IMPACT (reg. no. CZ.02.1.01/0.0/0.0/15_003/0000468), EUH2020 project LADIO no. 731970, and by Grant Agency of the CTU in Prague projects SGS16/230/OHK3/3T/13, SGS18/104/OHK3/1T/37.
Footnotes
 The variable is a function of and
References
 Agarwal, S., Furukawa, Y., Snavely, N., Simon, I., Curless, B., Seitz, S.M., Szeliski, R.: Building rome in a day. Communications of the ACM 54(10), 105–112 (2011)
 Agarwal, S., Mierle, K., Others: Ceres solver. http://ceressolver.org
 Albl, C., Kukelova, Z., Pajdla, T.: R6prolling shutter absolute camera pose. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 2292–2300 (2015)
 Anderson, E., Bai, Z., Dongarra, J., Greenbaum, A., McKenney, A., Du Croz, J., Hammarling, S., Demmel, J., Bischof, C., Sorensen, D.: Lapack: A portable linear algebra library for highperformance computers. In: Proceedings of the 1990 ACM/IEEE conference on Supercomputing. pp. 2–11. IEEE Computer Society Press (1990)
 Besl, P.J., McKay, N.D.: Method for registration of 3d shapes. In: Sensor Fusion IV: Control Paradigms and Data Structures. vol. 1611, pp. 586–607. International Society for Optics and Photonics (1992)
 Eves, H.W.: Elementary matrix theory. Courier Corporation (1966)
 Förstner, W.: Image Matching, vol. II, chap. 16, pp. 289–379. Addison Wesley (1993)
 Förstner, W.: Uncertainty and projective geometry. In: Handbook of Geometric Computing, pp. 493–534. Springer (2005)
 Förstner, W., Wrobel, B.P.: Photogrammetric Computer Vision. Springer (2016)
 Frahm, J.M., FiteGeorgel, P., Gallup, D., Johnson, T., Raguram, R., Wu, C., Jen, Y.H., Dunn, E., Clipp, B., Lazebnik, S., et al.: Building rome on a cloudless day. In: European Conference on Computer Vision. pp. 368–381. Springer (2010)
 Guennebaud, G., Jacob, B., et al.: Eigen v3.3. http://eigen.tuxfamily.org (2010)
 Hager, W.W.: Updating the inverse of a matrix. SIAM review 31(2), 221–239 (1989)
 Hartley, R., Zisserman, A.: Multiple view geometry in computer vision. Cambridge university press (2003)
 Heinly, J., Schönberger, J.L., Dunn, E., Frahm, J.M.: Reconstructing the World* in Six Days *(As Captured by the Yahoo 100 Million Image Dataset). In: Computer Vision and Pattern Recognition (CVPR) (2015)
 Ila, V., Polok, L., Solony, M., Istenic, K.: Fast incremental bundle adjustment with covariance recovery. In: International Conference on 3D Vision (3DV) (Oct 2017)
 Ila, V., Polok, L., Solony, M., Svoboda, P.: Slam++a highly efficient and temporally scalable incremental slam framework. The International Journal of Robotics Research 36(2), 210–230 (2017)
 Kaess, M., Dellaert, F.: Covariance recovery from a square root information matrix for data association. Robotics and autonomous systems 57(12), 1198–1210 (2009)
 Kanatani, K.i., Morris, D.D.: Gauges and gauge transformations for uncertainty description of geometric structure with indeterminacy. IEEE Transactions on Information Theory 47(5), 2017–2028 (2001)
 Lhuillier, M., Perriollat, M.: Uncertainty ellipsoids calculations for complex 3d reconstructions. In: Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006. pp. 3062–3069. IEEE (2006)
 Li, Y., Snavely, N., Huttenlocher, D.P.: Location recognition using prioritized feature matching. In: European conference on computer vision. pp. 791–804. Springer (2010)
 Maurer, M., Rumpler, M., Wendel, A., Hoppe, C., Irschara, A., Bischof, H.: Georeferenced 3d reconstruction: Fusing public geographic data and aerial imagery. In: Robotics and Automation (ICRA), 2012 IEEE International Conference on. pp. 3557–3558. IEEE (2012)
 Morris, D.D.: Gauge freedoms and uncertainty modeling for 3 D computer vision. Ph.D. thesis, Citeseer (2001)
 Moulon, P., Monasse, P., Marlet, R., Others: Openmvg. an open multiple view geometry library. https://github.com/openMVG/openMVG
 Nashed, M.Z.: Generalized Inverses and Applications: Proceedings of an Advanced Seminar Sponsored by the Mathematics Research Center, the University of WisconsinâMadison, October 810, 1973. No. 32, Elsevier (2014)
 Polic, M.: 3d scene analysis. http://cmp.felk.cvut.cz/~policmic
 Polic, M., Pajdla, T.: Camera uncertainty computation in large 3d reconstruction. In: International Conference on 3D Vision (2017)
 Polic, M., Pajdla, T.: Uncertainty computation in large 3d reconstruction. In: Scandinavian Conference on Image Analysis. pp. 110–121. Springer (2017)
 Polok, L., Ila, V., Smrz, P.: 3d reconstruction quality analysis and its acceleration on gpu clusters. In: 2016 24th European Signal Processing Conference (EUSIPCO). pp. 1108–1112 (Aug 2016). https://doi.org/10.1109/EUSIPCO.2016.7760420
 Rao, C.R., Rao, C.R., Statistiker, M., Rao, C.R., Rao, C.R.: Linear statistical inference and its applications, vol. 2. Wiley New York (1973)
 Schönberger, J.L., Frahm, J.M.: Structurefrommotion revisited. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016)
 Snavely, N., Seitz, S.M., Szeliski, R.: Photo tourism: exploring photo collections in 3d. In: ACM transactions on graphics (TOG). vol. 25, pp. 835–846. ACM (2006)
 Sweeney, C.: Theia multiview geometry library: Tutorial & reference. http://theiasfm.org
 Tian, Y.: The moorepenrose inverses of m n block matrices and their applications. Linear algebra and its applications 283(1), 35–60 (1998)
 Zhang, F.: The Schur Complement and Its Applications. Springer US (2005)