An Efficient Algebraic Solution to the Perspective-Three-Point Problem

TL;DR

This paper introduces a new algebraic method for solving the perspective-3-point problem, enabling accurate and robust camera pose estimation with lower computational cost by directly deriving attitude and position from geometric constraints.

Contribution

The paper presents an algebraic approach that directly computes camera attitude and position from three reference points, improving accuracy and efficiency over previous methods.

Findings

Higher numerical accuracy demonstrated through simulations

Reduced computational cost compared to existing methods

Robust performance in near-singular configurations

Abstract

In this work, we present an algebraic solution to the classical perspective-3-point (P3P) problem for determining the position and attitude of a camera from observations of three known reference points. In contrast to previous approaches, we first directly determine the camera's attitude by employing the corresponding geometric constraints to formulate a system of trigonometric equations. This is then efficiently solved, following an algebraic approach, to determine the unknown rotation matrix and subsequently the camera's position. As compared to recent alternatives, our method avoids computing unnecessary (and potentially numerically unstable) intermediate results, and thus achieves higher numerical accuracy and robustness at a lower computational cost. These benefits are validated through extensive Monte-Carlo simulations for both nominal and close-to-singular geometric configurations.

Tables

Table 1. Table 1: Nominal case: Pose mean errors.

	position	orientation
Kneip’s method	1.18E-05	1.02E-05
Masselli’s method	1.84E-08	4.89E-10
Proposed method	1.66E-10	5.30E-12
Proposed method+Root polishing	5.07E-11	1.53E-13

Table 2. Table 2: Singular case 1: Pose median errors.

	position	orientation
Kneip’s method	1.42E-14	1.34E-14
Masselli’s method	7.13E-15	6.15E-15
Proposed method	5.16E-15	3.73E-15

Table 3. Table 3: Singular case 2: Pose median errors.

	position	orientation
Kneip’s method	8.10E-14	8.85E-14
Masselli’s method	7.24E-14	6.07E-14
Proposed method	6.73E-14	1.75E-14

Equations

$${}^G{\mathbf{p}}_i$$

$${(}^G{\mathbf{p}}_1-{}^G{\mathbf{p}}_2{)}^T{}_C^G\mathbf{C}({}^C{\mathbf{b}}_1\times {}^C{\mathbf{b}}_2)$$

$${(}^G{\mathbf{p}}_1-{}^G{\mathbf{p}}_3{)}^T{}_C^G\mathbf{C}({}^C{\mathbf{b}}_1\times {}^C{\mathbf{b}}_3)$$

$${(}^G{\mathbf{p}}_2-{}^G{\mathbf{p}}_3{)}^T{}_C^G\mathbf{C}({}^C{\mathbf{b}}_2\times {}^C{\mathbf{b}}_3)$$

$${}_C^G\mathbf{C}=\mathbf{C}({\mathbf{k}}_1,{\theta }_1)\mathbf{C}({\mathbf{k}}_2,{\theta }_2)\mathbf{C}({\mathbf{k}}_3,{\theta }_3)$$

$${\mathbf{k}}_1\triangleq \frac{{}^G{\mathbf{p}}_1-{}^G{\mathbf{p}}_2}{\Vert {}^G{\mathbf{p}}_1-{}^G{\mathbf{p}}_2\Vert },{\mathbf{k}}_3\triangleq \frac{{}^C{\mathbf{b}}_1\times {}^C{\mathbf{b}}_2}{\Vert {}^C{\mathbf{b}}_1\times {}^C{\mathbf{b}}_2\Vert },{\mathbf{k}}_2\triangleq \frac{{\mathbf{k}}_1\times {\mathbf{k}}_3}{\Vert {\mathbf{k}}_1\times {\mathbf{k}}_3\Vert }$$

$${\mathbf{k}}_1^T\mathbf{C}({\mathbf{k}}_2,{\theta }_2){\mathbf{k}}_3$$

$$\mathbf{C}({\mathbf{k}}_2,{\theta }_2)=\cos {\theta }_2\mathbf{I}-\sin {\theta }_2\lfloor {\mathbf{k}}_2\rfloor +(1-\cos {\theta }_2){\mathbf{k}}_2{\mathbf{k}}_2^T$$

$${\theta }_2=\arccos ({\mathbf{k}}_1^T{\mathbf{k}}_3)\pm \frac{\pi }{2}$$

$${\mathbf{u}}_i^T\mathbf{C}({\mathbf{k}}_1,{\theta }_1)\mathbf{C}({\mathbf{k}}_2,{\theta }_2)\mathbf{C}({\mathbf{k}}_3,{\theta }_3){\mathbf{v}}_i=0$$

$${\mathbf{u}}_i\triangleq {}^G{\mathbf{p}}_i-{}^G{\mathbf{p}}_3,{\mathbf{v}}_i$$

$$\mathbf{C}({\mathbf{k}}_1,{\theta }_1)\mathbf{C}({\mathbf{k}}_2,{\theta }_2){\mathbf{C}}^T({\mathbf{k}}_1,{\theta }_1)=\mathbf{C}(\mathbf{C}({\mathbf{k}}_1,{\theta }_1){\mathbf{k}}_2,{\theta }_2)$$

$${\mathbf{u}}_i^T\mathbf{C}({\mathbf{k}}_1,{\theta }_1)\mathbf{C}(\mathbf{C}({\mathbf{k}}_2,{\theta }_2){\mathbf{k}}_3,{\theta }_3)\mathbf{C}({\mathbf{k}}_2,{\theta }_2){\mathbf{v}}_i$$

$$\Rightarrow {\mathbf{u}}_i^T\mathbf{C}({\mathbf{k}}_1,{\theta }_1)\mathbf{C}({\mathbf{k}}_3^{\prime },{\theta }_3){\mathbf{v}}_i^{\prime }$$

$${\mathbf{v}}_i^{\prime }$$

$${\mathbf{k}}_3^{\prime }$$

$${\mathbf{k}}_1^{\prime }\triangleq {\mathbf{C}}^T({\mathbf{k}}_2,{\theta }_2){\mathbf{k}}_1={\mathbf{k}}_3\times {\mathbf{k}}_2$$

$$\mathbf{C}(\mathbf{k},\theta )\mathbf{v}$$

$$=\left[\begin{array}{cc}-\lfloor \mathbf{k}{\rfloor }^2\mathbf{v}& -\lfloor \mathbf{k}\rfloor \mathbf{v}\end{array}\right]\left[\begin{array}{c}\cos \theta \\ \sin \theta \end{array}\right]+({\mathbf{k}}^T\mathbf{v})\mathbf{k}$$

$${\left(\left[\begin{array}{cc}-\lfloor {\mathbf{k}}_1{\rfloor }^2{\mathbf{u}}_i& \lfloor {\mathbf{k}}_1\rfloor {\mathbf{u}}_i\end{array}\right]\left[\begin{array}{c}\cos {\theta }_1\\ \sin {\theta }_1\end{array}\right]+({\mathbf{k}}_1^T{\mathbf{u}}_i){\mathbf{k}}_1\right)}^T$$

$$\cdot$$

$${\left[\begin{array}{c}\cos {\theta }_1\\ \sin {\theta }_1\end{array}\right]}^T\left[\begin{array}{cc}{\mathbf{u}}_i^T\lfloor {\mathbf{k}}_1{\rfloor }^2\lfloor {\mathbf{k}}_3^{\prime }{\rfloor }^2{\mathbf{v}}_i^{\prime }& {\mathbf{u}}_i^T\lfloor {\mathbf{k}}_1{\rfloor }^2\lfloor {\mathbf{k}}_3^{\prime }\rfloor {\mathbf{v}}_i^{\prime }\\ {\mathbf{u}}_i^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }{\rfloor }^2{\mathbf{v}}_i^{\prime }& {\mathbf{u}}_i^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }\rfloor {\mathbf{v}}_i^{\prime }\end{array}\right]\left[\begin{array}{c}\cos {\theta }_3\\ \sin {\theta }_3\end{array}\right]$$

$$+$$

$$=$$

$${\mathbf{u}}_1^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }\rfloor$$

$$=({}^G{\mathbf{p}}_1-{}^G{\mathbf{p}}_3{)}^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }\rfloor$$

$$=({}^G{\mathbf{p}}_1-{}^G{\mathbf{p}}_2+{}^G{\mathbf{p}}_2-{}^G{\mathbf{p}}_3{)}^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }\rfloor$$

$$=({}^G{\mathbf{p}}_2-{}^G{\mathbf{p}}_3{)}^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }\rfloor$$

$$={\mathbf{u}}_2^T{\mathbf{k}}_3^{\prime }{\mathbf{k}}_1^T={\mathbf{u}}_2^T\lfloor {\mathbf{k}}_1\rfloor \lfloor {\mathbf{k}}_3^{\prime }\rfloor$$

$$0$$

$$={\mathbf{u}}_i^T\mathbf{C}({\mathbf{k}}_1,{\theta }_1)\mathbf{C}({\mathbf{k}}_1,-\varphi )\mathbf{C}({\mathbf{k}}_1,\varphi )\mathbf{C}({\mathbf{k}}_3^{\prime },{\theta }_3){\mathbf{v}}_i^{\prime }$$

$$={\mathbf{u}}_i^T\mathbf{C}({\mathbf{k}}_1,{\theta }_1-\varphi )\mathbf{C}(\mathbf{C}({\mathbf{k}}_1,\varphi ){\mathbf{k}}_3^{\prime },{\theta }_3)\mathbf{C}({\mathbf{k}}_1,\varphi ){\mathbf{v}}_i^{\prime }$$

$$={\mathbf{u}}_i^T\mathbf{C}({\mathbf{k}}_1,{\theta }_1^{\prime })\mathbf{C}({\mathbf{k}}_3^{\prime \prime },{\theta }_3){\mathbf{v}}_i^{\prime \prime }$$

$${\theta }_1^{\prime }\triangleq {\theta }_1-\varphi ,{\mathbf{v}}_i^{\prime \prime }\triangleq \mathbf{C}({\mathbf{k}}_1,\varphi ){\mathbf{v}}_i^{\prime },{\mathbf{k}}_3^{\prime \prime }\triangleq \mathbf{C}({\mathbf{k}}_1,\varphi ){\mathbf{k}}_3^{\prime }$$

$$0$$

$$={\mathbf{u}}_1^T(\cos \varphi \mathbf{I}-\sin \varphi \lfloor {\mathbf{k}}_1\rfloor +(1-\cos \varphi ){\mathbf{k}}_1{\mathbf{k}}_1^T){\mathbf{k}}_3^{\prime }$$

$$=\cos \varphi {\mathbf{u}}_1^T{\mathbf{k}}_3^{\prime }-\sin \varphi {\mathbf{u}}_1^T\lfloor {\mathbf{k}}_1\rfloor {\mathbf{k}}_3^{\prime }$$

$$=\cos \varphi {\mathbf{u}}_1^T{\mathbf{k}}_3^{\prime }-\sin \varphi {\mathbf{u}}_1^T{\mathbf{k}}_2$$

$$\Rightarrow \left[\begin{array}{cc}\cos \varphi & \sin \varphi \end{array}\right]$$

Full text

1 Introduction

The Perspective-n-Point (PnP) is the problem of determining the 3D position and orientation (pose) of a camera from observations of known point features. The PnP is typically formulated and solved linearly by employing lifting (e.g., [1]), or as a nonlinear least-squares problem minimized iteratively (e.g., [10]) or directly (e.g., [12]). The minimal case of the PnP (for n=3) is often used in practice, in conjunction with RANSAC, for removing outliers [5].

The first solution to the P3P problem was given by Grunert [9] in 1841. Since then, several methods have been introduced, some of which [9, 4, 19, 5, 17, 8] were reviewed and compared, in terms of numerical accuracy, by Haralick et al. [11]. Common to these algorithms is that they employ the law of cosines to formulate a system of three quadratic equations in the features’ distances from the camera. They differ, however, in the elimination process followed for arriving at a univariate polynomial. Later on, Quan and Lan [21] and more recently Gao et al. [6] employed the same formulation but instead used the Sylvester resultant [3] and Wu-Ritz’s zero-decomposition method [23], respectively, to solve the resulting system of equations, and, in the case of [6], to determine the number of real solutions. Regardless of the approach followed, once the feature’s distances have been computed, finding the camera’s orientation, expressed as a unit quaternion [13] or a rotation matrix [14], often requires computing the eigenvectors of a $4\times 4$ matrix (e.g., [21]) or performing singular value decomposition (SVD) of a $3\times 3$ matrix (e.g., [6]), respectively, both of which are time-consuming. Furthermore, numerical error propagation from the computed distances to the rotation matrix significantly reduces the accuracy of the computed pose estimates.

To the best of our knowledge, the first method111Nister and Stewenius [20] also follow a geometric approach for solving the generalized P3P resulting into an octic univariate polynomial whose odd monomials vanish for the case of the central P3P. that does not employ the law of cosines in its P3P problem formulation is that of Kneip et al. [15], and later on that of Masselli and Zell [18]. Specifically, [15] and [18] follow a geometric approach for avoiding computing the features’ distances and instead directly solve for the camera’s pose. In both cases, however, several intermediate terms (e.g., tangents and cotangents of certain angles) need to be computed, which negatively affect the speed and numerical precision of the resulting algorithms.

Similar to [15] and [18], our proposed approach does not require first computing the features’ distances. Differently though, in our derivation, we first eliminate the camera’s position and the features’ distances to result into a system of three equations involving only the camera’s orientation. Then, we follow an algebraic process for successively eliminating two of the unknown 3-dof and arriving into a quartic polynomial. Our algorithm (summarized in Alg. 1) requires fewer operations and involves simpler and numerically more stable expressions, as compared to either [15] or [18], and thus performs better in terms of efficiency, accuracy, and robustness. Specifically, the main advantages of our approach are:

• •

Our algorithm’s implementation takes about 40% of the time required by the current state of the art [15]. 222Although Masselli and Zell [18] claim that their algorithm runs faster than Kneip et al.’s [15], our results (see Section 3) show the opposite to be true (by a small margin). The reason we arrive at a different conclusion is that our simulation randomly generates a new geometric configuration for each run, while Masselli employs only one configuration during their entire simulation, in which they save time due to caching.

• •

Our method achieves better accuracy than [15, 18] under nominal conditions. Moreover, we are able to further improve the numerical precision by applying root polishing to the solutions of the quartic polynomial while remaining faster than [15, 18].

• •

Our algorithm is more robust than [15, 18] when considering close-to-singular configurations (the three points are almost collinear or very close to each other).

The remaining of this paper is structured as follows. Section 2 presents the definition of the P3P problem, as well as our derivations for estimating first the orientation and then the position of the camera. In Section 3, we assess the performance of our approach against  [15] and [18] in simulation for both nominal and singular configurations. Finally, we conclude our work in Section 4.

2 Problem Formulation and Solution

2.1 Problem Definition

Given the positions, ${}^{\scriptscriptstyle{G}}\mathbf{p}_{i}$, of three known features $f_{i},\ i=1,2,3$, with respect to a reference frame $\{G\}$, and the corresponding unit-vector, bearing measurements, ${}^{\scriptscriptstyle{C}}\mathbf{b}_{i}$, $i=1,2,3$, our objective is to estimate the position, ${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}$, and orientation, i.e., the rotation matrix ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$, of the camera $\{C\}$.

2.2 Solving for the orientation

From the geometry of the problem (see Fig. 1), we have (for $i=1,2,3$):

[TABLE]

where $d_{i}\triangleq\|^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}-{}^{\scriptscriptstyle{G}}\mathbf{p}_{i}\|$ is the distance between the camera and the feature $f_{i}$.

In order to eliminate the unknown camera position, ${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}$, and feature distance, $d_{i},\ i=1,2,3$, we subtract pairwise the three equations corresponding to (1) for $(i,j)=(1,2),\ (1,3)$ and $(2,3)$, and project them on the vector ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}({}^{\scriptscriptstyle{C}}\mathbf{b}_{i}\times{}^{\scriptscriptstyle{C}}\mathbf{b}_{j})$ to yield the following system of 3 equations in the unknown rotation ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$:

[TABLE]

Next, and in order to compute one of the 3 unknown degrees of rotational freedom, we introduce the following factorization of ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$:

[TABLE]

where333$\mathbf{C}(\mathbf{k},\theta)$ denotes the rotation matrix describing the rotation about the unit vector, $\mathbf{k}$, by an angle $\theta$. Note that in the ensuing derivations, all rotation angles are defined using the left-hand rule.

[TABLE]

Substituting (5) in (2), yields a scalar equation in the unknown $\theta_{2}$:

[TABLE]

which we solve by employing Rodrigues’ rotation formula [16]:444$\lfloor\mathbf{k}\rfloor$ denotes the $3\times 3$ skew-symmetric matrix corresponding to $\mathbf{k}$ such that $\lfloor\mathbf{k}\rfloor\mathbf{a}=\mathbf{k}\times\mathbf{a}$, $\forall\ \mathbf{k},\mathbf{a}\in\mathbb{R}^{3}$. Note also that if $\mathbf{k}$ is a unit vector, then $\lfloor\mathbf{k}\rfloor^{2}=\mathbf{k}\mathbf{k}^{\scriptscriptstyle{T}}-\mathbf{I}$, while for two vectors $\mathbf{a}$, $\mathbf{b}$, $\lfloor\mathbf{a}\rfloor\lfloor\mathbf{b}\rfloor=\mathbf{b}\mathbf{a}^{\scriptscriptstyle{T}}-(\mathbf{a}^{\scriptscriptstyle{T}}\mathbf{b})\mathbf{I}$. Lastly, it is easy to show that $\lfloor\lfloor\mathbf{a}\rfloor\mathbf{b}\rfloor=\mathbf{b}\mathbf{a}^{\scriptscriptstyle{T}}-\mathbf{a}\mathbf{b}^{\scriptscriptstyle{T}}.$

[TABLE]

to get

[TABLE]

Note that we only need to consider one of these two solutions [in our case, we select $\theta_{2}=\arccos(\mathbf{k}_{1}^{\scriptscriptstyle{T}}\mathbf{k}_{3})-\frac{\pi}{2}$; see Fig. 2], since the other one will result in the same ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ (see Appendix 5.2 for a formal proof).

In what follows, we describe the process for eliminating $\theta_{3}$ from (3) and (4), and eventually arriving into a quartic polynomial involving a trigonometric function of $\theta_{1}$. To do so, we once again substitute in (3) and (4) the factorization of ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ defined in $\eqref{eq:ccc}$ to get (for $i=1,2$):

[TABLE]

where

[TABLE]

and employ the following property of rotation matrices

[TABLE]

to rewrite (10) in a simpler form as

[TABLE]

where

[TABLE]

The last equality in (13) is geometrically depicted in Fig. 2 and algebraically derived in Appendix 5.1. Analogously, it is straightforward to show that

[TABLE]

Next, by employing Rodrigues’ rotation formula [see (8)], for expressing the product of a rotation matrix and a vector as a linear function of the unknown $\begin{bmatrix}\cos\theta&\sin\theta\end{bmatrix}^{\scriptscriptstyle{T}}$, i.e.,

[TABLE]

in (12) yields (for $i=1,2$):

[TABLE]

Expanding (15) and rearranging terms, yields (for $i=1,2$)

[TABLE]

Notice that the term $\mathbf{u}_{i}^{\scriptscriptstyle{T}}\lfloor\mathbf{k}_{1}\rfloor\lfloor\mathbf{k}^{\prime}_{3}\rfloor$ appears three times in (2.2), and

[TABLE]

This motivates to rewrite (12) as (for $i=1,2$):

[TABLE]

where

[TABLE]

To simplify the equation analogous to (2.2) that will result from (18) [instead of (2.2)], we seek to find a $\phi$ (not necessarily unique) such that $\mathbf{u}_{1}^{\scriptscriptstyle{T}}\mathbf{k}^{\prime\prime}_{3}=0$, and hence, $\mathbf{u}_{i}^{\scriptscriptstyle{T}}\lfloor\mathbf{k}_{1}\rfloor\lfloor\mathbf{k}^{\prime\prime}_{3}\rfloor=0$ [see (17)], i.e.,

[TABLE]

[TABLE]

where

[TABLE]

and thus [from (19) using (8)]

[TABLE]

Now, we can expand (18) using (14) to get an equation analogous to (2.2):

[TABLE]

Substituting $\mathbf{u}_{i}^{\scriptscriptstyle{T}}\lfloor\mathbf{k}_{1}\rfloor\lfloor\mathbf{k}^{\prime\prime}_{3}\rfloor=0$ [see (17)] in (24) and renaming terms, yields (for $i=1,2$):

[TABLE]

where555The simplified expressions for the following terms, shown after the second equality, require lengthy algebraic derivations which we omit due to space limitations.

[TABLE]

For $i=1,2$, (25) results into the following system:

[TABLE]

Note that since $\bar{f}_{11}\bar{f}_{14}+\bar{f}_{12}\bar{f}_{15}=0$, we can further simplify (26) by introducing $\theta_{3}^{\prime}$, where

[TABLE]

Replacing $\theta_{3}$ by $\theta_{3}^{\prime}$ in (26), we have

[TABLE]

where

[TABLE]

From (28), we have

[TABLE]

Computing the norm of both sides of (2.2), results in

[TABLE]

which is a 4th-order polynomial in $\cos\theta_{1}^{\prime}$ that can be compactly written as:

[TABLE]

with

[TABLE]

We compute the roots of (38) in closed form to find $\cos\theta_{1}^{\prime}$. Similarly to [15] and [18], we employ Ferrari’s method [2] to attain the resolvent cubic of (38), which is subsequently solved by Cardano’s formula [2]. Once the (up to) four real solutions of (38) have been determined, an optional step is to apply root polishing following Newton’s method, which improves accuracy for minimal increase in the processing cost (see Section 3.2). Regardless, for each solution of $\cos\theta_{1}^{\prime}$, we will have two possible solutions for $\sin\theta_{1}^{\prime}$, i.e.,

[TABLE]

which, in general, will result in two different solutions for ${}^{\scriptscriptstyle{C}}_{\scriptscriptstyle{G}}\mathbf{C}$. Note though that only one of them is valid if we use the fact that $d_{i}>0$ (see Appendix 5.3).

Next, for each pair of $(\cos\theta_{1}^{\prime},\sin\theta_{1}^{\prime})$, we compute $\cos\theta_{3}^{\prime}$ and $\sin\theta_{3}^{\prime}$ from (2.2), which can be written as

[TABLE]

Lastly, instead of first computing $\theta_{1}$ from (19) and $\theta_{3}$ from (27) to find ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ using (5), we hereafter describe a faster method for recovering ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$. Specifically, from (5), (12) and (18), we have

[TABLE]

Since $\mathbf{k}_{1}$ is perpendicular to $\mathbf{k}_{3}^{\prime\prime}$, we can construct a rotation matrix $\mathbf{\bar{C}}$ such that

[TABLE]

and hence

[TABLE]

where

[TABLE]

Substituting (54) in (53), we have

[TABLE]

where

[TABLE]

The advantages of (55) are: (i) The matrix product $\mathbf{C}(\mathbf{e}_{1},\theta_{1}^{\prime})\mathbf{C}(\mathbf{e}_{2},\theta_{3}^{\prime})$ can be computed analytically; (ii) $\mathbf{\bar{C}},\mathbf{\bar{\bar{C}}}$ are invariant to the (up to) four possible solutions and thus, we only need to construct them once.

2.3 Solving for the position

Substituting in (1) the expression for $d_{3}$ from (62) and rearranging terms, yields

[TABLE]

Note that we only use (1) for $i=3$ to compute ${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}$ from ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$. Alternatively, we could find the position using a least-squares approach based on (1) for $i=1,2,3$ (see Appendix 5.4), if we care more for accuracy than speed. Lastly, the proposed P3P solution is summarized in Alg. 1.

3 Experimental results

Our algorithm is implemented666Our code is submitted along with the paper as supplemental material. in C++ using the same linear algebra library, TooN [22], as [15]. We employ simulation data to test our code and compare it to the solutions of [15] and [18]. For each single P3P problem, we randomly generate three 3D landmarks, which are uniformly distributed in a $0.4\times 0.3\times 0.4$ cuboid centered around the origin. The position of the camera is ${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}=\mathbf{e}_{3}$, and its orientation is ${}^{\scriptscriptstyle{C}}_{\scriptscriptstyle{G}}\mathbf{C}=\mathbf{C}(\mathbf{e}_{1},\pi)$.

3.1 Numerical accuracy

We generate simulation data without adding any noise or rounding error to the bearing measurements, and run all three algorithms on 50,000 randomly-generated configurations to assess their numerical accuracy. Note that the position error is computed as the norm of the difference between the estimate and the ground truth of ${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}$. As for the orientation error, we compute the rotation matrix that transforms the estimated ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ to the true one, convert it to the equivalent axis-angle representation, and use the absolute value of the angle as the error. Since there are multiple solutions to a P3P problem, we compute the errors for all of them and pick the smallest one (i.e., the root closest to the true solution).

The distributions and the means of the position and orientation errors are depicted in Fig.s 3 - 4 and Table 1. As evident, we get similar results to those presented in [18] for Kneip et al.’s [15] and Masselli and Zell’s methods [18], while our approach outperforms both of them by two orders of magnitude in terms of accuracy. This can be attributed to the fact that our algorithm requires fewer operations and thus exhibits lower numerical-error propagation.

Furthermore, and as shown in the results of Table 1, we can further improve the numerical precision by applying root polishing. Typically, two iterations of Newton’s method [24] lead to significantly better results, especially for the orientation, while taking only 0.01 $\mu$s per iteration, or about 4% of the total processing time.

3.2 Processing cost

We use a test program that solves 100,000 randomly generated P3P problems and calculates the total execution time to evaluate the computational cost of the three algorithms considered. We run it on a 2.0 GHz$\times$4 Core laptop and the results show that our code takes 0.54 $\mu$s on average (0.52 $\mu$s without root polishing) while [15] and [18] take 1.3 $\mu$s and 1.5 $\mu$s, respectively. This corresponds to a 2.5$\times$ speed up (or 40% of the time of [15]). Note also, in contrast to what is reported in [18], Masselli’s method is actually slower than Kneip’s. As mentioned earlier, Masselli’s results in [18] are based on 1,000 runs of the same features’ configuration, and take advantage of data caching to outperform Kneip.

3.3 Robustness

There are two typical singular cases that lead to infinite solutions in the P3P problem:

• •

Singular case 1: The 3 landmarks are collinear.

• •

Singular case 2: Any two of the 3 bearing measurements coincide.

In practice, it is almost impossible for these conditions to hold exactly, but we may still have numerical issues when the geometric configuration is close to these cases. To test the robustness of the three algorithms considered, we generate simulation data corresponding to small perturbations (uniformly distributed within $[-0.05~{}~{}0.05]$) of the features’ positions when in singular configurations. The errors are defined as in Section 3.1, while we compute the medians of them to assess the robustness of the three methods. For fairness, we do not apply root polishing to our code here. According to the results shown in Fig.s 5 - 8 and Tables 2 - 3, our method achieves the best accuracy in these two close-to-singular cases. The reason is that we do not compute any quantities that may suffer from numerical issues, such as cotangent and tangent in [15] and [18], respectively.

4 Conclusion and Future Work

In this paper, we have introduced an algebraic approach for computing the solutions of the P3P problem in closed form. Similarly to [15] and [18], our algorithm does not solve for the distances first, and hence reduces numerical-error propagation. Differently though, it does not involve numerically-unstable functions (e.g., tangent, or cotangent) and has simpler expressions than the two recent alternative methods [15, 18], and thus it outperforms them in terms of speed, accuracy, and robustness to close-to-singular cases.

As part of our ongoing work, we are currently extending our approach to also address the case of the generalized (non-central camera) P3P [20].

5 Appendix

5.1 Proof of $\mathbf{k}^{\prime}_{3}=\mathbf{k}_{2}\times\mathbf{k}_{1}$

First, note that $\mathbf{k}_{2}\times\mathbf{k}_{1}$ is a unit vector since $\mathbf{k}_{2}$ is perpendicular to $\mathbf{k}_{1}$. Also, from (13) and (7) we have

[TABLE]

Then, we can prove $\mathbf{k}^{\prime}_{3}=\mathbf{k}_{2}\times\mathbf{k}_{1}$ by showing that their inner product is equal to 1:

[TABLE]

5.2 Equivalence between the two solutions of $\theta_{2}$

When solving for $\theta_{2}$ [see (9)], we have two possible solutions $\theta_{2}^{(1)}=\arccos(\mathbf{k}_{1}^{\scriptscriptstyle{T}}\mathbf{k}_{3})-\frac{\pi}{2}$ and $\theta_{2}^{(2)}=\theta_{2}^{(1)}+\pi$. Next, we will prove that using $\theta_{2}^{(2)}$ to find ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ is equivalent to using $\theta_{2}^{(1)}$. First, note that (see Fig. 2)

[TABLE]

Then, we can write $\mathbf{C}(\mathbf{k}_{2},\theta_{2}^{(2)})$ as

[TABLE]

If we use $\theta_{2}^{(2)}$ to find ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$,

[TABLE]

Note that ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ in (60) is of the same form as that in (5), so any solutions of ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ computed using $\theta_{2}^{(2)}$ will be found by using $\theta_{2}^{(1)}$. Thus, we do not need to consider any other solutions for $\theta_{1}$ and $\theta_{3}$ beyond the ones found for ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$.

5.3 Determining the sign of $\sin\theta_{1}^{\prime}$

From (51), we have two solutions for $\sin\theta_{1}^{\prime}$, and thus for $\theta_{1}^{\prime}$, with $\theta_{1}^{\prime(2)}=-\theta_{1}^{\prime}$. This will also result into two solutions for $\theta_{3}^{\prime}$ [see (52)] and, hence, two solutions for $\theta_{3}$: $\theta_{3}$ and $\theta_{3}^{(2)}=\theta_{3}+\pi$. Considering these two options, we get two distinct solutions for ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$ [see (53)]:

[TABLE]

Then, notice that

[TABLE]

If $\mathbf{C}_{1}=\mathbf{C}_{2}$, this would require

[TABLE]

which cannot be true, hence $\mathbf{C}_{1}$ and $\mathbf{C}_{2}$ cannot be equal. Thus, there are always two different solutions of ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$.

If, however, we use the fact that $d_{i}\ (i=1,2,3)$ is positive, we can determine the sign of $\sin\theta_{1}^{\prime}$, and choose the valid one among the two solutions of ${}^{\scriptscriptstyle{G}}_{\scriptscriptstyle{C}}\mathbf{C}$. Subtracting (1) pairwise for ($i=3$) from ($i=1$), we have

[TABLE]

Multiplying both sides of (5.3) with ${\mathbf{k}_{3}^{\prime\prime}}^{\scriptscriptstyle{T}}\mathbf{C}(\mathbf{k}_{1},-\theta_{1}^{\prime})$ from the left, yields

[TABLE]

Using the fact that $d_{3}>0$ and $\delta>0$, we select the sign of $\sin\theta_{1}^{\prime}$ to be the same as that of $\mathbf{k}_{3}^{\scriptscriptstyle{T}}{}^{\scriptscriptstyle{C}}\mathbf{b}_{3}$.

5.4 Least-squares solution for the position

${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}$ can also be solved following a least-squares approach, which is slower but more accurate than (56). Specifically, (1) can result in the following system:

[TABLE]

Then, we only need to apply QR decomposition [7] and backsolve for ${}^{\scriptscriptstyle{G}}\mathbf{p}_{\scriptscriptstyle{C}}$ (i.e., we do not need to compute $d_{i},\ i=1,2,3$).

5.5 Derivation of $\bar{f}_{ij}$

[TABLE]

5.6 Derivation of $f_{ij}$ and $\bar{\bar{\mathbf{C}}}$

First, note that

[TABLE]

Let $\psi\triangleq\theta_{3}-\theta_{3}^{\prime}$, and thus

[TABLE]

From (63) and (27), we get

[TABLE]

[TABLE]

Then, from (26) and (28), we derive the expressions of $f_{ij}$:

[TABLE]

Additionally, we can derive the expression of $\mathbf{\bar{\bar{C}}}$, which is defined in (55):

[TABLE]

5.7 Comparison with the P3P code in OpenCV

We also compared the performance of our code with that in OpenCV (based on [6]), using the same setup as Sec. 3. The error distributions are showed in Fig. 9 and Fig. 10. It is obvious that the code in OpenCV has lower numerical accuracy comparing to ours. Also, it takes around 3 $\mu$s on average to compute P3P once, which is much slower than ours (0.52 $\mu$s according to Sec. 3).

In conclusion, our code performs much better than the one in OpenCV.