Fast and Accurate Camera Covariance Computation for Large 3D Reconstruction
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 sub-reconstructions.
Keywords:uncertainty, covariance propagation, Structure from Motion, 3D reconstruction
Three-dimensional reconstruction has a wide range of applications (e.g. virtual reality, robot navigation or self-driving cars), and therefore is an output of many algorithms such as Structure from Motion (SfM), Simultaneous location and mapping (SLAM) or Multi-view Stereo (MVS). Recent work in SfM and SLAM has demonstrated that the geometry of three-dimensional scene can be obtained from a large number of images ,,. Efficient non-linear refinement  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 , into the uncertainty o three-dimensional scene parameters thanks to fixing the first camera pose and scale. In SfM framework, however, we are often allowing for gauge freedom , and therefore practical computation of the uncertainty  is mostly missing in the state of the art pipelines ,,.
In SfM, reconstructions are in general obtained up to an unknown similarity transformation, i.e., a rotation, translation, and scale. The backward uncertainty propagation  (the propagation from detected feature points to parameters the of the reconstruction) requires the “inversion” of a Fischer information matrix, which is rank deficient , in this case. Naturally, we want to compute the uncertainty of the inner geometry  and ignore the infinite uncertainty of the free choice of the similarity transformation. This can be done by the Moore-Penrose (M-P) inversion of the Fisher information matrix ,,. However, the M-P 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 . We could use it for selecting the next best view  from a large collection of images ,, for detecting wrongly added cameras to existing partial reconstructions, for improving fitting to the control points , and for filtering the mostly unconstrained cameras in the reconstruction to speed up the bundle adjustment  by reducing the size of the reconstruction. It would also help to improve the accuracy of the iterative closest point (ICP) algorithm , by using the precision of the camera poses, and to provide the uncertainty of the points in 3D .
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 ,. Our approach builds on top of Gauss-Markov estimation with constraints by Rao . 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  and methods applied in SLAM, i.e., the block matrix inversion  and Woodbury matrix identity .
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 , we correctly apply the block matrix inversion, which has been done only approximately before . 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 sub-reconstructions. We show empirically that our approach is valid and practical.
Our algorithm is faster, more accurate and more stable than any previous method ,,. 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  and reconstruction pipelines like ,,. 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 ,,,. 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 . For the purpose of uncertainty propagation, a non-linear projection function is in practice often replaced by its first order approximation using its Jacobian matrix ,. For the propagation using higher order approximations of the projection function, as described in Förstner and Wrobel , 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 non-linear function in over-parameterized case  because of the projection function, which does not fully constrain the reconstruction parameters , 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 . Inner precision is invariant to changes of gauge, i.e. to similarity transformations of the cameras and the scene . 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 ,,. One way to do this is to use the Moore-Penrose (M-P) inversion  of the Fisher information matrix .
Recently, several works on speeding up the M-P inversion of the information matrix for SfM frameworks have appeared. Lhuillier and Perriollat  used the block matrix inversion of the Fisher information matrix. They performed M-P inversion of the Schur complement matrix  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 M-P inversion without fulfilling the rank additivity condition  and it was shown in  that approach  is not always accurate enough. Polic et al.  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 , which contains uncertainty formulation for Gauss-Markov 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. , 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 three-dimensional 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  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
Next, we define function and vector as a composition of all projection functions and related detection errors
This function is used in the non-linear least squares optimization (Bundle Adjustment )
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
To find the formula for uncertainty propagation, the non-linear projection functions can be linearized by the first order term of its Taylor expansion
which leads to the estimated correction of the parameters
Partial derivatives of the objective function must vanishing in the optimum
which defines the normal equation system
The normal equation system has seven degrees of freedom and therefore requires to fix seven parameters, called the gauge , 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 M-P inversion of Fisher information matrix or by Gauss-Markov Model with constraints . If we assume a constraints , which fix the scene scale, translation and rotation, we can write their derivatives, i.e. the nullspace , as
Using Lagrange multipliers , we are minimising the function
that has partial derivative with respect equal to zero in the optimum (as in Eqn. 9)
This constraints lead to the extended normal equations
and allow us to compute the inversion instead of M-P inversion
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 transformation111The variable is a function of and
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
Since it needs to hold for any , the matrix
is the nullspace of . Next, consider an order of parameters such that 3D point parameters follow the camera parameters
The cameras have parameters ordered as and the projection function equals
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 . Using Eqn. 17, we are getting for
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
The Jacobian and the nullspace can be written as
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
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
where is composed as a block-diagonal matrix from the red submatrices (see Fig. 1) of . The matrix is composed from red submatrices as
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
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
by diagonal matrices , which condition the columns of matrices , . Secondly, we permute the columns of to have point parameters followed by the camera parameters
where is an appropriate permutation matrix. The matrix is a full rank matrix which can be decomposed and inverted using a block matrix inversion
where is the symmetric Schur complement matrix of point parameters block
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., )
where is the left top submatrix of and is the corresponding sub-block of scale matrix .
6 Uncertainty for sub-reconstructions
The algorithm based on Gauss-Markov estimate with constraints, which is described in Section 5, works in principle properly for thousands of cameras. However, large-scale reconstructions with thousands cameras would require a large space, e.g. 131GB for Rome dataset , 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 sub-reconstruction comprising cameras and points in the vicinity of . Using sub-reconstructions, we can approximate the uncertainty computed from a complete reconstruction. The error of our approximation decreases with increasing size of a sub-reconstruction. If we add a camera to a reconstruction, we add at least four observations which influence the Fisher information matrix as
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
which has full rank. Using Woodbury matrix identity
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 sub-reconstruction produces an upper bound on the uncertainty for cameras involved in the sub-reconstruction. The accuracy of the upper bound depends on a particular decomposition of the complete reconstruction into sub-reconstructions. 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 sub-reconstructions are shown in Figs. 5, 5(a), 5(b). All experiments were performed on a single computer with one 2.6GHz Intel Core i7-6700HQ with 32GB RAM running a 64-bit Windows 10 operating system.
|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|
are listed in Table 2.
The standard way of computing the covariance matrix is by using the M-P 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 . 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 , which uses Eigen library  internally (runs 25.9 minutes on Daliborka dataset) and Alg. 3 is our Matlab implementation, which internally calls LAPACK library  (runs 0.45 seconds on Daliborka dataset). Further, we compared Lhuilier  and Polic  approaches, which approximate the uncertainty propagation, with our algorithm denoted as Nullspace bounding uncertainty propagation (NBUP).
|1.||M-P inversion of using Maple (Kanatani ) (Ground Truth)|
|2.||M-P inversion of using Ceres (Kanatani )|
|3.||M-P inversion of using Matlab (Kanatani )|
|4.||M-P inversion of Schur complement matrix with correction term (Lhuillier )|
|5.||TE inversion of Schur complement matrix with three points fixed (Polic )|
|6.||Nullspace bounding uncertainty propagation (NBUP)|
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 .
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 1-4 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
The error shows the differences between GT covariance matrices and the computed ones . The matrix
has dimension and normalises the error to percentages of the absolute magnitude of the original units. Symbol stands for element-wise 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.
of the algorithms is shown in Fig. 5. Note that the M-P inversion (i.e. Alg. 1-3) cannot be evaluated on medium and larger datasets 5-9 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 7-9 where our NBUP method is about 8 times faster.
on sub-reconstructions was tested on datasets 5-9. We decomposed reconstructions several times using a different number of cameras inside smaller sub-reconstructions, and measured relative and absolute errors of approximated covariances for cameras parameters. Fig. 6 shows the decrease of error for larger sub-reconstructions. There were 25 sub-reconstructions 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. 42. Fig. 5(b) shows that the absolute covariance error decreases significantly with increasing the number of cameras in a sub-reconstruction.
Fig. 5 shows the error of the simplest approximation of covariances used in practice. For every camera, one hundred of its neighbours using view-graph were used to get a sub-reconstruction 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 .
Current methods for evaluating of the uncertainty , 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 Moore-Penrose inversion , which cannot be used in case of medium and large-scale 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  to obtain a full-rank matrix  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 sub-reconstruction of roughly 100-300 cameras provides reliable estimates of the uncertainties for arbitrarily large scenes.
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), EU-H2020 project LADIO no. 731970, and by Grant Agency of the CTU in Prague projects SGS16/230/OHK3/3T/13, SGS18/104/OHK3/1T/37.
-  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://ceres-solver.org
-  Albl, C., Kukelova, Z., Pajdla, T.: R6p-rolling 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 high-performance 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 3-d 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., Fite-Georgel, 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.: Geo-referenced 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 8-10, 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.: Structure-from-motion 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://theia-sfm.org
-  Tian, Y.: The moore-penrose 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)