Radially Distorted Homographies, Revisited

TL;DR

This paper introduces a unified, fast, and accurate method for estimating homographies with radial lens distortion in images, addressing three distinct distortion configurations simultaneously, and demonstrating superior speed over existing methods.

Contribution

It presents a novel unified approach and minimal solvers for radially distorted homographies applicable to three different distortion scenarios, improving speed while maintaining accuracy.

Findings

Solvers are faster than state-of-the-art methods.

Maintains similar accuracy to existing approaches.

Effective on fisheye camera images.

Abstract

Homographies are among the most prevalent transformations occurring in geometric computer vision and projective geometry, and homography estimation is consequently a crucial step in a wide assortment of computer vision tasks. When working with real images, which are often afflicted with geometric distortions caused by the camera lens, it may be necessary to determine both the homography and the lens distortion-particularly the radial component, called radial distortion-simultaneously to obtain anything resembling useful estimates. When considering a homography with radial distortion between two images, there are three conceptually distinct configurations for the radial distortion; (i) distortion in only one image, (ii) identical distortion in the two images, and (iii) independent distortion in the two images. While these cases have been addressed separately in the past, the present paper provides a novel and unified approach to solve all three cases. We demonstrate how the proposed approach can be used to construct new fast, stable, and accurate minimal solvers for radially distorted homographies. In all three cases, our proposed solvers are faster than the existing state-of-the-art solvers while maintaining similar accuracy. The solvers are tested on well-established benchmarks including images taken with fisheye cameras. A reference implementation of the proposed solvers is made available as part of HomLib (https://github.com/marcusvaltonen/HomLib).

Tables

Table 1. Table 1 : Properties of the evaluated state-of-the-art solvers including the median runtime for 100 000 100\,000 random problem instances.

Solver	Pts.	$Sols.$	Min.	Time
\SetRowht=1.1 Two-sided with independent $\lambda$ and ${\lambda }^{\prime }$
cvpr/2015/kukelova_etal	5	$5$	$✓$	\qty34.8\micro
cvpr/2015/kukelova_etal	6	$2$		\qty12.3\micro
Proposed	5	$5$	$✓$	\qty8.1\micro
\SetRowht=1.1 Two-sided with ${\lambda }^{\prime }=\lambda$
cvpr/2001/fitzgibbon	5	$18$		\qty64.0\micro
cvpr/2015/kukelova_etal	5∗	$5$		\qty25.6\micro
cvpr/2015/kukelova_etal	6∗	$2$		\qty11.2\micro
Proposed	4.5	$6$	$✓$	\qty3.3\micro
\SetRowht=1.1 One-sided, i.e. ${\lambda }^{\prime }=0$
cvpr/2001/fitzgibbon	5∗∗	$6$		\qty17.6\micro
icpr/2024/nakano	4.5	$3$	$✓$	\qty4.5\micro
Proposed	4.5	$3$	$✓$	\qty0.6\micro

Table 2. Table 2 : Results on HPatches . mAA (AUC score) when allowing the solvers to run up to \qty 1000 \milli (two-sided independent), \qty 250 \milli (two-sided equal), \qty 30 \milli (one-sided). Best results are marked in bold.

Solver	mAA (AUC)
	@\qty1px	$@\text{qty}3px$	@\qty5px
\SetRowht=1.1 Two-sided with independent $\lambda$ and ${\lambda }^{\prime }$
cvpr/2015/kukelova_etal (5pt)	0.0983	$0.3343$	0.5079
cvpr/2015/kukelova_etal (6pt)	0.1704	$0.3693$	0.5276
Proposed (5pt)	0.1261	0.3820	0.5379
\SetRowht=1.1 Two-sided with ${\lambda }^{\prime }=\lambda$
cvpr/2001/fitzgibbon (5pt)	0.3092	$0.5133$	0.6310
cvpr/2015/kukelova_etal (5pt)∗	0.1127	$0.3756$	0.5379
cvpr/2015/kukelova_etal (6pt)∗	0.2186	$0.4415$	0.5867
Proposed (4.5pt)	0.3227	0.5192	0.6329
\SetRowht=1.1 One-sided, i.e. ${\lambda }^{\prime }=0$
cvpr/2001/fitzgibbon (5pt)∗∗	0.2542	$0.5996$	0.7285
icpr/2024/nakano (4.5pt)	0.2537	$0.5992$	0.7269
Proposed (4.5pt)	0.2566	0.6019	0.7298

Table 3. Table 3 : Results on Grossmünster Church . Mean error and mean execution time of the solvers in a LOMSAC framework. Best results are marked in bold (across both categories).

Solver	Error	$Exec.time$
\SetRowht=1.1 Two-sided with independent $\lambda$ and ${\lambda }^{\prime }$
cvpr/2015/kukelova_etal (5pt)	0.039	$\text{qty}2657\text{milli}$
cvpr/2015/kukelova_etal (6pt)	0.040	$\text{qty}330\text{milli}$
Proposed (5pt)	0.038	$\text{qty}150\text{milli}$
\SetRowht=1.1 Two-sided with ${\lambda }^{\prime }=\lambda$
cvpr/2001/fitzgibbon (5pt)	0.040	$\text{qty}233\text{milli}$
cvpr/2015/kukelova_etal (5pt)∗	0.040	$\text{qty}2053\text{milli}$
cvpr/2015/kukelova_etal (6pt)∗	0.039	$\text{qty}321\text{milli}$
Proposed (4.5pt)	0.039	\qty84\milli

Equations

$${\alpha }_j{\mathbf{x}}_j^{\prime }=\mathbf{H}{\mathbf{x}}_j,{\alpha }_j\in \mathbb{R},j=1,\dots ,n.$$

$${\mathbf{x}}_4={\gamma }_1{\mathbf{x}}_1+{\gamma }_2{\mathbf{x}}_2+{\gamma }_3{\mathbf{x}}_3,$$

$${\mathbf{x}}_4^{\prime }={\gamma }_1^{\prime }{\mathbf{x}}_1^{\prime }+{\gamma }_2^{\prime }{\mathbf{x}}_2^{\prime }+{\gamma }_3^{\prime }{\mathbf{x}}_3^{\prime }.$$

$$\mathbf{\Xi }=\left(\begin{array}{ccc}{\mathbf{x}}_1& {\mathbf{x}}_2& {\mathbf{x}}_3\end{array}\right)\text{and}{\mathbf{\Xi }}^{\prime }=\left(\begin{array}{ccc}{\mathbf{x}}_1^{\prime }& {\mathbf{x}}_2^{\prime }& {\mathbf{x}}_3^{\prime }\end{array}\right).$$

$$\mathbf{\Gamma }={\mathbf{\Xi }}^{-1}{\mathbf{x}}_4\text{and}{\mathbf{\Gamma }}^{\prime }=({\mathbf{\Xi }}^{\prime }{)}^{-1}{\mathbf{x}}_4^{\prime }.$$

$$\mathbf{H}={\mathbf{\Xi }}^{\prime }(\mathrm{diag}{\mathbf{\Gamma }}^{\prime })(\mathrm{diag}\mathbf{\Gamma }{)}^{-1}{\mathbf{\Xi }}^{-1}.$$

$$\mathbf{M}\mathrm{adj}\mathbf{M}=\det (\mathbf{M})\mathbf{I}$$

$$\mathrm{adj}\mathbf{\Xi }=\left(\begin{array}{ccc}{v}_2{w}_3-v_3{w}_2& u_3{w}_2-u_2{w}_3& u_2{v}_3-u_3{v}_2\\ v_3{w}_1-v_1{w}_3& u_1{w}_3-u_3{w}_1& u_3{v}_1-u_1{v}_3\\ v_1{w}_2-v_2{w}_1& u_2{w}_1-u_1{w}_2& u_1{v}_2-u_2{v}_1\end{array}\right),$$

$$\mathbf{\Gamma }=\mathrm{adj}(\mathbf{\Xi }){\mathbf{x}}_4,{\mathbf{\Gamma }}^{\prime }=\mathrm{adj}({\mathbf{\Xi }}^{\prime }){\mathbf{x}}_4^{\prime },$$

$$\mathbf{H}={\mathbf{\Xi }}^{\prime }(\mathrm{diag}{\mathbf{\Gamma }}^{\prime })(\mathrm{diag}\mathbf{\Gamma }{)}^{-1}\mathrm{adj}(\mathbf{\Xi }).$$

$$w_j$$

$$w_j^{\prime }$$

$$\mathrm{adj}(\mathrm{diag}{\mathbf{\Gamma }}^{\prime })\mathrm{adj}({\mathbf{\Xi }}^{\prime }){\mathbf{x}}_5^{\prime }\sim \mathrm{adj}(\mathrm{diag}\mathbf{\Gamma })\mathrm{adj}(\mathbf{\Xi }){\mathbf{x}}_5,$$

$$\mathbf{N}(\lambda )\sim {\mathbf{N}}^{\prime }(0)⟺\mathbf{N}(\lambda )\times {\mathbf{N}}^{\prime }(0)=0,$$

$$\mathbf{N}(\lambda )\sim {\mathbf{N}}^{\prime }(\lambda )⟺\mathbf{N}(\lambda )\times {\mathbf{N}}^{\prime }(\lambda )=0,$$

$$\mathbf{N}(\lambda )\sim {\mathbf{N}}^{\prime }({\lambda }^{\prime })⟺\mathbf{N}(\lambda )\times {\mathbf{N}}^{\prime }({\lambda }^{\prime })=0,$$

Full text

(i) distortion in only one image, (ii) identical distortion in the two images, and (iii) independent distortion in the two images.

While these cases have been addressed separately in the past, the present paper provides a novel and unified approach to solve all three cases. We demonstrate how the proposed approach can be used to construct new fast, stable, and accurate minimal solvers for radially distorted homographies. In all three cases, our proposed solvers are faster than the existing state-of-the-art solvers while maintaining similar accuracy. The solvers are tested on well-established benchmarks including images taken with fisheye cameras. A reference implementation of the proposed solvers is made available as part of HomLib111https://github.com/marcusvaltonen/HomLib.

1 Introduction

Homographies constitute an important class of transformations that are often used in algorithms for solving computer vision tasks such as precalibration [cvpr/1999/sturm_maybank, pami/2000/zhang, pami/2007/hartley_kang], autocalibration [eccv/1998/triggs], stereo and multi-camera calibration [iccv/2003/ueshiba_tomita, ao/2005/hu_tan, 3dv/2016/zhu_etal, ao/2017/guan_etal], colour calibration [pami/2019/finlayson_etal], metric rectification [cvpr/1998/liebowitz_zisserman, ijcv/2000/criminisi_etal], stereo rectification [cvpr/1993/hartley_gupta, ijcv/1999/hartley, cvpr/1999/loop_zhang, prl/2010/kumar_etal], ego-motion estimation [icra/2002/liang_pears, jfr/2015/zienkiewicz_davison, icip/2016/wadenback_etal, icpram-selected/2019/valtonen-ornhag_wadenback, icpram/2020/valtonen-ornhag, wacv/2021/valtonen-ornhag_etal], panoramic stitching and mosaicing [cga/1996/szeliski, ijcv/2007/brown_lowe, cvpr/2007/brown_etal, jfr/2009/nicosevici_etal, cvpr/2011/gao_etal, cvpr/2013/zaragoza_etal], scene understanding [cvprw/2006/perera_etal, its/2010/arrospide_etal, iccv/2011/flint_etal], visual servoing [ijrr/2007/benhimane_malis, tac/2009/hu_etal], augmented reality [isar/2000/simon_etal, isar/2000/seo_hong, cga/2002/prince_etal], and more.

In most of these situations, the invocation of a homography is justified by the underlying geometry being considered, as a homography describes central projection of one projective plane onto another [book/2004/hartley_zisserman]. For example, consider a perspective pinhole camera viewing a plane in the scene. In this case, the mapping of homogeneous coordinates from the scene plane to the image plane will be a homography, cf. [pami/2000/zhang]. Alternatively, considering two such views of a plane in the scene, the mapping between the two views will, again, be a homography, for the image region that contains the plane [book/2004/hartley_zisserman].

Algebraically, a homography is represented by a (unique up to scale) non-singular matrix $\mathbf{H}\in\mathbb{R}^{3\times 3}$ that maps homogeneous coordinates $\mathbf{x}$ in one plane to $\mathbf{x}^{\prime}$ in the other plane according to $\mathbf{x}^{\prime}\sim\mathbf{H}\mathbf{x}$, where ‘$\sim$’ denotes equality up to scale [book/2011/richter-gebert, book/2004/hartley_zisserman]. To fixate the eight degrees of freedom of $\mathbf{H}$, it is sufficient to specify four point correspondences, i.e. four points and their respective images under $\mathbf{H}$, as long as no three points (or their images under $\mathbf{H}$) are collinear. If this non-collinearity condition is met, the points are said to be in general position. Given at least four such correspondences, one can use the well-known Direct Linear Transformation (DLT) to solve for $\mathbf{H}$ linearly [book/2004/hartley_zisserman]. Since the pinhole perspective camera does not model lenses, it is typically used together with a separate distortion model such as the Brown–Conrady model (see e.g. [book/2022/szeliski]) or the one-parameter division model [dagm/1987/lenz, cvpr/2001/fitzgibbon, ipol/2014/aleman-flores_etal] when applied to real images. For this reason, homography estimation methods based on DLT but which also consider radial distortion have been proposed in the literature [cvpr/2001/fitzgibbon, cvpr/2015/kukelova_etal, icpr/2024/nakano].

Methods for automatically establishing correspondences generally cannot be guaranteed to only produce correct matches, so in practice the set of correspondences will contain outliers in the form of spurious correlations. Such outliers have the potential to completely invalidate any subsequent model estimation, and it is important to reject them before proceeding. This is often done using RANSAC [acm/1981/fischler_bolles] or one of its more refined incarnations [iccv/2003/nister, pami/2012/raguram_etal, cviu/2000/torr_zisserman, cvpr/2005/chum_matas, dagm/2003/chum_etal], where a consensus set of inliers is found by repeatedly fitting the model to small random samples and then checking which of the other data also support this model. The fewer data that are used for fitting the model in RANSAC, the higher is the chance of drawing an uncontaminated sample that only contains inliers, and thus hopefully in the process gains support from most of the other inliers. For this reason, considerable effort has been directed to create so called minimal solvers for different geometric estimation problems, i.e. algorithms that use as few data as possible while ensuring a finite number of solutions. Minimal solvers can provide both likely inlier sets and reasonable initialisations for further refinement such as bundle adjustment [iwva/1999/triggs_etal].

Our contributions.

Inspired primarily by earlier work by cvpr/2015/kukelova_etal and more recent work by icpr/2024/nakano, we revisit the important problem of estimating radially distorted homographies, and demonstrate that by replacing the DLT with a closed-form expression for $\mathbf{H}$ [book/1910/veblen_young, book/1962/seidenberg, book/2011/richter-gebert, prl/2019/guo], we obtain numerically stable minimal solvers that are faster than the current state-of-the-art solvers for all possible configurations, i.e.

(i) the one-sided case, (ii) the two-sided equal case, and (iii) the two-sided independent case.

An illustration of the problem geometry for the three cases is shown in Fig. 2.

Organisation of the paper.

The rest of the paper is organised as follows. In Sec. 2 we discuss relevant related work, including current state-of-the-art methods for radially distorted homographies. We formulate the problem mathematically and recount the classical closed-form formula for homographies in Sec. 3. Our proposed solvers are derived in Sec. 4. The solvers are evaluated in Sec. 5, and Sec. 6 concludes the paper.

2 Related Work

As mentioned in the introduction, homographies are a crucial component of many camera calibration algorithms. For example, cvpr/1999/sturm_maybank use homographies from scene planes to provide constraints on the intrinsic camera parameters. Similarly, pami/2000/zhang computes homographies from a planar calibration target to the image plane using DLT, but in contrast to [cvpr/1999/sturm_maybank] also incorporates radial distortion in a final non-linear refinement step, which forms the basis for the widely used calibration method in OpenCV.

The first work to include radial distortion already at the homography computation step is due to cvpr/2001/fitzgibbon, who extended the DLT equations into a quadratic eigenvalue problem to capture lens distortion modelled by the one-parameter division model [dagm/1987/lenz, cvpr/2001/fitzgibbon].

Subsequently, solvers that make use of additional assumptions have been proposed for some notable special cases, e.g. for panoramic stitching [cvpr/2007/brown_etal, cvpr/2008/jin, bmvc/2009/byrod_etal] where images are assumed to be captured with coinciding camera centres. Instead of eliminating the translation, another option is to assume that the rotation is known, which can make sense in drone applications where an IMU can provide sufficiently accurate rotation estimates [wacv/2021/valtonen-ornhag_etal]. The present paper, however, does not make any such simplifying assumptions for the scene or pose parameters.

We consider three different cases for radially distorted homographies, which we will refer to as the one-sided case, the two-sided equal case, and the two-sided independent case (see Fig. 2). Each of these cases has been considered separately in previous work, but we take a novel unified approach for all three cases. Starting chronologically, cvpr/2001/fitzgibbon solved the *two-sided equal case * by deriving a quadratic eigenvalue problem with up to 18 real solutions. cvpr/2015/kukelova_etal proposed two solvers for the two-sided independent case, i.e. for radially distorted homographies with independent distortion coefficients $\lambda$ and $\lambda^{\prime}$. Their minimal 5 point solver uses a Gröbner basis approach, and requires performing Gauss–Jordan elimination on a \numproduct16 x 21 template matrix. Their non-minimal approach uses 6 points, and is computationally significantly cheaper, but while less sensitive to noise is more sensitive to outliers due to the inclusion of the additional point. Recently, icpr/2024/nakano proposed a DLT-based solution to the one-sided case $\lambda^{\prime}=0$, finding the parameters of the transformation from the null-space to a \numproduct10 x 12 design matrix. For an overview comparison of some of the key properties (including timings) of the existing solvers and our proposed solvers, see Tab. 1.

3 Problem formulation

Let us start by considering the problem without radial distortion. Given a number of point correspondences, $\mathbf{x}_{j}^{\prime}\leftrightarrow\mathbf{x}_{j}$ for $j=1,\ldots,n$, expressed using homogeneous coordinates $\mathbf{x}_{j}=(u_{j},v_{j},w_{j})$ and $\mathbf{x}_{j}^{\prime}=(u_{j}^{\prime},v_{j}^{\prime},w_{j}^{\prime})$, we want to find a non-singular $\mathbf{H}\in\mathbb{R}^{3\times 3}$ such that

[TABLE]

In the computer vision literature, this problem is traditionally solved using the Direct Linear Transformation (DLT) [book/2004/hartley_zisserman, book/2022/szeliski], and this route was used as a foundation for the solvers with radial distortion introduced in [cvpr/2001/fitzgibbon], [cvpr/2015/kukelova_etal] and [icpr/2024/nakano]. We will instead use a much older approach, which we shall discuss next, as a springboard to constructing our proposed solvers.

3.1 The classical closed-form solution

While largely overlooked by the computer vision community, the classical closed-form solution for a homography has a long history and is well established in the literature [book/1910/veblen_young, book/1962/seidenberg, book/2011/richter-gebert]. Its core idea is the fact that it is easy to find the mapping to (or from) the four intermediary points $\{(1,0,0),\,(0,1,0),\,(0,0,1),\,(1,1,1)\}$. More specifically, it works as follows. Let $n=4$ and assume that the points are in general position. This means that we can find non-zero scalars $\gamma_{1},\gamma_{2},\gamma_{3}$ such that

[TABLE]

and likewise, non-zero scalars $\gamma_{1}^{\prime},\gamma_{2}^{\prime},\gamma_{3}^{\prime}$ such that

[TABLE]

Indeed, if one of these scalars were to be zero, that would contravene our assumption that the points are in general position, i.e., three of the points would need to be collinear on at least one side of the transformation.

Now, the homography $\mathbf{U}=\begin{pmatrix}\gamma_{1}\mathbf{x}_{1}&\gamma_{2}\mathbf{x}_{2}&\gamma_{3}\mathbf{x}_{3}\end{pmatrix}$ maps the intermediary points to $\mathbf{x}_{1},\ldots,\mathbf{x}_{4}$, and in the same way, $\mathbf{U}^{\prime}=\begin{pmatrix}\gamma_{1}^{\prime}\mathbf{x}_{1}^{\prime}&\gamma_{2}^{\prime}\mathbf{x}_{2}^{\prime}&\gamma_{3}^{\prime}\mathbf{x}_{3}^{\prime}\end{pmatrix}$ maps the intermediary points to $\mathbf{x}_{1}^{\prime},\ldots,\mathbf{x}_{4}^{\prime}$. Hence, the sought homography is obtained by reversing the first mapping and applying the mappings in succession, yielding the composite mapping $\mathbf{H}=\mathbf{U}^{\prime}\mathbf{U}^{-1}$.

We may write this procedure more compactly by denoting $\mathbf{\Gamma}=(\gamma_{1},\gamma_{2},\gamma_{3})$ and $\mathbf{\Gamma}^{\prime}=(\gamma_{1}^{\prime},\gamma_{2}^{\prime},\gamma_{3}^{\prime})$, as well as

[TABLE]

With this notation, we now have

[TABLE]

Then, finally, $\mathbf{H}$ can be computed in closed-form as

[TABLE]

There are some small and obvious improvements that we shall wish to apply to the basic method described above. First of all, since

[TABLE]

for any square matrix $\mathbf{M}$ [book/1944/aitken], the non-diagonal inverses in 5 and 6 can be replaced with adjugates, as the scalar factor is irrelevant when working with homogeneous coordinates. Recalling the explicit formula for \numproduct3 x 3 adjugates (see e.g. Ch. 3 in [book/1944/aitken]), we may write

[TABLE]

and similarly for $\operatorname{adj}(\mathbf{\Xi}^{\prime})$. These are now plugged into 5 and 6 to give

[TABLE]

and

[TABLE]

Secondly, if all points are finite and normalised to have their last entry equal to one, the resulting expression 10 becomes a little simpler and requires fewer arithmetic operations to evaluate [prl/2019/guo].

3.2 Including radial distortion

Assuming the standard one-parameter division model for radial distortion [dagm/1987/lenz, cvpr/2001/fitzgibbon], we have

[TABLE]

Three distinct cases are of interest here (cf. Fig. 2):

(i) the one-sided case $\lambda^{\prime}=0$, (ii) the two-sided equal case $\lambda^{\prime}=\lambda$, and (iii) the two-sided independent case (where $\lambda$ and $\lambda^{\prime}$ are independent).

4 Deriving the solvers

For the sake of generality, let us initially assume that $\lambda$ and $\lambda^{\prime}$ are independent and potentially both non-zero. Then, for a fifth point correspondence $\mathbf{x}_{5}^{\prime}\leftrightarrow\mathbf{x}_{5}$, we can use 10 to write

[TABLE]

where the left hand side is a vector valued cubic polynomial $\mathbf{N}^{\prime}(\lambda^{\prime})\in\mathbb{R}^{3}[\lambda^{\prime}]$ and the right hand side will be a cubic polynomial $\mathbf{N}(\lambda)\in\mathbb{R}^{3}[\lambda]$ (the expressions for the coefficients will be identical up to toggling of the primes).

Using the explicit formula for the adjugate of a \numproduct3 x 3 matrix again, we can expand both sides as necessary and read off the coefficients of $\mathbf{N}(\lambda)$ and $\mathbf{N}^{\prime}(\lambda^{\prime})$.

4.1 The one-sided case $\lambda^{\prime}=0$

This is the case considered by icpr/2024/nakano, and it arises naturally in connection with camera calibration. Here, the goal is to find the homography from the scene plane (e.g. a planar calibration target) to the rectified scene plane, together with the distortion coefficient $\lambda$. An alternative, equivalent interpretation, is that of finding a homography between a distorted and a rectified image.

Setting $\lambda^{\prime}=0$ in 12 and taking the cross product yields

[TABLE]

which corresponds to three scalar cubic equations in $\lambda$, which will have up to three real-valued solutions. We can solve any one of them using standard techniques, e.g. Cardano’s formula [book/1944/turnbull], the trigonometric method [book/1944/turnbull], the companion matrix method [sjmaa/2015/aurentz_etal] or Sturm sequences [book/1944/turnbull]. For its speed, and to avoid complex numbers, we use the trigonometric method.

4.2 The two-sided equal case

This is the case originally considered by cvpr/2001/fitzgibbon, and arises naturally in situations with two images taken by the same camera (and without changing the optical configuration), for example when taking shots for a panorama. Setting $\lambda^{\prime}=\lambda$ in 12 yields

[TABLE]

which corresponds to three scalar sextic equations in $\lambda$. We can again solve any one of them using a standard technique of our choice, and it will result in up to six real-valued solutions. To avoid computations involving complex numbers, we use Sturm sequences [book/1944/turnbull].

4.3 The two-sided independent case

This is the case considered by cvpr/2015/kukelova_etal, and arises naturally e.g. when stitching images coming from different cameras (or if the optical configuration has changed, e.g. due to zoom or focus adjustments). From 12 we get

[TABLE]

which corresponds to three scalar polynomial equations of degree six in $\lambda$ and $\lambda^{\prime}$.

A polynomial system of equations of this type can be solved e.g. using resultant-based methods or Gröbner basis methods. We use the automatic Gröbner basis generator proposed in [cvpr/2017/larsson_etal] to create a numerical solver. The system has in general nine complex-valued solutions. Attentive readers note that the solver by cvpr/2015/kukelova_etal only had five potential solutions—the four extra solutions obtained from our solver are introduced due to the adjugate matrix computations and can therefore be discarded by checking the determinant of the corresponding matrix. By using the basis heuristic proposed in [cvpr/2018/larsson_etal], an elimination template of size \numproduct9 x 18 was found, which combined with Sturm sequence root finding generates a numerically stable and computationally efficient solver.

5 Experiments

We compare our solvers against the existing state-of-the-art solvers, i.e. icpr/2024/nakano for the one-sided case, cvpr/2001/fitzgibbon for the two-sided equal case, and cvpr/2015/kukelova_etal for the two-sided independent case. All solvers are implemented in C++ using Eigen, and timings are shown in Tab. 1. The results show the median value of $100\,000$ random problem instances, executed on a standard laptop equipped with an 11th Gen Intel® Core™ i5-1145G7 @\qty2.6\giga CPU. While this reporting of the runtime gives some general sense of the speed of the solver, additional aspects, e.g. the sample size and the typical number of solutions, will influence the effective runtime.

5.1 Numerical stability

To evaluate the numerical stability of the solvers, we have generated $10\,000$ scenes with 3D points distributed randomly on an unknown plane. A pair of cameras with a depth of 0.1–10 to the scene plane, focal length = 1000, and distortion coefficients in the range $[-0.2,-0.01]$ were generated such that the field-of-view for both cameras is 70 degrees and the scene points are in front of the cameras. This captures a low to medium level of distortion, for which the one-parameter division model is known to be a good approximation for typical physical cameras.

As can be seen in Fig. 3, all the considered methods are numerically stable. For the solvers with two distortion coefficients, the algebraic error is computed as the geometric mean of the individual ones, i.e. $k_{\text{err}}=\sqrt{k_{1,\text{err}}k_{2,\text{err}}}$.

5.2 Sensitivity to noise

In addition to the numerical stability, we want to investigate the sensitivity of the proposed solvers with respect to Gaussian noise affecting the data. Gaussian noise of varying standard deviation is therefore applied to the distorted (normalised) point correspondences and the corresponding median errors over $10\,000$ random problem instances are reported, see Fig. 4.

Overall, the evaluated solvers all have comparable noise sensitivity to the other solvers for the same case. The strongest deviations are seen for cvpr/2001/fitzgibbon which is somewhat more sensitive with respect to the distortion coefficient in the one-sided case, and for the 6 point version of cvpr/2015/kukelova_etal in the two-sided independent case, which is less sensitive (this is somewhat expected, as this method uses an additional point, cf. [pami/2014/pham_etal]). In situations with both high noise and a high inlier ratio, this could be advantageous, but conversely, if the noise and inlier ratio are low, the inclusion of the extra point may instead be an encumbrance.

5.3 Integration in a robust framework

We integrate the polynomial solvers under evaluation into a RANSAC framework to get a better understanding of the performance. While all solvers are fast, RANSAC also requires evaluating a hypothesis for each model estimated, which affects the total time. We study the cumulative number of of inliers found versus execution time in Fig. 5, where each measurement shows the average over $1000$ random problem instances. Settings in previous sections apply, and the standard deviation of the noise is set to \qty0.5px, which is realistic for traditional descriptors such as SIFT [ijcv/2004/lowe], and an inlier threshold of \qty5px is used. We find that our proposed solvers for the one-sided and two-sided equal cases consistently find more inliers for different outlier ratios. For the two-sided independent case, the 6 point solver by cvpr/2015/kukelova_etal has a clear advantage due to it only having two solutions to evaluate, even though the margin decreases for higher outlier ratios.

5.4 Evaluation on real data

Next, we evaluate the proposed solvers on real images by integrating them in a LOMSAC-framework [bmvc/2012/lebeda_etal] applying a local optimisation step to the so-far-the-best model and non-linear refinement.

Evaluation on HPatches.

For homography estimation, it is common to evaluate on the HPatches dataset [cvpr/2017/balntas_etal], is a Homography benchmark consisting of $116$ planar scenes and comes with an accurate reference in the form of manually annotated correspondences that are often used as ground truth. Unfortunately, HPatches does not have noticeable radial distortion, and the dataset does not contain a reference estimate for this. As a sanity check, we nonetheless report the mAA (AUC scores) on HPatches for thresholds of \qty1px, \qty3px and \qty5px in Tab. 2.

We follow the evaluation settings in LightGlue [iccv/2023/lindenberger_etal] and resize the shorter side of the image to $480$ pixels, and use a total of $1024$ keypoints with SuperPoint [cvprw/2018/detone_etal]. SuperPoint + LightGlue is used for establishing the correspondences.

Detailed scores for the maximum runtime is shown in Tab. 2. Here we allowed the solvers to run up to \qty1 (two-sided independent), \qty250\milli (two-sided equal), and \qty30\milli (one-sided), respectively. mAA (AUC score) as a function of runtime is shown in Fig. 6. In general our proposed solvers yield higher mAA than previous solvers.

Evaluation on Grossmünster Church.

Due to the lack of real-life homography benchmarks with radial distortion we use the Grossmünster Church dataset [eccv/2020/larsson_etal] which contains real images with fisheye distortion. In contrast to the HPatches dataset distortion is present, and is not perfectly approximated by the one-parameter division model. This dataset was originally used in RadialSfM [eccv/2020/larsson_etal], and we use it to generate a 3D reconstruction of the scene, see Fig. 1. In order to extract ground truth homographies, we use the acquired poses, and compute the relative pose $\mathbf{R}=\mathbf{R}_{2}\mathbf{R}_{1}^{\top}$ and relative translation $\mathbf{t}=\mathbf{t}_{2}-\mathbf{R}\mathbf{t}_{1}$. As suggested in [cvpr/2023/barath_etal], we select only images with a single dominant plane—in our case we chose the wall containing the two church towers—and fitted a plane directly on the 3D model. Among the images in the dataset, seven images were selected for which it was possible to extract six consecutive image pairs with more than $100$ tentative image correspondences (including outliers), see Fig. 7. The ground truth homography was then computed as $\mathbf{H}=\mathbf{K}(\mathbf{R}-\mathbf{t}(\mathbf{n}^{\prime})^{\top}/d)\mathbf{K}^{-1}$, where $\mathbf{K}$ is the shared intrinsic parameters across all views, $\mathbf{n}^{\prime}=\mathbf{R}_{1}\mathbf{n}$ is the relative plane normal, and $d$ the corresponding depth.

For the experiment, we extracted SIFT keypoints and applied Lowe’s ratio test. The final correspondences for one image pair are shown in Fig. 8. Since radial distortion is present in both images, only the two-sided solvers are used in this comparison. The homography error is computed as $\epsilon=\lVert\mathbf{H}-\mathbf{H}_{\textup{est}}\rVert_{\textup{F}}$, where both $\mathbf{H}$ and $\mathbf{H}_{\textup{est}}$ are normalised to have unit Frobenius norm and positive determinant.

All solvers are configured to run LOMSAC [bmvc/2012/lebeda_etal] with a minimum of $500$ iterations and final non-linear least squares fitting is applied. The average homography error is almost identical for all solvers; however, the execution times differ, see Tab. 3.

6 Conclusion

In this paper, we have revisited the fundamental problem of simultaneously estimating a homography and radial distortion. Basing our analysis on a classical closed-form solution instead of the DLT, we have taken a unified approach to derive three new, fast, stable and accurate solvers for the possible cases of radially distorted homographies, i.e.

(i) the one-sided case, (ii) the two-sided equal case, and (iii) the two-sided independent case.

We have evaluated our proposed solvers against existing state-of-the-art solvers for each of the three cases, and our new solvers either outperform or compare favourably to the state-of-the-art in these evaluations.

Acknowledgements

This work was supported by the Wallenberg Artificial Intelligence, Autonomous Systems and Software Program (WASP), funded by the Knut and Alice Wallenberg Foundation, as well as by the strategic research environment ELLIIT, funded by the Swedish government.