A linear method for camera pair self-calibration and multi-view reconstruction with geometrically verified correspondences

TL;DR

This paper presents a linear self-calibration method for camera pairs that improves 3D reconstruction accuracy by verifying correspondences and integrating multiple pair reconstructions into a robust pipeline.

Contribution

A novel linear self-calibration approach for camera pairs with unknown focal lengths, including a correspondence verification method and an integrated 3D reconstruction pipeline.

Findings

Performs comparably to Kruppa equations in rotation error

Improves point correspondence robustness in architectural scenes

Introduces a method to average focal length estimates across camera pairs

Abstract

We examine 3D reconstruction of architectural scenes in unordered sets of uncalibrated images. We introduce a linear method to self-calibrate and find the metric reconstruction of a camera pair. We assume unknown and different focal lengths but otherwise known internal camera parameters and a known projective reconstruction of the camera pair. We recover two possible camera configurations in space and use the Cheirality condition, that all 3D scene points are in front of both cameras, to disambiguate the solution. We show in two Theorems, first that the two solutions are in mirror positions and then the relations between their viewing directions. Our new method performs on par (median rotation error $\Delta R = 3.49^{\circ}$) with the standard approach of Kruppa equations ($\Delta R = 3.77^{\circ}$) for self-calibration and 5-Point algorithm for calibrated metric reconstruction of a camera pair. We reject erroneous image correspondences by introducing a method to examine whether point correspondences appear in the same order along $x, y$ image axes in image pairs. We evaluate this method by its precision and recall and show that it improves the robustness of point matches in architectural and general scenes. Finally, we integrate all the introduced methods to a 3D reconstruction pipeline. We utilize the numerous camera pair metric recontructions using rotation-averaging algorithms and a novel method to average focal length estimates.

Tables

Table 1. Table 1 : A summary of notation, with references to uses in the text

subscripts
M	Metric Reconstruction, e.g. $P_M$
P	Projective Reconstruction, e.g. $P_P$
1(or 2)	Refers to Camera 1(or 2) in camera pair, e.g $P_{P1}$
GT	Ground Truth
superscripts
1(or 2)	Refers to solution 1(or 2) for the second camera, e.g ${𝐱}^{\mathrm{𝟏}}$
Accents, as in ${𝐩}^{\prime }$ and $𝐩$	Discriminate between the 2 solutions for camera 2
P matrix representations
	Metric reconstruction
	Metric reconstruction
	Metric reconstruction
${𝐦}_{𝐢}$	i-row vector of left $3\times 3$ P matrix block
	Projective Reconstruction in canonical form
Simplifications
${𝐂}^{𝐢}\triangleq {𝐂}_{𝐢}$	Appear in some Lemmas
$P^i\equiv P_i$
${\omega }^{\ast i}\triangleq {\omega }_i^{\ast }$

Table 2. Table 2 : Performance of the 5P algorithm and our method in recovering the relative rotations R i ​ j subscript 𝑅 𝑖 𝑗 R_{ij} . We initialize BA with each of the aforementioned methods and do 20 20 20 BA iterations. In BA, the internal parameter matrices K i subscript 𝐾 𝑖 K_{i} are held constant. Dataset castle-P30   [ 43 ] .

	Median $\Delta R{(}^{\circ })$	(pairs)	(pairs)
Proposed Method	$3.49$	$123$	$106$
5P	$3.77$	$128$	$104$

Table 3. Table 3 : Precision and Recall for geometric verification of tentative correspondences in Dataset  [ 45 ] . In small italics we give results on the complete set and, in the usual font, on the subset of architectural scenes. We applied the recursive verification method with Threshold T 𝑇 T as a percentage of image size. We give the percentage values α 𝛼 \alpha for T 𝑇 T , following Eq. ( 22 )

Performance Measure	$\alpha =0.02$	$\alpha =0.04$	$\alpha =0.06$	$\alpha =0.08$	$\alpha =0.10$	$\alpha =0.15$	$\alpha =0.20$
Precision	0.99 0.98	0.98 0.97	0.97 0.96	0.96 0.96	0.95 0.95	0.91 0.91	0.87 0.87
Recall	0.80 0.64	0.89 0.72	0.93 0.76	0.94 0.78	0.96 0.80	0.95 0.80	0.96 0.81

Table 4. Table 4 : Precision and Recall for geometric verification of tentative correspondences in Dataset  [ 43 ] . We applied the recursive verification method with Threshold T 𝑇 T as a percentage of image size. We give the percentage values α 𝛼 \alpha for T 𝑇 T , following Eq. ( 22 )

Performance Measure	$\alpha =0.02$	$\alpha =0.04$	$\alpha =0.06$	$\alpha =0.08$	$\alpha =0.10$	$\alpha =0.15$	$\alpha =0.20$
Precision	0.60	0.70	0.76	0.79	0.82	0.86	0.88
Recall	0.85	0.85	0.85	0.85	0.85	0.85	0.84

Table 5. Table 5 : Focal length averaging on castle-P30   [ 42 ] . We did 300 300 300 RANSAC iterations and used the 8-point algorithm for fundamental matrix computation

Method	Median	Confidence count	Joint confidence count
Mean $\Delta f$ error	$0.28$	$0.17$	$0.07$

Table 6. Table 6 : Mean errors in 𝐭 , R , f 𝐭 𝑅 𝑓 \mathbf{t},R,f estimation in the final reconstruction. Datasets from  [ 42 ] .

Dataset	$\Delta 𝐭{(}^{\circ })$	$\Delta R{(}^{\circ })$	$\Delta f$
castle-P30	3.10	1.06	0.0389
castle-P19	7.35	4.17	0.0586
entry-P10	4.62	4.67	0.2118
Herz-Jesu-P8	1.00	0.68	0.0266
Herz-Jesu-P25	0.41	0.31	0.0049
Fountain-P11	0.44	0.41	0.0095

Equations

$$Q_{\infty }^{\ast }=\left[\begin{array}{cc}{I}_{3\times 3}& 0_3\\ 0_3^{\mathbf{T}}& 0\end{array}\right]$$

$$C^{\ast }=P{Q}^{\ast }{P}^T$$

$${\pi }^{\prime }=H^{-T}\pi$$

$$F$$

$$\hat{\mathbf{a}}$$

$$\hat{A}$$

$$K_i=\left[\begin{array}{ccc}{f}_i& 0& 0\\ 0& f_i& 0\\ 0& 0& 1\end{array}\right]$$

$$P_{M1}$$

$$P_{M2}$$

$$\left[\begin{array}{ccc}{f}_2^2& 0& 0\\ 0& f_2^2& 0\\ 0& 0& 1\end{array}\right]={\omega }_2^{\ast }=P_{P2}H{Q}_{\infty }^{\ast }{H}^T{P}_{P2}^T$$

$$H=\left[\begin{array}{cc}{K}_1& 0\\ {\mathbf{v}}^{\mathbf{T}}& \sigma \end{array}\right]$$

$${\pi }_{\infty ,\mathbf{P}}\triangleq {\left(\begin{array}{cc}\mathbf{p}& 1\end{array}\right)}^T\triangleq {\left(\begin{array}{cccc}{p}_1& p_2& p_3& 1\end{array}\right)}^T$$

$$H=\left[\begin{array}{cc}{K}_1& 0\\ -{\mathbf{p}}^{\mathbf{T}}{K}_1& 1\end{array}\right]$$

$${\omega }_2^{\ast }=P_{P2}\left[\begin{array}{cc}{K}_1{K}_1^T& -K_1{K}_1^T\mathbf{p}\\ -{\mathbf{p}}^{\mathbf{T}}{K}_1{K}_1^T& {\mathbf{p}}^{\mathbf{T}}{K}_1{K}_1^T\mathbf{p}\end{array}\right]P_{P2}^T$$

$$P_{P2}\triangleq \left[\begin{array}{cccc}{p}_{11}& p_{12}& p_{13}& p_{14}\\ p_{21}& p_{22}& p_{23}& p_{24}\\ p_{31}& p_{32}& p_{33}& p_{34}\end{array}\right]$$

$${\mathbf{x}}_{\mathbf{o}}\triangleq \left(\begin{array}{c}{f}_1^2\\ f_2^2\\ f_1^2{p}_1^2+f_1^2{p}_2^2+p_3^2\\ p_3\\ f_1^2{p}_1\\ f_1^2{p}_2\end{array}\right)$$

$$A{\mathbf{x}}_{\mathbf{o}}=\mathbf{b}$$

$$\left(\begin{array}[]{rrrrrrr}p_{21}^{2}+p_{22}^{2}&-1&p_{24}^{2}&-2\,p_{23}p_{24}&-2\,p_{21}p_{24}&-2\,p_{22}p_{24}&-p_{23}^{2}\\ p_{21}p_{31}+p_{22}p_{32}&0&p_{24}p_{34}&-p_{23}p_{34}-p_{24}p_{33}&-p_{21}p_{34}-p_{24}p_{31}&-p_{22}p_{34}-p_{24}p_{32}&-p_{23}p_{33}\\ p_{11}p_{31}+p_{12}p_{32}&0&p_{14}p_{34}&-p_{13}p_{34}-p_{14}p_{33}&-p_{11}p_{34}-p_{14}p_{31}&-p_{12}p_{34}-p_{14}p_{32}&-p_{13}p_{33}\\ p_{11}^{2}+p_{12}^{2}&-1&p_{14}^{2}&-2\,p_{13}p_{14}&-2\,p_{11}p_{14}&-2\,p_{12}p_{14}&-p_{13}^{2}\\ p_{11}p_{21}+p_{12}p_{22}&0&p_{14}p_{24}&-p_{13}p_{24}-p_{14}p_{23}&-p_{11}p_{24}-p_{14}p_{21}&-p_{12}p_{24}-p_{14}p_{22}&-p_{13}p_{23}\\ p_{31}^{2}+p_{32}^{2}&0&p_{34}^{2}&-2\,p_{33}p_{34}&-2\,p_{31}p_{34}&-2\,p_{32}p_{34}&-p_{33}^{2}+1\end{array}\right)$$

$$A_5{\mathbf{x}}_{\mathbf{o}}={\mathbf{b}}_5$$

$$\begin{bmatrix}1&0&0&0&0&{\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}0}&b_{1}\\ 0&1&0&0&0&{\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}0}&b_{2}\\ 0&0&1&0&0&{\color[rgb]{0.0600000000000001,0.46,1}\definecolor[named]{pgfstrokecolor}{rgb}{0.0600000000000001,0.46,1}\pgfsys@color@cmyk@stroke{0.94}{0.54}{0}{0}\pgfsys@color@cmyk@fill{0.94}{0.54}{0}{0}0}&b_{3}\\ 0&0&0&1&0&c&b_{4}\\ 0&0&0&0&1&d&b_{5}\end{bmatrix}$$

$$P_{P1}=\left[\begin{array}{cc}I& 0\end{array}\right],P_{P2}=\left[\begin{array}{cc}{[\mathbf{a}]}_{x}F& \mathbf{a}\end{array}\right]$$

$$f_1^2$$

$$f_2^2$$

$$f_1^2{p}_1^2+f_1^2{p}_2^2+p_3^2$$

$$p_3+c{f}_1^2{p}_2$$

$$f_1^2{p}_1+d{f}_1^2{p}_1$$

$${\mathbf{x}}_{\mathbf{o}}^1$$

$${\mathbf{x}}_{\mathbf{o}}^2$$

$$H_k=\left[\begin{array}{cc}{K}_1^{-1}& 0\\ 0^{\mathbf{T}}& 1\end{array}\right]$$

$$F_P$$

$$F_{GT}^{\prime }$$

$$F_P=ϵF_{GT}^{\prime }$$

$$\hat{\mathbf{a}}$$

$$\hat{A}$$

Full text

††journal: arxiv

1 Introduction

Multi-view geometry (mvg) is a Computer Vision (CV) subfield that attempts to understand the structure of the 3D world given a collection of its images [1]. As the binocular human vision is naturally 3D, the same underlying principles allow the recovery of the 3D world structure in mvg reconstruction methods. However, a prerequisite is to have calibrated cameras, an assumption that is violated in unordered image sets. In this paper we focus on self-calibration and multi-view reconstruction using relations between camera pairs.

Assuming a camera pair with unknown and different focal lengths as the only unknown internal parameters, a standard approach to self-calibration and metric reconstruction first applies the 7 point algorithm [1] inside a RANSAC [2] procedure to find the fundamental matrix. In this projective framework, the Kruppa Equations [3] are used to determine the unknown focal lengths. Next, applying the 5 point algorithm [4] inside a RANSAC procedure, leads to a metric reconstruction. Since focal lengths are recovered in a projective framework, only epipolar geometry constraints may be used to check the solution plausibility. Solving self-calibration and metric reconstruction problems simultaneously permits the application of the more intuitive and restrictive geometric arguments of the metric framework.

Self-calibration methods are derived from relations on the dual absolute conic (DAC) $Q^{*}_{\infty}$ and the dual image of the absolute conic (DIAC) $\omega_{\infty}^{*}$ [5, 6, 3]. However, existing methods require three or four images to provide a solution [5, 6], use numerical methods to determine DAC [5], provide an initial DAC estimate that violates the rank-2 condition [5, 6] and do not examine the relations between the recovered putative solutions [5]. In mvg reconstructions additional assumptions have been made to determine focal lengths, as availability of EXIF tags [7, 8], equality of focal lengths across all images [9, 10] and vanishing points correspondences [11].

Towards a multi-view reconstruction, camera pairs have been utilized. Different estimates for a rotation matrix $R$ can be combined with a rotation averaging algorithm [12] and reconstructions of pairs of images can be combined with rotation registration methods [13, 14, 15] to initialize an instance of Structure-from-Motion with known rotations [16, 7] and produce a multiple-view reconstruction.

Erroneous solutions in mvg problems are directly caused by erroneous or noisy image correspondences. Two complementary approaches, applying RANSAC procedures to repeateadly sample minimal point sets and verifying the initial point coresspondences, have been utilised to improve the validity of the recovered solution [17, 18].

In this paper, we derive a linear method for the self-calibration and metric reconstruction of camera pairs with unknown and different focal lengths, unifying two problems that were previously solved independently, to a single system of equations. We further disambiguate the two solutions recovered by our method through the derivation of two theorems about the solutions’ relations. We improve the robustness and applicability of this method by introducing a procedure to verify tentative point correspondences between images, using the Longest Common/Increasing Subsequence (LCS/LIS) problem [19]. The verification method is tailored for outdoors scenes of buildings, and is based on enforcing expected geometric properties of such scenes. We integrate our afforementioned methods to a multi-view reconstruction pipeline, utilizing $L_{\infty}$-norm algorithms and introducing a method to average different estimates of a single focal length $f_{i}$, which uses the structure of the problem, specifically that each estimate for $f_{i}$ comes from a pair of images $i,j$ and is so paired with a second estimate $f_{j}$.

The rest of this article is organized as follows: Section 2 provides background on the reconstruction problem and the verification of image correspondences. Section 3 introduces our method for self-calibration and metric reconstruction. Section 4 presents our method for correspondences verification between the images, first Section 4.1 presents a method which is reduced to the LCS problem and then Section 4.2 presents the final practical algorithm. In Section 5, we integrate our methods to a reconstruction pipeline. In doing so, we develop novel averaging methods for $R,f$ estimates recovered from image pairs. Results for camera pair reconstruction, correspondences verification, focal length averaging and mv reconstruction are given in Section 6.

2 Background & Related Work

In the following bold font (e.g. $\mathbf{v}$ ) is used for vectors and capital case normal font (e.g. $K$) is reserved for matrices.

2.1 Elements of multiple view geometry

In this section we summarize basic notions about the projection of 3D scenes to 2D planes [1, 20]. In a metric reconstruction parallel world lines converge at the plane at infinity $\mathbf{\pi}_{\infty}$: $\begin{pmatrix}0&0&0&1\end{pmatrix}^{T}$. The absolute conic $\Omega_{\infty}$ is a conic on $\mathbf{\pi}_{\infty}$ which satisfies $X_{1}^{2}+X_{2}^{2}+X_{3}^{2}=0$, $X_{4}=0$, where $\mathbf{X}=\begin{pmatrix}X_{1}&X_{2}&X_{3}&X_{4}\end{pmatrix}^{T}$ is the homogeneous representation of world points.

By taking all the planes tangent to $\Omega_{\infty}$, we construct $Q_{\infty}^{*}$, which is the dual surface of $\Omega_{\infty}$. $Q_{\infty}^{*}$ is described in a metric reconstruction by the $4\times 4$ matrix

[TABLE]

Now, considering projective reconstructions of 3-space and projections to image plane we have the following Results [1]:

Result 1**.**

The projection of $Q^{*}$ by projection matrix $P$ in the image plane is the dual conic

[TABLE]

Result 2**.**

If the 3-space is transformed by homography $H$, that is $\mathbf{X^{\prime}}=H\mathbf{X}$, then planes of 3-space are transformed according to

[TABLE]

Result 3**.**

If $H$ is a $4\times 4$ matrix representing a projective transformation of 3-space, then the fundamental matrices corresponding to the pairs of camera matrices $\{P,P^{\prime}\}$, $\{PH,P^{\prime}H\}$ are the same.

Result 4**.**

Suppose the $\operatorname{rank}2$ matrix $F$ can be decomposed in two different ways as

[TABLE]

then

[TABLE]

for some non-zero constant $\kappa$ and 3-vector $\mathbf{v}$

Using the preceding Results, we formulate the equations to solve the camera self-calibration problem and to determine $\mathbf{\pi_{\infty}}$ position in a projective reconstruction.

2.2 Notation

We summarize our notation in Table 1.

2.3 Verifying point correspondences

In standard approaches to find image correspondences regions of interest are described by local feature descriptors [21]. Consequently, erroneous matches occur between similar in local appearance image regions.

A way to reject erroneous matches is using arguments about the geometry of the depicted scenes (geometric verification). In SIFT features, Hough transform was used to acquire the orientation of the detected features [21]. Another common approach applies a rudimentary transform (e.g affine, similarity) between the images, to reject some correspondences before fitting the full model [22, 23, 24, 25]. For these methods we mention that:

• •

Most require specific image features to be extracted, from which special parameters are used to fit the transform

• •

Result in rejection of a large number of correspondences

• •

When we tried using a similarity transform to geometrically verify matches in the reconstruction problem, our results did not improve

Another direction is to improve the covariance of local feature descriptors [24, 17, 26]. The regions of interest can be first transformed before extracting a feature descriptor [24] or ellipses may be matched instead of points [17].

Finally, the neighborhoods of putative matched points are examined in some verification methods. Such approaches include counting the number of correspondences between the neighborhoods of two tentative point matches [27] or examining the order of matched features between the neighborhoods and counting the number of features out of order [28]. The geometric verification method we propose uses properties concerning the order of matched points as well.

2.4 Approaches to multiple view reconstruction

In a reconstruction pipeline, initially Structure from Motion (SfM) is solved to get $P,\quad\mathbf{X}$ assuming image point correspondences and self-calibrated cameras. The fundamental method to solve SfM is Bundle Adjustment (BA) [29], an iterative, numerical algorithm to minimize the reprojection error of the recovered solution.

In standard approaches to SfM a sequence of SfM sub-problems are solved (sequential SfM) [8, 18, 30]. In each iteration, more, possibly uncalibrated, cameras and world points are added to the SfM problem which is solved using BA. However such methods are sensitive to the initial camera pair selection, solve a large number of optimization problems numerically and optimize an objective function with possibly multiple local minima.

A different approach has been developed for solving the SfM with known Rotations problem within the framework of optimal algorithms in multiple-view geometry (mvg) and $L_{\infty}$ mvg algorithms [31, 32, 16, 7, 33, 34]. In this formulation, the camera rotation matrices $R$ are given. SfM is formulated as a convex-optimization problem, for which a unique global minimum exists. For the actual solution of SfM with known rotations, either a sequence of Second-order cone programs are solved to arrive at an exact solution, or approximate solutions are recovered by solving SOCP or Linear programs [9, 35, 7, 11]. BA may still be applied as a last fine-tuning of the solution.

A SfM solution, allows the reconstruction of a low number of 3D points (sparse point cloud), limited by the number of image points correspondences. Multi-view stereo (mvs) algorithms can be used at this point to produce a dense point cloud, which contains a much larger number of 3D points [36]. Finally, surface reconstruction algorithms can be used to produce a 3D surface [37].

3 A method for Metric Reconstruction in pairs of Uncalibrated Images

3.1 Formulation of System Equations

Let us consider two cameras $P_{1},P_{2}$ and further that $P_{1}$ coordinate system is aligned with the world coordinate system. Let us further assume, that the corresponding image coordinate systems are selected so that the internal parameters of each camera $K_{i}$ can be written as

[TABLE]

where $f_{i}$ is the focal length. The previous assumptions are routinely employed in multiple view geometry and are thoroughly discussed in the literature [1].

We start from a projective reconstruction of the 2 cameras, given by $P_{P1},P_{P2}$, which is related to the metric reconstruction by a world (3D) homography $H$ as in

[TABLE]

Using Result 1, Eq.(1), we project $Q_{\infty}^{*}$ to the image plane of camera 2. For this projection, $\omega_{2}^{*}$ we have:

[TABLE]

To introduce the unknowns in Eq. (3), we use the canonical representation of the projective reconstruction, so that $P_{P1}=\begin{bmatrix}I&\mathbf{0}\end{bmatrix}$. From Eq. (2), we have for the homography

[TABLE]

where $\mathbf{v}$ is yet undetermined and the scale factor $\sigma$ can be ignored ($\sigma=1$).

To fully determine $H$, we turn to the plane at infinity

[TABLE]

Using Result 2 we arrive at

[TABLE]

Substituting $H$ from Eq. (4) to Eq. (3) we get

[TABLE]

Eq. (5) comprise a non-linear system with respect to the five unknowns (plane at infinity coordinates and focal lengths) we pursue to determine to acquire a metric reconstruction of the scene. We note that $\omega^{*}_{\infty}$ is symmetric by definition, and is also homogeneous, thus it provides five independent equations.

3.2 Linearization

In Eq. (5), we substitute

[TABLE]

We group the unknowns in the following complexes

[TABLE]

The augmented matrix $[A|\mathbf{b}]$ for the linear system

[TABLE]

is then given by

[TABLE]

We derived the above equations (in order of appearance) from elements $\omega_{2}^{*}(2,2)$, $\omega_{2}^{*}(2,3)$, $\omega_{2}^{*}(1,3)$, $\omega_{2}^{*}(1,1)$, $\omega_{2}^{*}(1,2)$, $\omega_{2}^{*}(3,3)$ of $\omega_{2}^{*}$. In the following, we use the first five equations as explained in Section 3.4.

The matrix of Eq. (8) is rank deficient. Thus, we presented a linear system of five (in the best case) linearly-independent equations, in six unknowns. To solve it, we turn to the polynomial relations between the coordinates of $\mathbf{x_{o}}$.

3.3 Recovering the solutions

Taking five of Eqs. (7) we have the linear system

[TABLE]

Applying Gaussian elimination to (9), we bring the augmented matrix to the form

[TABLE]

where:

The elements in default font, are in the usual form expected when we apply Gaussian elimination in the general case 2. 2.

The elements in green font, are a result of the problem’s structure, that is of the special relations in Eq. (10) 3. 3.

Finally, the element in blue font, is as given when we use the canonical representation for the projective reconstruction, which is:

[TABLE]

Where $\mathbf{a}$ is the left null vector of $F$, $F^{T}\mathbf{a}=\mathbf{0}$. By using the canonical pair, the leftmost $3\times 3$ block in $P_{P2}$ is rank $2$, and consequently has linearly-dependent row-vectors

The derivation of Eq. (10) is given in the Supplementary Material.

To solve for the focal lengths ($f_{1},f_{2}$) and $\mathbf{\pi_{\infty}}$ ($p_{1},p_{2},p_{3}$), we now have from (10)

[TABLE]

We substitute $p_{1}$, from (16), and $p_{3}$, from (15), to Eq. (14), and obtain a second-order equation with respect to $p_{2}$. Thus, we determine $f_{1},f_{2}$ uniquely and $p_{1},p_{2},p_{3}$ with a two-way ambiguity. We refer to those two solutions as

[TABLE]

3.4 The effect of homogeneous representation on the derived equations

In homogeneous coordinate systems representations equal up to a multiplicative constant refer to the same entity. We explore here how this ambiguity affects the formulation of Eq. (7).

Let

• •

$\{P_{GT1}^{\prime},P_{GT2}^{\prime}\}\triangleq\{\begin{bmatrix}I&\mathbf{0}\end{bmatrix},\begin{bmatrix}A&\mathbf{a}\end{bmatrix}\}$, be the ground truth camera matrices we aim to recover

• •

$\{P_{P1},P_{P2}\}\triangleq\{\begin{bmatrix}I&\mathbf{0}\end{bmatrix},\begin{bmatrix}\hat{A}&\hat{\mathbf{a}}\end{bmatrix}\}$, be the starting projective reconstruction

$P_{GT1}^{\prime}$ is related to $P_{GT1}$ by the homography

[TABLE]

Thus, we get from $P_{P}$ to $P_{GT}$ by homography $H$, from $P_{P}$ to $P_{GT}^{\prime}$ by $H^{\prime}$ and from $P_{GT}$ to $P_{GT}^{\prime}$ by $H_{k}$.

For the camera pairs, we have the Fundamental matrices

[TABLE]

From Result 3, since $P_{GT}^{\prime},P_{P}$ reconstructions are related by $H^{\prime}$, the reconstructions share a common Fundamental matrix. Since Fundamental matrices are homogeneous entities, we have

[TABLE]

Now, we turn to Result 3 and get

[TABLE]

We write the previous equations in matrix form to get the projective transformation $H^{\prime}$

[TABLE]

$H^{\prime}$ satisfies

[TABLE]

Now, we get $H$ from $H^{\prime},H_{k}^{-1}$

[TABLE]

We set the bottom-right element to $1$, as we disregard the true scale of the reconstruction, and get the final form of $H$

[TABLE]

One should compare the homographies of Eq. (19) and Eq. (4). Using Eq. (4) instead of Eq. (19), we get from $P_{P}$ to

[TABLE]

We observe that the translation direction $\mathbf{a}$ is correct but the left-most $3\times 3$ block of camera $2$ is multiplied by a constant $\mu$.

To see how the constants in Eq. (19) affect Eq. (7), we substitute Eq. (19) in Eq. (3) and get for $\omega_{2}^{*}$ the expression

[TABLE]

To avoid the determination of additional unknowns in Eq. (7), we have

• •

All equations derived from $\omega_{2}^{*}$ elements off the diagonal are of the form $\mathbf{ax_{o}=0}$, thus the constant $(\kappa\epsilon)^{-1}$ can be eliminated

• •

The equation derived from element $\omega_{2}^{*}(3,3)$ cannot be used without determining additional constants. So, we may only use the rest five of the six original equations of (7)

The complete method to solve the metric reconstruction and self calibration problem follows:

We solve the system (7), keeping five equations and discarding the equation derived from $\omega_{2}^{*}(3,3)$ 2. 2.

In the previous step (1), we recovered $c^{-1}f_{2}^{2}$. To fully determine $f_{2}$, many different approaches are possible. We propose to repeat step 1, putting camera 2 at the origin of the coordinate system (in place of camera 1). This can be done by transposing the Fundamental matrix $F$ for the camera pair. Following this approach, we may additionally determine the constant $\kappa\epsilon$ 3. 3.

Using the homography of Eq. (4) or Eq. (19), we recover the metric reconstruction $P_{M1},P_{M2}$. Depending on the homography used, one camera matrix ($P_{M2}$ for Eq. (4) or $P_{M1}$ for Eq. (19)) will have the left-most $3\times 3$ block multiplied by a constant. This has no effect on the correctness of the representation, and the image points are the same in each case

3.5 Solution disambiguation and geometric relations of the two solutions

We use the Cheirality condition (Corollary 1) to determine the valid solution of Eq. (7). Whenever the two recovered solutions represent cameras with divergent viewing directions, Cheirality condition is more likely to identify the valid solution. We explore in Theorems 1 and 2 the geometric relations between the two solutions, aiming to vizualize solutions’ relations and disambiguation. Proofs of Theorems 1 and 2 are outlined in Figs. 2 and 3. Full proofs are given in the Supplemental Material.

Theorem 1**.**

Let

[TABLE]

denote the reconstructions derived from (17). Then, cameras $P_{m2}^{1},P_{m2}^{2}$ are in mirror positions with respect to the origin (position of $P_{m1}^{1},P_{m1}^{2}$). The centers of projection $\mathbf{C_{m2}^{1},C_{m2}^{2}}$ satisfy

[TABLE]

Theorem 2**.**

Let camera $1$ be positioned on the origin of the world coordinate system, with a viewing direction aligned to $z$ axis. We denote $\mathbf{v_{m2}^{1},v_{m2}^{2}}$ the viewing directions of $P_{m2}^{1},P_{m2}^{2}$ and $\mathbf{C_{m2}^{1},C_{m2}^{2}}$ the position vectors of the corresponding centers of projection. Then, $\mathbf{C_{m2}^{1},C_{m2}^{2}}$ bisect the angles formed by $\mathbf{v_{m2}^{1},v_{m2}^{2}}$, in the plane defined by $\mathbf{v_{m2}^{1},v_{m2}^{2}}$. Thus, we have:

[TABLE]

From Theorems 1,2, we easily deduce that

[TABLE]

Corollary 1**.**

The correct one of solutions (17) can be identified by requiring all world points that are visible from camera $2$ to be in the space in front of camera $2$.

4 Geometric verification of tentative image correspondences

4.1 Reduction to Longest Common Subsequence problem

The geometric property we pursued to enforce in tentative image correspondences is the order of imaged points with respect to the horizontal and vertical image directions. We have:

• •

If a point, A, is imaged to the left of a point, B, in the first image, then A should be to the left of B in the second image as well. We call this property Consistency-x

• •

Similarly, if a point, A, is below another point, B, in the first image, then A should be below B in the second image as well. We call this property Consistency-y

To see how we can arrive at the LCS/LIS problem we examine each one of the two Consistency properties independently. We present here the analysis concerning Consistency-x.

We start with a formal definition of Consistency-x. A set of correspondences $S=\{(p_{1}^{i},p_{2}^{i})\}$, where $p_{j}^{i}$ is the x-coordinate of a point ($i$) in image $j$, has Consistency-x, if for all points $p_{1}^{i}$ of image $1$ in $S$:

All points in image $1$ that are in the Consistent-x set and are to the left of $p_{1}^{i}$ match in image $2$ with points that are to the left of $p_{2}^{i}$:

[TABLE] 2. 2.

All points in image $1$ that are in the Consistent-x set and are to the right of $p_{1}^{i}$ match in image $2$ with points that are to the right of $p_{2}^{i}$:

[TABLE]

We seek the most populous set $S$ of correspondences which is Consistent-x. We can reduce the Consistency-x problem to LIS in the following way: We sort points $p_{1}^{i}$ in image 1 with respect to the x-axis ($x_{i}$). This sorting is a permutation in the sequence of correspondences. We apply this same permutation to $p_{2}^{i}$ and get a sequence from the ordinates $x_{i}$ of points $p_{2}^{i}$. We seek the LIS of this last sequence.

The LCS/LIS problems are efficiently solved with complexity $\mathcal{O}(n\log n)$ [38, 39, 19] or even $\mathcal{O}(n\log\log n)$ [40] if special data structures are implemented. To solve LCS/LIS, we used the patience sorting algorithm [38].

4.1.1 Perplexities of the combined Consistency-x,y problem and an efficient approximate method

In the combined Consistency-x,y problem we seek to find the largest subset of correspondences which are consistent in both x and y axes. The relation “Consistency-x and Consistency-y” is not transitive. We can observe that easily with a counter-example. Thus, the Consistency-x,y relation is not a partial order, a condition sufficient to rule out reduction to LCS/LIS [19].

Formally, in the Consistency-x,y problem, we seek a set $S_{xy}$ of image correspondences so that:

• •

$S_{xy}$ has the Consistency-x property

• •

$S_{xy}$ has the Consistency-y property

• •

The number of elements ($n$) of the set $S_{xy}$ is maximized

We propose an approximate solution by the following method:

We find the largest Consistent-x subset, $S_{x}$, solving an LIS problem 2. 2.

We find the largest Consistent-y subset of $S_{x}$, solving again an LIS problem

In our suboptimal solution Consistency-x,y holds, thus our primary aim to reject erroneous matches is achieved. Nevertheless, some true matches are rejected.

4.2 A practical verification method

The consistency properties we introduced, depend on assumptions on the geometric structure of the scene. In photos of architectural scenes, usually the $y-$axis in camera coordinate system aligns with the perpendicular to the floor vector, leading to the assumption $R_{z}=R_{x}=0$. In special cases, as in photographs of houses on a street, the camera $x-$axis may also be aligned between photographs.

Such assumptions may be violated between different views. The effects of camera rotations on a scene are illustrated in Fig. 4. We observe that:

• •

Lines parallel to $x-$ or $y-$ axis in one image may appear tilted in another, if the camera coordinate systems are not aligned. The same effect is caused by scene depth variation

• •

The relative order of points may change between two images. Moreover, it is more likely for two points to change order with respect to the $x$-axis, if those points are close in $x-$ axis but distant in $y$ axis, $\Delta x\ll\Delta y$.

Still, in the case of photographs of architectural scenes, we can assume small rotations around the $x,y$ axes, as the photographer’s position in space is constrained. In-plane ($R_{z}$) rotations are uncommon and can nevertheless be fixed automatically [41].

4.2.1 Approximations to Consistency properties

As Consistency properties are violated by projective phenomena, enforcing them leads to the rejection of many true correspondences. Thus, we relax Consistency properties to arrive at a practical verification method. We describe the method concerned with the order of points on the $x-$axis. Similar modifications apply to the Consistency-y property.

First, we introduce a threshold value ($T$) to allow violations in the order of points that remain within a predefined distance range. So, two consecutive sequence points $s_{i}$, $s_{i+1}$ are considered in correct order if

[TABLE]

where $s_{i}$ is the $x-$coordinate of the i-th point in the ordered sequence. We set $T$ as a fraction of the maximum distance in the $y-$ axis, of any two points in the image we examine, that matched to points in the paired image:

[TABLE]

In the following we refer to this process as “setting the threshold as a percentage of image size”.

We propose to use a recursive method, acting on image subregions of different size. We solve a sequence of LIS problems, each one with a different $T$ value, set as a percentage of image size:

We solve an LIS problem using a threshold $T$ as a percentage of image size. The result is a Consistent-x set of correspondences 2. 2.

We split the image in two subregions, each with equal number of correspondences $n_{i}$. The split is done on the $y-$axis 3. 3.

(Recursion): We repeat the process on each of the two subregions. We terminate if the region size is smaller than a predefined constant $c$ (we used $c=200$ pixels)

The recursive method has the advantage of allowing for larger violations in the order in the $x-$axis for points that are distant in the $y-$axis, as explained in Section 4.2. Concerning the computational complexity, we have:

[TABLE]

where $L$ is the number of recursion steps. $L$ depends on the initial image size and $c$. Consequently, the recursive method adds no significant computational burden to the initial LIS problem formulation.

Finally, we remark that other approaches, as dropping recursion or fitting a simple transform to map lines between the images to estimate the $T$ value, produced worse results than the proposed recursive method.

5 An application to the multiple-view reconstruction problem

We integrate our methods for the geometric verification of image correspondences and the pair-based estimation of $R,f$, in existing pipelines to solve the multiple-view reconstruction problem and produce a 3D-model of a scene.

Our approach is outlined in Fig. 5. The final reconstruction is done using the non-sequential SfM with known rotations formulation of [7], which we modify extensively, using the methods of the preceeding sections as well as the $f,R$ averaging algorithms we describe in the following.

5.1 Averaging pair-based solutions for $f,R$

In this paper we introduced computationally efficient methods for $R_{ij},f_{i}$ estimation, which we apply in randomly sampled minimal correspondences sets, in a way that resembles RANSAC procedures [2]. The multiple $f_{i}$, $R_{ij}$ estimates, one from each minimal sample, are then averaged, to produce the final solutions.

In $f_{i}$, case, we introduce a novel averaging method. In the case of pairwise rotations $R_{ij}$, we apply the Weiszfeld algorithm [12, 15], which converges to the median ($L_{1}$-average) rotation. We also use a form of the Weiszfeld algorithm (multiple rotation averaging) in the rotation registration problem to get the final camera rotation matrices $R_{i}$ (Section 5.1.2).

5.1.1 Focal length estimates

The distribution of $f_{i}$ estimates collected from all the possible image pairs $i,j$ can be skewed or multimodal (Fig. 6), in which case the mean or median estimate will not correctly determine $f_{i}$ value.

We introduce new measures to evaluate the fit of focal length estimates. We initially introduce the Confidence count (cc) and then modify cc using the problem structure to introduce the Joint confidence count (Jcc). We assume that in each image pair that contains image $i$, we receive a number of correct and a number of erroneous estimates for $f_{i}$, and that erroneous estimates originating from different image pairs vary significantly in value, whereas correct ones aggregate.

We visualize cc computation in Fig. 7. Simplifying aspects of the computation, we can describe it as a binning procedure, where the bin range is adapted to contain all estimates within $\beta\%$ deviation:

We collect all $f_{i}^{n}$ estimates of $f_{i}$, originating from all the different images we have matched with image $i$ 2. 2.

For each $f_{i}^{n}$, we count the number of estimates, $f_{i}^{k}$, within a $\beta\%$ error range. This sum is the confidence count ${cc}_{i}^{n}$ for estimate $f_{i}^{n}$ 3. 3.

We normalize ${cc}_{i}^{n}$ values to $0\dots 1$ range. This step is critical for Jcc computation

To further improve the $f$ estimation, we introduce Jcc (Fig. 8). Since each estimate $f_{i}^{n}$ is paired with some estimate $f_{j}^{n}$ (the estimates were computed in an image pair), we expect that if $f_{i}^{n}$ is a good estimate then $f_{j}^{n}$ will be accurate too. To compute Jcc, we follow a similar to cc procedure, but this time each estimate $f_{i}^{k}$ in $\beta\%$ range contributes a different amount to Jcc sum. This amount is proportional to $cc_{j}^{k}$ of estimate $f_{j}^{k}$ that is paired with $f_{i}^{k}$. Good $f_{j}^{k}$ estimates have higher confidence counts, and contribute more to Jcc.

In greater detail, to compute the ${Jcc}_{i}^{n}$ of estimate $f_{i}^{n}$ about image $i$, we have:

Let $k=1\dots m$ be the $m$ images we matched with image $i$. For each image $k$ we have:

• •

From all estimates within $\beta\%$ range of $f_{i}^{n}$, we pick the $L$ ones that originate from pair $i,k$.

• •

Since every $f_{i}$ estimate originating from $i,k$ pair is matched to an $f_{k}$ estimate, from the $L$ estimates of the previous step we get the corresponding $L$ estimates of $f_{k}$

• •

For each of the $L$ estimates of $f_{k}$, we have a confidence count ${cc}_{k}^{n}$. We get their mean. We do not use the direct sum, to diminish the influence of a large sum (large $L$) of low cc’s. 2. 2.

${Jcc}_{i}^{n}$ is the sum of the previous $m$ mean values.

5.1.2 Rotation estimates

In this section, we summarize rotation averaging using the Weiszfeld algorithm [12, 15]. Weiszfeld algorithm returns the $L_{1}$-mean in a set of points in space $\mathcal{R}^{n}$. Many different metrics have been defined for rotation matrices [15]. We limit our analysis here to

[TABLE]

Weiszfeld algorithm is a gradient-descent method and is guaranteed to converge to the true $L_{1}$-mean in the case of single rotation averaging, as averaging of pairwise rotation estimates $R_{ij}$.

The $L_{1}$-mean of $R_{i}$ estimates of a single rotation is the rotation $R_{y}$ that minimizes:

[TABLE]

In this case of rotation registration, the convergence of Weiszfeld algorithm is not guaranteed.

In detail, we applied Weiszfeld algorithm to weight the estimates $R_{ij}$ of the pairwise rotations we acquired through random sampling of minimal point sets ($8$ points) yielding a $R_{ij}$ solution.

In the rotation registration problem we applied the Weiszfeld algorithm in the following manner:

We construct the rotations graph, with one node for every image and an edge $e_{ij}$ between nodes $i,j$ if we know the relative rotation $R_{ij}$ between the respective images. We take a spanning tree in this graph, and using $R_{j}=R_{ij}R_{i}$ we get the initial $R_{j}$ estimates 2. 2.

For every node $i$ in the graph, we use all available estimates $R_{ij}$ to get inconsistent estimates $R_{i}^{k}$, $k=1,2,\dots$ through $R_{i}=R_{ji}R_{j}$. We average estimates $R^{k}_{i}$ with one iteration of Weiszfeld algorithm 3. 3.

We repeat the previous step $n$ times ($n=20$)

In all our experiments we set $\beta=10\%$.

6 Results & Discussion

6.1 Metric Reconstruction in Pairs of Images

We implemented Kruppa equations, a well-studied and popular method for camera self calibration, and used it as reference method for the estimation of internal camera parameters. To compare the methods, we used synthetic camera projection matrices and image correspondences. We added Gaussian noise to the image points positions, and not to world points or other entities, to simulate actual noisy correspondences.

We observed that our method (Sec. 3) and the Kruppa method produce identical $f_{i}$ estimates. In rare cases with extremely noise-corrupted correspondences, our method failed and the Kruppa method produced largely erroneous focal length estimates.

We conclude that the two methods are equivalent concerning the self-calibration problem. Still our method is advantageous in additionally providing a metric reconstruction.

Next, we evaluate camera pair reconstruction. We compared our method to the 5-Point(5P) algorithm [4]. We used both approaches as initialization to BA [29] and evaluated the quality of the final reconstructions (Tab. 2). We used a multiple-view dataset [43] and determined the relative positions of all camera pairs with point correspondences. The same focal length estimates were used in both compared approaches. $f$ estimates were obtained by the method we introduced in Sec. 3. The two methods we compared were:

• •

Initialize BA with our method: We randomly sampled minimal subsets of correspondences and averaged the acquired solutions with rotation averaging [12]. We allowed for 20 BA iterations

• •

Initialize BA with 5P algorithm: We used a RANSAC procedure to sample minimal 5P subsets and to pick the solution. We allowed for 20 BA iterations

To quantify the reconstruction error, we used the angle ($\Delta R$) between the relative rotation estimate and the true relative rotation, $R_{ij}$, between two paired views $i,j$.

The initialization of BA is important, to improve convergence and to reduce the computational cost. We observe that both the 5P algorithm and our method can be used as BA initialization with similar performance (Tab. 2). This result implies that to further reduce the reconstruction error, we should improve other problem parameters as image correspondences and focal length estimates.

6.2 Geometric verification of tentative correspondences

In geometric verification, correspondences are classified as correct or erroneous. We evaluate this classifier using precision and recall. Two different datasets were used:

• •

Dataset [43]: The set contains outdoor scenes of landmark buildings. The ground truth camera matrices, $P$, are provided, from which we can separate correct and erroneous correspondences. In detail, from given $P_{i},P_{j}$ we recover the fundamental matrix $F_{ij}$ and then evaluate Sampson’s approximation to geometric error for each tentative correspondence [44]

• •

Dataset [45]: This set contains both architectural scenes and scenes with objects. We give performance results on each of those two subsets independently. In this set, point correspondences between images are provided and labeled as correct or erroneous

The results are presented in Tables 3, 4. Precision of the classifier is more important than its recall, as it is more important to have an oulier-free set of correspondences than to recover all true correspondences. Furthermore, the recursive verification method discards erroneous matches with very high precision. This result supports our argument that points more distant in the one, e.g $x$, axis are more likely to violate order with respect to the other, e.g $y$, axis. Finally, performance varies with scene type. The development of our method was based on scene properties found in architectural scenes. In scenes composed of objects, differences in scene geometry and the increased freedom in viewer’s position cause more violations in the Consistency properties. In Dataset [43], the performance in scenes of one main building, and consequently of a single main horizontal and vertical direction, as Fountain-P11, entry-P10, Herz-Jesu-P8, reaches almost flawless precision ($1$). In more complex scenes, which include more buildings, as in castle-P30, the achieved precision degrades to values in the range $0.7-0.8$.

6.3 Improving focal length estimation in multi-view reconstructions

We show in Tab. 5 the improved $f$ estimates we get with cc. Further improvement is achieved by Jcc measure. To quantify the error in $f$ estimation we use $\Delta f$ [46, 47, 48]:

[TABLE]

In Fig. 10 we present an non-ambiguous Jcc distribution from which $f$ can be correctly determined, in contrast to $f$ estimates distributions displayed in Fig. 6.

6.4 Multi-view reconstruction in unordered image sets

In Tab. 6, we provide quantitative performance measures for multi-view reconstructions that were acquired applying the proposed pipeline (Fig. 5), on unordered image datasets and with no other input apart from the scene photographs. To quantify reconstruction error in camera translation, we used the angle ($\Delta\mathbf{t}$) between the translation estimate and the available true translation $\mathbf{t}$. In Fig. 11 we qualitatively display the results of the proposed reconstruction pipeline. The results in Tab. 6 and Fig. 11 demonstrate that the introduced methods can be used in unordered image sets to produce quality reconstructions of the photographed scenes.

7 Conclusions

Using the DIAC, we developed a linear self-calibration and metric reconstruction method. Two theorems describe the relative configuration of the two recovered solutions and provide support to use the Cheirality condition for solution disambiguation. Comparisons to Kruppa equations and the 5P algorithm revealed that our method performs similarly to these standard approaches. Subsequently we show that the large number of $f,R$ estimates that are produced by our self-calibration and metric reconstruction method can be utilised through averaging methods, shifting our focus from choosing the best solution to finding, eg as in finding an optimised $f$ estimate prior to self-calibration, the best solution averaging method. We also developed a general method to verify point matches between images, which can be solved by reduction to LCS. The corresponding verification method can be used in any problem with image correspondences input. The verification method successfully rejected outliers in both architectural and general scenes, with more success in the former category. All our methods were integrated to a full multiple-view reconstruction pipeline to produce visually high-quality reconstructions on both standard datasets and image sets we shot using a conventional camera. Multi-view reconstructions were obtained combining camera pair reconstructions using rotation averaging algorithms and a novel approach to average focal length estimates.