Homography from two orientation- and scale-covariant features

TL;DR

This paper introduces a novel homography estimation method leveraging orientation- and scale-covariant features like SIFT, deriving new constraints that enable minimal two-point estimation, reducing computational effort and improving stability.

Contribution

It presents new geometric constraints on scales and rotations, leading to a minimal two-point homography solver that enhances efficiency and robustness in geometric estimation tasks.

Findings

Requires only two feature correspondences for homography estimation

Reduces RANSAC iterations significantly compared to four-point methods

Demonstrates effectiveness on synthetic and real-world datasets

Abstract

This paper proposes a geometric interpretation of the angles and scales which the orientation- and scale-covariant feature detectors, e.g. SIFT, provide. Two new general constraints are derived on the scales and rotations which can be used in any geometric model estimation tasks. Using these formulas, two new constraints on homography estimation are introduced. Exploiting the derived equations, a solver for estimating the homography from the minimal number of two correspondences is proposed. Also, it is shown how the normalization of the point correspondences affects the rotation and scale parameters, thus achieving numerically stable results. Due to requiring merely two feature pairs, robust estimators, e.g. RANSAC, do significantly fewer iterations than by using the four-point algorithm. When using covariant features, e.g. SIFT, the information about the scale and orientation is given at no cost. The proposed homography estimation method is tested in a synthetic environment and on publicly available real-world datasets.

Tables

Table 1. Table 1 : The theoretical computational complexity of the solvers. The operations in the solvers ( 1 1 1 st row – steps), the computational complexity of one estimation ( 2 2 2 nd – 1 iter), the correspondence number required for the estimation ( 3 3 3 rd – m 𝑚 m ), possible outlier ratios ( 4 4 4 th – 1 − μ 1 𝜇 1-\mu ), the iteration number required for RANSAC with the confidence set to 0.95 0.95 0.95 ( 5 5 5 th – # iters), and computation complexity of the full procedure ( 6 6 6 th – # comps).

\cellcolorblack!10	\cellcolorblack!102SIFT				\cellcolorblack!103ORI [3]				\cellcolorblack!104PT [15]
steps	$6\times 9$ SVD + $6\times 6$ QR + $4\times 4$ EIG				$6\times 9$ SVD				$8\times 9$ SVD
1 iter	$6\ast 9^2+6^3+4^3=766$				$6\ast 9^2=486$				$8\ast 9^2=649$
$m$	2				3				4
1 - $\mu$	0.25	0.50	0.75	0.90	0.25	0.50	0.75	0.90	0.25	0.50	0.75	0.90
# iters	6	16	71	458	8	34	292	4603	12	71	1177	46 049
# comps	4 596	12 256	54 386	350 828	3 888	16 524	141 912	2 237 058	7 788	46 079	763 873	29 885 801

Table 2. Table 2 : Homography estimation on the AdelaideRMF ( 18 18 18 pairs; 43 43 43 planes) and Multi-H ( 4 4 4 pairs; 33 33 33 planes) and Strecha datasets ( 852 852 852 planes) by GC-RANSAC  [ 7 ] combined with minimal methods. Each column reports the results of a method. The required confidence was set to 0.95 0.95 0.95 . The reported properties are the mean re-projection error ( ϵ italic-ϵ \epsilon , in pixels); the number of samples drawn by GC-RANSAC (# iters.); and the processing time in seconds. Average of 100 100 100 runs on each image pair.

		\cellcolorGray 2SIFT	3ORI [3]	4PT [15]
	$ϵ$ (px)	\cellcolorGray 1.57	1.97	1.61
AdelaideRMF	x# iters.	\cellcolorGray 877	9 772	26 082
($43$#)	time (s)	\cellcolorGray 0.092	0.918	2.989
	$ϵ$ (px)	\cellcolorGray 1.90	3.41	1.87
Multi-H	x# iters.	\cellcolorGray 80 031	458 800	410 781
($33$#)	time (s)	\cellcolorGray 57.921	213.900	300.645
	$ϵ$ (px)	\cellcolorGray 1.42	1.51	1.25
Strecha	x# iters.	\cellcolorGray 4 718	17 414	60 973
($852$#)	time (s)	\cellcolorGray 1.435	3.180	10.246

Equations

$$\displaystyle\begin{array}[]{lllllll}a_{1}&=&\frac{\partial u_{2}}{\partial u_{1}}=\frac{h_{1}-h_{7}u_{2}}{s},&&a_{2}&=&\frac{\partial u_{2}}{\partial v_{1}}=\frac{h_{2}-h_{8}u_{2}}{s},\\[5.69054pt] a_{3}&=&\frac{\partial v_{2}}{\partial u_{1}}=\frac{h_{4}-h_{7}v_{2}}{s},&&a_{4}&=&\frac{\partial v_{2}}{\partial v_{1}}=\frac{h_{5}-h_{8}v_{2}}{s},\end{array}$$

$$u_1{h}_1+v_1{h}_2+h_3-u_1{u}_2{h}_7-v_1{u}_2{h}_8-u_2{h}_9$$

$$u_1{h}_4+v_1{h}_5+h_6-u_1{v}_2{h}_7-v_1{v}_2{h}_8-v_2{h}_9$$

$$h_1-\left(u_2+a_1{u}_1\right)h_7-a_1{v}_1{h}_8-a_1{h}_9$$

$$h_2-\left(u_2+a_2{v}_1\right)h_8-a_2{u}_1{h}_7-a_2{h}_9$$

$$h_4-\left(v_2+a_3{u}_1\right)h_7-a_3{v}_1{h}_8-a_3{h}_9$$

$$h_5-\left(v_2+a_4{v}_1\right)h_8-a_4{u}_1{h}_7-a_4{h}_9$$

$$\text{A}={\mathbf{J}}_2{\mathbf{J}}_1^{-1},$$

$${\text{J}}_i={\text{R}}_i{\text{U}}_i=\left[\begin{array}{cc}\cos ({\alpha }_i)& -\sin ({\alpha }_i)\\ \sin ({\alpha }_i)& \cos ({\alpha }_i)\end{array}\right]\left[\begin{array}{cc}{q}_{u,i}& w_i\\ 0& q_{v,i}\end{array}\right],$$

$${\text{J}}_i^{-1}=\frac{1}{c_i^2{q}_{u,i}{q}_{v,i}+s_i^2{q}_{u,i}{q}_{v,i}}\left[\begin{array}{cc}{s}_i{w}_i+c_i{q}_{v,i}& s_i{q}_{v,i}-c_i{w}_i\\ -s_i{q}_{u,i}& c_i{q}_{u,i}\end{array}\right].$$

$$a_1=\frac{c_2{q}_{u,2}(s_1{w}_1+c_1{q}_{v,1})-s_1{q}_{u,1}(c_2{w}_2-s_2{q}_{v,2})}{q_{u,1}{q}_{v,1}}$$

$$a_2=\frac{c_2{q}_{u,2}(s_1{q}_{v,1}-c_1{w}_1)+c_1{q}_{u,1}(c_2{w}_2-s_2{q}_{v,2})}{q_{u,1}{q}_{v,1}}$$

$$a_3=\frac{s_2{q}_{u,2}(s_1{w}_1+c_1{q}_{v,1})-s_1{q}_{u,1}(s_2{w}_2+c_2{q}_{v,2})}{q_{u,1}{q}_{v,1}}$$

$$a_4=\frac{s_2{q}_{u,2}(s_1{q}_{v,1}-c_1{w}_1)+c_1{q}_{u,1}(s_2{w}_2+c_2{q}_{v,2})}{q_{u,1}{q}_{v,1}}$$

$$q_i=\det {\mathbf{J}}_i=q_{u,i}{q}_{v,i}.$$

$$q_1^2{a}_2{a}_3-q_1^2{a}_1{a}_4+q_1{q}_2=0,$$

$$c_1{s}_2{q}_1{a}_1+s_1{s}_2{q}_1{a}_2-c_1{c}_2{q}_1{a}_3-c_2{s}_1{q}_1{a}_4=0.$$

$$h_8{u}_2{s}_1{s}_2+h_7{u}_2{s}_2{c}_1-h_8{v}_2{s}_1{c}_2-h_7{v}_2{c}_1{c}_2+$$

$$-h_2{s}_1{s}_2-h_1{s}_2{c}_1+h_5{s}_1{c}_2+h_4{c}_1{c}_2=0,$$

$$h_7^2{u}_1^2{q}_2+2h_7{h}_8{u}_1{v}_1{q}_2+h_8^2{v}_1^2{q}_2+h_5{h}_7{u}_2{q}_1+$$

$$-h_4{h}_8{u}_2{q}_1-h_2{h}_7{v}_2{q}_1+h_1{h}_8{v}_2{q}_1+2h_7{h}_9{u}_1{q}_2+$$

$$2h_8{h}_9{v}_1{q}_2+h_2{h}_4{q}_1-h_1{h}_5{q}_1+h_9^2{q}_2=0.$$

$$\mathbf{M}\mathbf{h}=0,$$

$$\mathbf{H}=x{\mathbf{H}}_1+y{\mathbf{H}}_2+{\mathbf{H}}_3,$$

$$\hat{\mathbf{A}}={\mathbf{T}}_2\left[\begin{array}{cc}\mathbf{A}& 0\\ 0& 1\end{array}\right]{\mathbf{T}}_1^{-1},$$

$$t^2{q}_1^2{\hat{a}}_2{\hat{a}}_3-t^2{q}_1^2{\hat{a}}_1{\hat{a}}_4+q_1{q}_2=0,$$

$$h_7^2{u}_1^2{q}_2{t}^2+2h_7{h}_8{u}_1{v}_1{q}_2{t}^2+h_8^2{v}_1^2{q}_2{t}^2+h_5{h}_7{u}_2{q}_1+$$

$$-h_4{h}_8{u}_2{q}_1-h_2{h}_7{v}_2{q}_1+h_1{h}_8{v}_2{q}_1+2h_7{h}_9{u}_1{q}_2{t}^2+$$

$$2h_8{h}_9{v}_1{q}_2{t}^2+h_2{h}_4{q}_1-h_1{h}_5{q}_1+h_9^2{q}_2{t}^2=0.$$

$${\text{p}}_i+\mathbf{\Delta }\left[\begin{array}{c}{\Pi }_x(x,y,z)\\ {\Pi }_y(x,y,z)\end{array}\right]+\left[\begin{array}{cc}\frac{\partial {\Pi }_x^i(x,y,z)}{\partial u}& \frac{\partial {\Pi }_x^i(x,y,z)}{\partial v}\\ \frac{\partial {\Pi }_y^i(x,y,z)}{\partial u}& \frac{\partial {\Pi }_y^i(x,y,z)}{\partial v}\end{array}\right]\left[\begin{array}{c}\Delta u\\ \Delta v\end{array}\right],$$

$${\mathbf{J}}_i=\left[\begin{array}{cc}\frac{\partial {\Pi }_x^i(x,y,z)}{\partial u}& \frac{\partial {\Pi }_x^i(x,y,z)}{\partial v}\\ \frac{\partial {\Pi }_y^i(x,y,z)}{\partial u}& \frac{\partial {\Pi }_y^i(x,y,z)}{\partial v}\end{array}\right],$$

$${\left[\begin{array}{cc}\Delta x& \Delta y\end{array}\right]}^{\text{T}}\approx {\text{J}}_i{\left[\begin{array}{cc}\Delta u& \Delta v\end{array}\right]}^{\text{T}}.$$

$$\frac{\partial {\Pi }_x^i(x,y,z)}{\partial u}=\frac{\partial {\Pi }_x^i(x,y,z)}{\partial x}\frac{x}{\partial u}+$$

$$\frac{\partial {\Pi }_x^i(x,y,z)}{\partial x}\frac{y}{\partial u}+\frac{\partial {\Pi }_x^i(x,y,z)}{\partial x}\frac{z}{\partial u}=\nabla ({\mathbf{\Pi }}_x^i{)}^{\text{T}}{\mathbf{s}}_u,$$

$$\frac{\partial {\Pi }_x^i}{\partial v}=\nabla ({\mathbf{\Pi }}_x^i{)}^{\text{T}}{\text{s}}_v,\frac{\partial {\Pi }_y^i}{\partial u}=\nabla ({\mathbf{\Pi }}_y^i{)}^{\text{T}}{\text{s}}_u,\frac{\partial {\Pi }_y^i}{\partial v}=\nabla ({\mathbf{\Pi }}_y^i{)}^{\text{T}}{\text{s}}_v,$$

$${\mathbf{J}}_i=\left[\begin{array}{c}\nabla ({\mathbf{\Pi }}_x^i{)}^{\text{T}}\\ \nabla ({\mathbf{\Pi }}_y^i{)}^{\text{T}}\end{array}\right]\left[\begin{array}{cc}{\text{s}}_u& {\text{s}}_v\end{array}\right].$$

$$\left[\begin{array}{c}\Delta x_2\\ \Delta y_2\end{array}\right]={\mathbf{J}}_2{\mathbf{J}}_1^{-1}\left[\begin{array}{c}\Delta x_1\\ \Delta y_1\end{array}\right]=\mathbf{A}\left[\begin{array}{c}\Delta x_1\\ \Delta y_1\end{array}\right].$$

Full text

1 Introduction

This paper addresses the problem of interpreting, in a geometrically justifiable manner, the rotation and scale which the orientation- and scale-covariant feature detectors, e.g. SIFT [22] or SURF [10], provide. Then, by exploiting these new constraints, we involve all the obtained parameters of the SIFT features (i.e. the point coordinates, angle, and scale) into the homography estimation procedure. In particular, we are interested in the minimal case, to estimate a homography from solely two correspondences.

Nowadays, a number of algorithms exist for estimating or approximating geometric models, e.g. homographies, using affine-covariant features. A technique, proposed by Perdoch et al. [29], approximates the epipolar geometry from one or two affine correspondences by converting them to point pairs. Bentolila and Francos [11] proposed a solution for estimating the fundamental matrix using three affine features. Raposo et al. [31] and Barath et al. [6] showed that two correspondences are enough for estimating the relative camera motion. Moreover, two feature pairs are enough for solving the semi-calibrated case, i.e. when the objective is to find the essential matrix and a common unknown focal length [9]. Also, homographies can be estimated from two affine correspondences [17], and, in case of known epipolar geometry, from a single correspondence [5]. There is a one-to-one relationship between local affine transformations and surface normals [17, 8]. Pritts et al. [30] showed that the lens distortion parameters can be retrieved using affine features.

Affine correspondences encode higher-order information about the scene geometry. This is the reason why the previously mentioned algorithms solve geometric estimation problems exploiting fewer features than point correspondence-based methods. This implies nevertheless their major drawback: obtaining affine features accurately (e.g. by Affine SIFT [28], MODS [26], Hessian-Affine, or Harris-Affine [24] detectors) is time-consuming and, thus, is barely doable in time-sensitive applications.

Most of the widely-used feature detectors provide parts of the affine feature. For instance, there are detectors obtaining oriented features, e.g. ORB [32], or there are ones providing also the scales, e.g. SIFT [22] or SURF [10]. Exploiting this additional information is a well-known approach in, for example, wide-baseline matching [23, 26]. Yet, the first papers [1, 2, 3, 25, 4] involving them into geometric model estimation were published just in the last few years. In [25], the feature orientations are involved directly in the essential matrix estimation. In [1], the fundamental matrix is assumed to be a priori known and an algorithm is proposed for approximating a homography exploiting the rotations and scales of two SIFT correspondences. The approximative nature comes from the assumption that the scales along the axes are equal to the SIFT scale and the shear is zero. In general, these assumptions do not hold. The method of [2] approximates the fundamental matrix by enforcing the geometric constraints of affine correspondences on the epipolar lines. Nevertheless, due to using the same affine model as in [1], the estimated epipolar geometry is solely an approximation. In [3], a two-step procedure is proposed for estimating the epipolar geometry. First, a homography is obtained from three oriented features. Finally, the fundamental matrix is retrieved from the homography and two additional correspondences. Even though this technique considers the scales and shear as unknowns, thus estimating the epipolar geometry instead of approximating it, the proposed decomposition of the affine matrix is not justified theoretically. Therefore, the geometric interpretation of the feature rotations is not provably valid. A recently published paper [4] proposes a way of recovering full affine correspondences from the feature rotation, scale, and the fundamental matrix. Applying this method, a homography is estimated from a single correspondence in case of known epipolar geometry. Still, the decomposition of the affine matrix is ad hoc, and is, therefore, not a provably valid interpretation of the SIFT rotations and scales. Moreover, in practice, the assumption of the known epipolar geometry restricts the applicability of the method.

The contributions of this paper are: (i) we provide a geometrically valid way of interpreting orientation- and scale-covariant features approaching the problem by differential geometry. (ii) Building on the derived formulas, we propose two general constraints which hold for covariant features. (iii) These constraints are then used to derive two new formulas for homography estimation and (iv), based on these equations, a solver is proposed for estimating a homography matrix from two orientation- and scale-covariant feature correspondences. This additional information, i.e. the scale and rotation, is given at no cost when using most of the widely-used feature detectors, e.g. SIFT or SURF. It is validated both in a synthetic environment and on more than $10\;000$ publicly available real image pairs that the solver accurately recovers the homography matrix. Benefiting from the number of correspondences required, robust estimation, e.g. by GC-RANSAC [7], is two orders of magnitude faster than by combining it with the standard techniques, e.g. four-point algorithm [15].

2 Theoretical background

Affine correspondence $(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{A})$ is a triplet, where $\mathbf{p}_{1}=[u_{1}\;v_{1}\;1]^{\text{T}}$ and $\mathbf{p}_{2}=[u_{2}\;v_{2}\;1]^{\text{T}}$ are a corresponding homogeneous point pair in two images and $\mathbf{A}$ is a $2\times 2$ linear transformation which is called local affine transformation. Its elements in a row-major order are: $a_{1}$, $a_{2}$, $a_{3}$, and $a_{4}$. To define $\mathbf{A}$, we use the definition provided in [27] as it is given as the first-order Taylor-approximation of the $\text{3D}\to\text{2D}$ projection functions. For perspective cameras, the formula for $\mathbf{A}$ is the first-order approximation of the related homography matrix as follows:

[TABLE]

where $u_{i}$ and $v_{i}$ are the directions in the $i$th image ($i\in\{1,2\}$) and $s=u_{1}h_{7}+v_{1}h_{8}+h_{9}$ is the projective depth. The elements of $\mathbf{H}$ in a row-major order are: $h_{1}$, $h_{2}$, …, $h_{9}$.

The relationship of an affine correspondence and a homography is described by six linear equations. Since an affine correspondence involves a point pair, the well-known equations (from $\mathbf{H}\mathbf{p}_{1}\sim\mathbf{p}_{2}$) hold [15]. They are as follows:

[TABLE]

After re-arranging (3), four additional linear constraints are obtained from $\mathbf{A}$ which are the following.

[TABLE]

Consequently, an affine correspondence provides six linear equations for the elements of the related homography.

3 Affine transformation model

In this section, the interpretation of the feature scales and rotations are discussed. Two new constraints that relate the elements of the affine transformation to the feature scale and rotation are derived. These constraints are general, and they can be used for estimating different geometric models, e.g. homographies or fundamental matrices, using orientation- and scale-covariant features. In this paper, the two constraints are used to derive a solver for homography estimation from two correspondences. For the sake of simplicity, we use SIFT as an alias for all the orientation- and scale-covariant detectors. The formulas hold for all of them.

3.1 Interpretation of the SIFT output

Reflecting the fact that we are given a scale $q_{i}\in\mathbb{R}$ and rotation $\alpha_{i}\in[0,2\pi)$ independently in each image ($i\in\{1,2\}$; see Fig. 1), the objective is to define affine correspondence $\mathbf{A}$ as a function of them. For this problem, approaches were proposed in the recent past [3, 4]. None of them were nevertheless proven to be a valid interpretation.

To understand the SIFT output, we exploit the definition of affine correspondences proposed in [8]. In [8], $\mathbf{A}$ is defined as the multiplication of the Jacobians of the projection functions in the two images as follows:

[TABLE]

where $\mathbf{J}_{1}$ and $\mathbf{J}_{2}$ are the Jacobians of the 3D $\to$ 2D projection functions. Proof is in Appendix A. For the $i$th Jacobian, the following is a possible decomposition:

[TABLE]

where angle $\alpha_{i}$ is the rotation in the $i$th image, $q_{u,i}$ and $q_{v,i}$ are the scales along axes $u$ and $v$, and $w_{i}$ is the shear ($i\in\{1,2\}$). Let us use the following notation: $c_{i}=\cos(\alpha_{i})$ and $s_{i}=\sin(\alpha_{i})$. The equation for the inverse matrix becomes

[TABLE]

The denominator can be formulated as follows: $(c_{i}^{2}+s_{i}^{2})q_{u,i}q_{v,i}$, where $c_{i}^{2}+s_{i}^{2}$ is a trigonometric identity and equals to one. After multiplying the matrices in (6), the following equations are given for the affine elements:

[TABLE]

These formulas show how the affine elements relate to $\alpha_{i}$, the scales along axes $u$ and $v$ and shears $w_{i}$.

In case of having orientation- and scale-covariant features, e.g. SIFT, the known parameters are the rotation $\alpha_{i}$ of the feature in the $i$th image and a uniform scale $q_{i}$. It can be easily seen that the scale $q_{i}$ is interpreted as follows:

[TABLE]

Therefore, our goal is to derive constraints that relate affine elements of $\mathbf{A}$ to the orientations $\alpha_{i}$ and scales $q_{i}$ of the features in the first and second images. We will derive such constraints by eliminating the scales along axes $q_{u,i}$ and $q_{v,i}$ and the shears $w_{i}$ from equations (8)-(11). To do this, we use an approach based on the elimination ideal theory [13]. Elimination ideal theory is a classical algebraic method for eliminating variables from polynomials of several variables. This method was recently used in [21] for eliminating unknowns from equations that are not dependent on input measurements. Here, we use the method in a slightly different way. We first create the ideal $I$ [13] generated by polynomials (8)-(11), polynomial (12) and trigonometric identities $c_{i}^{2}+s_{i}^{2}=1$ for $i\in\{1,2\}$. Note that here we consider all elements of these polynomials, including $c_{i}$ and $s_{i}$, as unknowns. Then we compute generators of the elimination ideal $I_{1}=I\cap\mathbb{C}[a_{1},a_{2},a_{3},a_{4},q_{1},q_{2},s_{1},c_{1},s_{2},c_{2}]$ [13]. The generators of $I_{1}$ do not contain $q_{u,i}$, $q_{v,i}$ and $w_{i}$. The elimination ideal $I_{1}$ is generated by two polynomials:

[TABLE]

Generators (13)-(14) can be computed using a computer algebra system, e.g. Macaulay2 [14]. The input code for Macaulay2 is in the supplementary material. The new constraints relate the elements of $\mathbf{A}$ to the scales and rotations of the features in both images. Note that both these equations can be divided by $q_{1}\neq 0$. After this simplification, (13) corresponds to $\det\mathbf{A}=q_{2}/q_{1}=q$ and equation (14) relates the rotations of the features to the elements of $\mathbf{A}$. The two new constraints are general, and they can be used for estimating different geometric models, e.g. homographies or fundamental matrices, using orientation- and scale-covariant detectors. Next, we use (13)-(14) to derive new constraints on a homography.

4 Homography from two correspondences

In this section, we derive new constraints that relates $\mathbf{H}$ to the feature scales and rotations in the two images. Then a solver is proposed to estimate $\mathbf{H}$ from two SIFT correspondences based on these new constraints. Finally, we discuss how the widely-used normalization of the point correspondences [16] affects the output of orientation- and scale-covariant detectors and subsequently the new constraints.

4.1 Homography and covariant features

First, we derive constraints that relates the homography $\mathbf{H}$ to the scales and rotations of the features in the first and second images. To do this, we combine constraints (13) and (14) derived in previous section with the constraints on the homography matrix (5).

Constraints (13) and (14) cannot be directly substituted into (5). However, we can use a similar approach as in the previous section for deriving (13) and (14). First, ideal $J$ generated by six polynomials (5), (13) and (14) is constructed. Then the unknown elements of the affine transformation $\mathbf{A}$ are eliminated from the generators of $J$. We do this by computing the generators of $J_{1}=J\cap\mathbb{C}[h_{1},\dots,h_{9},u_{1},v_{1},u_{2},v_{2},q_{1},q_{2},s_{1},c_{1},s_{2},c_{2}]$. The elimination ideal $J_{1}$ is generated by two polynomials:

[TABLE]

The input code for Macaulay2 used to compute these generators is provided as supplementary material.

Polynomials (15) and (16) are new constraints that relate the homography matrix to the scales and rotations of the features in the first and second images. These constraints will help us for recovering $\mathbf{H}$ from two orientation- and scale-covariant feature correspondences.

4.2 2-SIFT solver

Constraint (15) is linear in the elements of $\mathbf{H}$. For two SIFT correspondences, two such equations are given, which, together with the four equations for point correspondences (4), result in six homogeneous linear equations in the nine elements of $\mathbf{H}$. In matrix form, these equations are:

[TABLE]

where $\mathbf{M}$ is a ${6}\times{9}$ coefficient matrix and $\mathbf{h}$ is a vector of 9 elements of homography matrix $\mathbf{H}$. For two SIFT correspondences in two views, coefficient matrix $\mathbf{M}$ has a three-dimensional null space. Therefore, the homography matrix can be parameterized by two unknowns as

[TABLE]

where $\mathbf{H}_{1},\mathbf{H}_{2},\mathbf{H}_{3}$ are created from the 3D null space of $\mathbf{M}$ and $x$ and $y$ are new unknowns. Now we can plug the parameterization (18) into constraint (16). For two SIFT correspondences, this results in two quadratic equations in two unknowns. Such equations have four solutions and they can be easily solved using e.g. the Gröbner basis or the resultant based method [13]. Here, we use the solver based on Gröbner basis method that can be created using the automatic generator [19]. This solver performs Gauss-Jordan elimination of a ${6}\times{10}$ template matrix which contains just monomial multiples of the two input equations. Then the solver extracts solutions to $x$ and $y$ from the eigenvectors of a ${4}\times{4}$ multiplication matrix that is extracted from the template matrix. Finally, up to four real solutions to $\mathbf{H}$ are computed by substituting solutions for $x$ and $y$ to (18).

Note that we do not know any degeneracies of the proposed solver which can occur in real life. For instance, the degeneracy of the four-point algorithm, i.e. the points are co-linear, is not a degenerate case for the 2SIFT solver.

4.3 Normalization of the affine parameters

The normalization of the point coordinates is a crucial step to increase the numerical stability of $\mathbf{H}$ estimation [16]. Suppose that we are given a $3\times 3$ normalizing transformation $\mathbf{T}_{i}$ transforming the center of gravity of the point cloud in the $i$th image to the origin and its average distance from it to $\sqrt{2}$. The formula for normalizing $\mathbf{A}$ is as follows [6]:

[TABLE]

where $\widehat{\mathbf{A}}$ is the normalized affinity. Matrix $\mathbf{T}_{i}$ transforms the points by translating them (last column) and applying a uniform scaling (diagonal). Due to the fact that the last column of $\mathbf{T}_{i}$ has no effect on the top-left $2\times 2$ sub-matrix of the normalized affinity, the equation can be rewritten as follows: $\widehat{\mathbf{A}}=\text{diag}(t_{2},t_{2})\;\mathbf{A}\;\text{diag}(1/t_{1},1/t_{1})=t_{2}/t_{1}\mathbf{A},$ where $t_{1}$ and $t_{2}$ are the scales of the normalizing transformations in the two images. Thus, for normalizing the affine transformation, it has to be multiplied by $t_{2}/t_{1}$.

The scaling factor affects constraint (13) which, for $\widehat{\mathbf{A}}$, has the form

[TABLE]

where $t={t_{1}}/{t_{2}}$ and $\widehat{a_{i}}$ are elements of $\widehat{\mathbf{A}}$. Consequently constraint (16) for the normalized coordinates has the form

[TABLE]

Note that this normalization does not affect the structure of the derived 2SIFT solver. The only difference is that, for the normalized coordinates, the coefficients in the template matrix are multiplied by scale factor $t$ as in (21).

5 Experimental results

In this section, we compare the proposed solver (2SIFT) with the widely-used normalized four-point (4PT) algorithm [15] and a method using three oriented features [3] (3ORI) for estimating the homography.

5.1 Computational complexity

First, we compare the computational complexity of the competitor algorithms, see Table 1. The first row consists of the major steps of each solver. For instance, $6\times 9$ SVD + $6\times 6$ QR + $4\times 4$ EIG means, that the major steps are: the SVD decomposition of a $6\times 9$ matrix, the QR decomposition of a $6\times 6$ matrix and the eigendecomposition of a $4\times 4$ matrix. In the second row, the implied computational complexities are summed. In the third one, the number of correspondences required for the solvers are written. The fourth row lists example outlier ratios in the data. In the fifth one, the theoretical number of iterations of RANSAC [15] is written for each outlier ratio with confidence set to $0.99$. The last row shows the computational complexity, i.e. the complexity of one iteration multiplied by the number of iteration, of RANSAC combined with the minimal methods. It can be seen that the proposed method leads to significantly smaller computational complexity. Moreover, we believe that by designing a specific solver to our two quadratic equations in two unknowns, similarly as in [20], the computational complexity of our solver can be even reduced.

5.2 Synthesized tests

To test the accuracy of the homographies obtained by the proposed method, first, we created a synthetic scene consisting of two cameras represented by their $3\times 4$ projection matrices $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$. They were located randomly on a center-aligned sphere. A plane with random normal was generated in the origin and ten random points, lying on the plane, were projected into both cameras. The points were at most one unit far from the origin. To get the ground truth affine transformations, we calculated homography $\mathbf{H}$ by projecting four random points from the plane to the cameras and applying the normalized DLT [15] algorithm. The local affine transformation of each correspondence was computed from the ground truth homography by (3). Note that $\mathbf{H}$ could have been calculated directly from the plane parameters. However, using four points promised an indirect but geometrically interpretable way of noising the affine parameters: adding noise to the coordinates of the four projected points. To simulate the SIFT orientations and scales, $\mathbf{A}$ was decomposed to $\mathbf{J}_{1}$, $\mathbf{J}_{2}$. Since the decomposition is ambiguous, $\alpha_{1}$, $q_{u,1}$, $q_{v,1}$, $w_{1}$ were set to random values. $\mathbf{J}_{1}$ was calculated from them. Finally, $\mathbf{J}_{2}$ was calculated as $\mathbf{J}_{2}=\mathbf{A}\mathbf{J}_{1}$. Zero-mean Gaussian-noise was added to the point coordinates, and, also, to the coordinates which were used to estimate the affine transformations.

Fig. 3 reports the numerical stability of the methods in the noise-free case. The frequencies (vertical axis), i.e. the number of occurrences in $100\;000$ runs, are plotted as the function of the $\log_{10}$ average transfer error (in px; horizontal) computed from the estimated homography and the not used correspondences. It can be seen that all tested solvers are numerically stable. Fig. 4 plots the $||\textbf{H}_{\text{est}}-\textbf{H}_{\text{gt}}||_{\text{F}}$ errors as the function of image noise level $\sigma$ (vertical axis) and the ratio (horizontal) of the camera distance, i.e. the radius of the sphere on which the cameras lie, and the object size. The homographies were normalized. The proposed 2SIFT algorithm (left) is less sensitive to the choice of both parameters than the 3ORI (middle) and 4PT (right) methods.

Fig. 5 reports the re-projection error (vertical; in pixels) as the function of the image noise $\sigma$ with additional noise added to the SIFT orientations (left) and scales (right) besides the noise coming from the noisy affine transformations. In the top row, the error is plotted as the function of the image noise $\sigma$. The curves show the results on different noise levels in the orientations and scales. In the bottom row, the error is plotted as the function of the orientation (left plot) and scale (right) noise. The noise in the point coordinates was set to $1.0$ px. The scale noise for the left plot was set to $1\%$. The orientation noise for the right one was set to [math]. It can be seen that, even for large noise in the scale and orientation, the new solver performs reasonably well.

5.3 Real world tests

To test the proposed method on real-world data, we downloaded the AdelaideRMF111cs.adelaide.edu.au/ hwong/doku.php?id=data, Multi-H222web.eee.sztaki.hu/ dbarath, Malaga333www.mrpt.org/MalagaUrbanDataset and Strecha444https://cvlab.epfl.ch/ datasets. AdelaideRMF and Multi-H consist of image pairs of resolution from $455\times 341$ to $2592\times 1944$ and manually annotated (assigned to a homography or to the outlier class) correspondences. Since the reference point sets do not contain rotations and scales, we detected points applying the SIFT detector. The correspondences provided in the datasets were used to estimate ground truth homographies. For each homography, we selected the points out of the detected SIFT correspondences which are closer than a manually set inlier-outlier threshold, i.e. $2$ pixels. As robust estimator, we chose GC-RANSAC [7] since it is state-of-the-art and its implementation is available555https://github.com/danini/graph-cut-ransac. GC-RANSAC is a locally optimized RANSAC with PROSAC [12] sampling. For fitting to a minimal sample, GC-RANSAC used one of the compared methods, e.g. the proposed one. For fitting to a non-minimal sample, the normalized 4PT algorithm was applied.

Given an image pair, the procedure to evaluate the estimators on AdelaideRMF and Multi-H is as follows: first, the ground truth homographies, estimated from the manually annotated correspondence sets, were selected one by one. For each homography: (i) The correspondences which did not belong to the selected homography were replaced by completely random correspondences to reduce the probability of finding a different plane than what was currently tested. (ii) GC-RANSAC was applied to the point set consisting of the inliers of the homography and outliers. (iii) The estimated homography is compared to the ground truth one estimated from the manually selected inliers.

The Strecha dataset consists of image sequences of buildings. All images are of size $3072\times 2048$. The methods were applied to all possible image pairs in each sequence. The Malaga dataset was gathered entirely in urban scenarios with a car equipped with several sensors, including a high-resolution camera and five laser scanners. $15$ video sequences are provided and we used every $10$th image from each sequence. The ground truth projection matrices are provided for both datasets. To get a reference correspondence set for each image pair in the Strecha and Malaga datasets, first, calculated the fundamental matrix from the ground truth camera poses provided in the datasets. SIFT detector was applied. Correspondences were selected for which the symmetric epipolar distance was smaller than $1.0$ pixel. RANSAC was applied to the filtered correspondences finding the most dominant homography with a threshold set to $1.0$ pixel and confidence to $0.9999$. The inliers of this homography were considered as a reference set. In case of having less then $50$ reference points, the pair was discarded from the evaluation. In total, $852$ image pairs were tested in the Strecha dataset and $9\;064$ pairs in the Malaga dataset.

Example results are shown in Fig. 2. The inliers of the homography estimated by estimated by 2SIFT are drawn. Also, the number of iteration required for 2SIFT and 4PT and the ground truth inlier ratios are reported. In all cases, 2SIFT made significantly fewer iterations than 4PT.

Table 2 reports the results on the AdelaideRMF (rows 2–4), Multi-H (5–7) and Strecha (8–10) datasets. The names of the datasets are written into the first column and the numbers of planes are in brackets. The names of the tested techniques are written in the first row. Each block, consisting of three rows, shows the mean re-projection error computed from the manually annotated correspondences and the estimated homographies ($\epsilon$; in pixels; avg. of $100$ runs on each pair); the number of samples drawn by the outer loop of GC-RANSAC ($\#$ iters.); and the processing time (in secs). The RANSAC confidence was set to $0.95$ and the inlier-outlier threshold to $2$ pixels. It can be seen that the proposed method has similar errors to that of the 4PT algorithm, but 2SIFT leads to 1–2 orders of magnitude speedup compared to 4PT.

The results on the Malaga dataset are shown in Figure 6. The confidence of GC-RANSAC was set to $0.95$ and the inlier-outlier threshold to $2.0$ pixels. The reported properties are the average re-projection error (left; in pixels), processing time (middle; in seconds) and the average number of iterations (right). It can be seen that the re-projection errors of 4PT and 2SIFT are fairly similar However, 2SIFT is significantly faster in all cases due to making much fewer iterations than 4PT.

6 Conclusion

We proposed a theoretically justifiable interpretation of the angles and scales which the orientation- and scale-covariant feature detectors, e.g. SIFT or SURF, provide. Building on this, two new general constraints are proposed for covariant features. These constraints are then exploited to derive two new formulas for homography estimation. Using the derived equations, a solver is proposed for estimating the homography from two correspondences. The new solver is numerically stable and easy to implement. Moreover, it leads to results superior in terms of geometric accuracy in many cases. Also, it is shown how the normalization of the point correspondences affects the rotation and scale parameters. Due to requiring merely two feature pairs, robust estimators, e.g. RANSAC, do significantly fewer iterations than by using the four-point algorithm. The method is tested in a synthetic environment and on publicly available real-world datasets consisting of thousands of image pairs. The source code is uploaded as supplementary material.

Appendix A Proof the affine decomposition

We prove that decomposition $\mathbf{A}=\mathbf{J}_{2}\mathbf{J}_{1}^{-1}$, where $\mathbf{J}_{i}$ is the Jacobian of the projection function w.r.t. the directions in the $i$th image, is geometrically valid. Suppose that a three-dimensional point $\mathbf{P}=\begin{bmatrix}x&y&z\end{bmatrix}^{\text{T}}$ lying on a continuous surface $S$ is given. Its projection in the $i$th image is $\textbf{p}_{i}=\begin{bmatrix}u_{i}&v_{i}\end{bmatrix}^{\text{T}}$. The projected coordinates, $u_{i}$ and $v_{i}$, are determined by the projection functions $\mathbf{\Pi}_{u},\mathbf{\Pi}_{v}:\mathbb{R}^{3}\to\mathbb{R}$ as follows: $u_{i}=\mathbf{\Pi}^{i}_{u}(x,y,z)$, $v_{i}=\mathbf{\Pi}^{i}_{v}(x,y,z),$ where the coordinates of the surface point are written in parametric form as $x=\mathcal{X}(u,v)$, $y=\mathcal{Y}(u,v)$, $z=\mathcal{Z}(u,v)$. It is well-known in differential geometry [18] that the basis of the tangent plane at point P is written by the partial derivatives of $S$ w.r.t. the spatial coordinates. The surface normal n is expressed by the cross product of the tangent vectors $\textbf{s}_{u}$ and $\textbf{s}_{v}$ where $\textbf{s}_{u}=\begin{bmatrix}\frac{\partial\mathcal{X}(u,v)}{\partial u}&\frac{\partial\mathcal{Y}(u,v)}{\partial u}&\frac{\partial\mathcal{Z}(u,v)}{\partial u}\end{bmatrix}^{\text{T}},$ and $\textbf{s}_{v}$ is calculated similarly. Finally, $\textbf{n}=\textbf{s}_{u}\times\textbf{s}_{v}$. Locally, around point P, the surface can be approximated by the tangent plane, therefore, the neighboring points in the $i$th image are written as the first-order Taylor-series as follows:

[TABLE]

where $[\Delta v,\Delta u]^{\text{T}}$ is the translation on surface $S$, and $\Delta x$, $\Delta y$ are the coordinates of the implied translation added to $\textbf{p}_{i}$. It can be seen that transformation $\mathbf{J}_{i}$ mapping the infinitely close vicinity around point $\textbf{p}_{i}$ in the $i$th image is given as

[TABLE]

thus

[TABLE]

The partial derivatives are reformulated using the chain rule. As an example, the first element it is as

[TABLE]

where $\mathbf{\nabla\Pi}_{x}^{i}$ is the gradient vector of $\mathbf{\Pi}_{x}$ w.r.t. coordinates $x$, $y$ and $z$. Similarly,

[TABLE]

Therefore, $\mathbf{J}_{i}$ can be written as

[TABLE]

Local affine transformation $\mathbf{A}$ transforming the infinitely close vicinity of point $\textbf{p}_{1}$ in the first image to that of $\textbf{p}_{2}$ in the second one is as follows:

[TABLE]