<figure><img src="image_1.png" /></figure>

HAL
open science

# Complete Singularity Analysis for the Perspective-Four-Point Problem

Beatriz Pascual-Escudero, Abhilash Nayak, Sébastien Briot, Olivier Kermorgant, Philippe Martinet, Mohab Safey El Din, François Chaumette

► To cite this version:

Beatriz Pascual-Escudero, Abhilash Nayak, Sébastien Briot, Olivier Kermorgant, Philippe Martinet, et al.. Complete Singularity Analysis for the Perspective-Four-Point Problem. International Journal of Computer Vision, 2021, 129 (4), pp.1217-1237. 10.1007/s11263-020-01420-0 . hal-03070525

HAL Id: hal-03070525

https://hal.science/hal-03070525

Submitted on 15 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Complete Singularity Analysis for the Perspective-Four-Point Problem

Beatriz Pascual-Escudero · Abhilash Nayak* · Sébastien Briot · Olivier Kermorgant · Philippe Martinet · Mohab Safey El Din · François Chaumette

Received: date / Accepted: date

**Abstract** This paper is concerned with pose estimation and visual servoing from four points. We determine the configurations for which the corresponding Jacobian matrix becomes singular, leading to inaccurate and unstable results. Using an adequate representation and algebraic geometry, it is shown that, for any orientation between the camera and the object,

B. Pascual-Escudero
École Centrale Nantes, Laboratoire des Sciences du Numérique de Nantes (LS2N), UMR CNRS 6004, Nantes, France
E-mail: beatriz.pascual.escudero@gmail.com
ORCID Id: 0000-0001-5957-0882

*Corresponding author: A. Nayak
Centre national de la Recherche Scientifique, Laboratoire des Sciences du Numérique de Nantes (LS2N), UMR CNRS 6004, Nantes, France
E-mail: Abhilash.Nayak@ls2n.fr
ORCID Id: 0000-0001-6228-203X

S. Briot
Centre national de la Recherche Scientifique, Laboratoire des Sciences du Numérique de Nantes (LS2N), UMR CNRS 6004, Nantes, France
E-mail: Sebastien.Briot@ls2n.fr
ORCID Id: 0000-0002-0419-6042

O. Kermorgant
École Centrale Nantes, Laboratoire des Sciences du Numérique de Nantes (LS2N), UMR CNRS 6004, Nantes, France
E-mail: Olivier.Kermorgant@ls2n.fr
ORCID Id: 0000-0002-1265-4138

P. Martinet
Inria Sophia Antipolis, Sophia Antipolis, France
E-mail: Philippe.Martinet@inria.fr
ORCID Id: 0000-0001-5827-0431

M. Safey El Din
Sorbonne Université, CNRS, LIP6, Equipe PolSys, Paris, France
E-mail: Mohab.Safey@lip6.fr
ORCID Id: 0000-0001-9463-1257

F. Chaumette
Inria, Univ. Rennes, CNRS, IRISA
Rennes, France
E-mail: Francois.Chaumette@inria.fr
ORCID Id: 0000-0002-1238-4385

there are always two to six singular locations of the camera in the generic case where the points are not coplanar, corresponding to the intersection of four cylinders. The particular case where the four points are coplanar is also characterized. Furthermore, some realistic example configurations are considered to substantiate the theory and to demonstrate failure cases in pose estimation and image-based visual servoing when the camera approaches a singularity.

**Keywords** Pose estimation · Visual servoing · Singularity · PnP

## 1 Introduction

Pose estimation is a classical problem in computer vision with many applications in augmented reality (Marchand et al, 2016) and robotics, especially in visual servoing (Hutchinson et al, 1996). When image points serve as image measurements and are matched with their corresponding 3D points, this problem is refered as PnP (Perspective-from-n-Points) (Horaud et al, 1989; Gao et al, 2003). While many efficient methods are available for a long time (Lowe, 1991; DeMenthon and Davis, 1995; Lu et al, 2010; Lepetit et al, 2009; Kneip et al, 2014), considering the failure cases did not receive a lot of attention from the community.

A classical method for solving the PnP consists in minimizing the reprojection error (Lowe, 1991):

$$ \hat{\mathbf{p}} = \arg \min_{\mathbf{p}} \sum_{i=1}^{n} (\mathbf{x}_i - \mathbf{\Pi} \, ^{\mathrm{c}}\mathbf{M}_o(\mathbf{p}) \mathbf{X}_i)^2 \quad (1) $$

where $\mathbf{p}$ are the six pose parameters representing the change of frame $^{c}\mathbf{M}_o$ between the object frame and the camera frame, $\mathbf{x}_i$ are the coordinates of the $n$ image points measured, and $\mathbf{X}_i$ are the known coordinates of their corresponding 3D-points expressed in the object frame, which are, after the change of frames $^{c}\mathbf{M}_o$, projected in the image using the perspective
projection matrix **Π**. The minimization of this non linear reprojection error is usually performed using iterative methods such as Levenberg-Marquardt or Newton-like descent, involving the Jacobian matrix $\frac{\partial x_i}{\partial p}$.

When the variation of the pose parameters are expressed in $sE_3$, the Lie algebra (tangent space at the identity) of the special Euclidean group $SE_3$, this Jacobian matrix relates the image point velocity to the relative camera-object velocity, leading to the basic *motion field* equations (Longuet-Higgins and Prazdny, 1980). In that case, the Jacobian matrix is equivalent to the *interaction matrix* at the basis of any image-based visual servoing scheme (Chaumette and Hutchinson, 2008). Using this form allows simple derivations and avoids any pose representation ambiguity (Marchand and Chaumette, 2002).

Determining when the interaction matrix becomes singular is crucial in order to avoid:

- accuracy issues when estimating the pose of the object. We will see that it occurs even if the pose is determined from methods that do not explicitly minimize (1), such as (DeMenthon and Davis, 1995; Lepetit et al, 2009; Kneip et al, 2014), (Persson and Nordberg, 2018).

- controllability issues in visual servoing due to the loss of rank of the interaction matrix. As a result, the singularities limit the workspace inside a convex space free of singularities.

These singularities can in theory be found by analyzing when the interaction matrix encounters a rank loss, but computations are highly heavy in practice. In the past, a geometric interpretation of the configurations leading to singularity cases for the P3P has been found in (Thompson, 1966). In (Michel and Rives, 1993), an appropriate decomposition of the interaction matrix allowed factorizing its determinant. Then, after heavy mathematical derivations, it was proven that singularities occur if the three 3-D points are aligned or if the optical center lies on the cylinder whose base is the circle that passes through the three points (see Fig. 1). The same result was also exhibited in (Papanikolopoulos, 1995) and more recently in (Tatsambon Fomena et al, 2011).

It is well known that four distinct solutions may exist for the P3P, which may lead the system to converge to a wrong pose in visual servoing. That is why, for both pose estimation and visual servoing, the common approach consists in observing additional features, i.e., $n \ge 4$ in (1). As far as we know, it was believed during many years that using at least four points in non degenerate configurations avoided singularities of the system. We will see that it was a wrong statement.

Let us also note that using other representations of the image points, such as the cylindrical coordinates instead of the Cartesian ones (Iwatsuki and Okiyama, 2005), the redundant spherical coordinates (Hamel and Mahony, 2002;

<figure><img src="image_2.png" /><figcaption>Fig. 1: Cylinder of singularities when observing 3 points.</figcaption></figure>

Burschka et al, 2003; Corke, 2010), or other angular representations (Liu et al, 2010), does not modify the problem. It is also the case for non linear combinations of the image points, such as moments (Tahri and Chaumette, 2005). Indeed, the interaction matrix of the new representation is nothing but the product of the transformation Jacobian with the interaction matrix of the original points. The singularities are thus the same whatever the representation used.

Recently, a tool named the “hidden robot concept” was introduced in (Briot and Martinet, 2013), in which it was proven to be efficient to determine:

- the singularity cases of vision-based controllers dedicated to parallel robots (Rosenzveig et al, 2014; Briot et al, 2015, 2016).

- the singularity cases of more general classes of visual servoing and 3-D localization problems related to the observation of *n* coplanar points (Briot et al, 2017a) and of three image lines (Briot et al, 2017b).

By using the hidden robot concept, it is possible to prove that the equations used in the vision-based controller or in the 3-D localization problem are equivalent to the kinematic equations of virtual robot architectures. From this correlation, it was possible to determine a simplified basis of the rows of the interaction matrix. Then, by using advanced tools such as the Grassmann-Cayley algebra (Ben-Horin and Shoham, 2006) and/or the Grassmann geometry (Merlet, 2006), we were able to find the singularity loci when *n* coplanar points were observed by a camera (Briot et al, 2017a). Our investigations showed that,

- for the observation of more than three coplanar points, singularities still occur, but only when all points share the same circle and the camera center is in this circle too, leading to the fact that the image points are all aligned.

- for the observation of more than three non planar points, a necessary condition of singularity is when all singularity cylinders associated with all subset of three points intersect in at least a single point. Our conclusion was
that, even for four observed points, this condition seemed unlikely to appear.

Our last assertion was questioned in (Hamel and Samson, 2017). It was shown in this work, through a rather complex mathematical proof, that singularities may exist in the observation of four non coplanar points. However, neither their relative configurations nor the singularity loci were given. The goal of the present paper is to fill this lack and to determine the singularity loci of the camera when observing four points.

The paper is organized as follows. Section 2 makes necessary recalls on the computation of the interaction matrix involved when observing n image points and explains how to compute a basis of the rows of this interaction matrix, basis which leads further to the simplification of the calculation of the singularity conditions. Section 3 gives complete results on the singularities when observing four points. Concepts from computational algebraic geometry are used to consider a simpler set of equations to derive the singular locations of the camera. Also, the particular case for which the four points are coplanar is discussed separately. Section 4 puts forth some realistic examples and shows the behavior of classical pose estimation methods and visual servoing when the camera approaches a singularity. Section 5 draws conclusions.

## 2 Computation of a basis for the rows of the interaction matrix when observing image points

In this section, we first make necessary recalls on:

- the computation of the interaction matrix corresponding to the observation of n image points,

- the computation of a basis for the rows of the interaction matrix based on Plücker coordinates, using a new interpretation of the results of (Briot et al, 2017a).

### 2.1 Interaction matrix related to the observation of points

We use the standard pin hole model with focal length equal to 1 and z-axis oriented along the optical axis. However, any other model based on projective geometry could be used (Michel and Rives, 1993).

A 3-D point $M_i$ with coordinates $[X_i Y_i Z_i]^T$ in the camera frame $\mathcal{F}_C$ is projected in the image plane on a 2-D point $m_i$ with coordinates $[x_i y_i]^T$ given by:

$$x_i = \frac{X_i}{Z_i}, \quad y_i = \frac{Y_i}{Z_i} \qquad (2)$$

The relation between the image point velocities $\dot{x}_i$ and $\dot{y}_i$ and the relative camera-object velocity $\tau_c$ expressed in the camera frame is obtained by differentiation of (2) and is

given by (Longuet-Higgins and Prazdny, 1980; Chaumette and Hutchinson, 2008)

$$\begin{bmatrix} \dot{x}_i \\ \dot{y}_i \end{bmatrix} = \mathbf{L}_i \tau_c$$

where

$$\mathbf{L}_i = \begin{bmatrix} -\frac{1}{Z_i} & 0 & \frac{x_i}{Z_i} & x_i y_i & -(1+x_i^2) & y_i \\ 0 & -\frac{1}{Z_i} & \frac{y_i}{Z_i} & 1+y_i^2 & -x_i y_i & -x_i \end{bmatrix} \quad (3)$$

Then, considering the observation of n points ($n \ge 3$) $M_1, M_2, \dots, M_n$, the interaction matrix linking the velocities of the points $m_i$ ($i = 1, 2, \dots, n$) grouped in the vector $\dot{\mathbf{s}} = [\dot{x}_1 \dot{y}_1 \dot{x}_2 \dot{y}_2 \dots \dot{x}_n \dot{y}_n]^T$ to the camera velocity $\tau_c$ by the relation

$$\dot{\mathbf{s}} = \mathbf{L}_{(n)} \tau_c \qquad (4)$$

is given by

$$\mathbf{L}_{(n)} = \begin{bmatrix} \mathbf{L}_1 \\ \mathbf{L}_2 \\ \vdots \\ \mathbf{L}_n \end{bmatrix} = \begin{bmatrix} -\frac{1}{Z_1} & 0 & \frac{x_1}{Z_1} & x_1 y_1 & -(1+x_1^2) & y_1 \\ 0 & -\frac{1}{Z_1} & \frac{y_1}{Z_1} & 1+y_1^2 & -x_1 y_1 & -x_1 \\ -\frac{1}{Z_2} & 0 & \frac{x_2}{Z_2} & x_2 y_2 & -(1+x_2^2) & y_2 \\ 0 & -\frac{1}{Z_2} & \frac{y_2}{Z_2} & 1+y_2^2 & -x_2 y_2 & -x_2 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ -\frac{1}{Z_n} & 0 & \frac{x_n}{Z_n} & x_n y_n & -(1+x_n^2) & y_n \\ 0 & -\frac{1}{Z_n} & \frac{y_n}{Z_n} & 1+y_n^2 & -x_n y_n & -x_n \end{bmatrix} \quad (5)$$

We recall that this matrix is equivalent to the Jacobian matrix $\frac{\partial s}{\partial p}$ involved in the minimization of (1).

Singularities appear when the matrix $\mathbf{L}_{(n)}$ is rank deficient, that is, when $\text{rank}(\mathbf{L}_{(n)}) < 6$. It amounts to the vanishing of all $(6 \times 6)$ minors of $\mathbf{L}_{(n)}$ and it turns out to be computationally expensive. In what follows, it will be shown that there is a computationally cheaper approach to determine these singularities.

### 2.2 Basis of the interaction matrix

In this section, we explain how to compute a basis for the rows of the interaction matrix associated with n points ($n \ge 1$). We show below that this basis can be represented by a system of Plücker lines. Let us first recall that the Plücker representation of a line $\mathcal{L}$ passing through a point M takes the following form:

$$\mathcal{L}: \begin{bmatrix} \mathbf{f} \\ \mathbf{m}(Q) \end{bmatrix} = \begin{bmatrix} \mathbf{f} \\ \overrightarrow{QM} \times \mathbf{f} \end{bmatrix} \qquad (6)$$

in which $\mathbf{f}$ is the direction of the line, and $\mathbf{m}(Q) = \overrightarrow{QM} \times \mathbf{f}$ is the moment of the line with respect to point Q. Now, let us analyze the expression of the interaction matrix $\mathbf{L}_i$ given in (3). It is obtained by first differentiating the perspective projection equations (2):

$$\begin{bmatrix} \dot{x}_i \\ \dot{y}_i \end{bmatrix} = \begin{bmatrix} (\dot{X}_i Z_i - X_i \dot{Z}_i)/Z_i^2 \\ (\dot{Y}_i Z_i - Y_i \dot{Z}_i)/Z_i^2 \end{bmatrix} = \frac{1}{Z_i^2} \begin{bmatrix} Z_i & 0 & -X_i \\ 0 & Z_i & -Y_i \end{bmatrix} \begin{bmatrix} \dot{X}_i \\ \dot{Y}_i \\ \dot{Z}_i \end{bmatrix}, \quad (7)$$
and then by using the kinematics equation

$$ \begin{bmatrix} \dot{X}_i \\ \dot{Y}_i \\ \dot{Z}_i \end{bmatrix} = -v_C - \omega_C \times \overrightarrow{CM_i} \qquad (8) $$

where $v_C$ and $\omega_C$ are respectively the translational and rotational components of the camera velocity $\tau_C = [v_C^T \omega_C^T]^T$. This relation can be rewritten under matrix form:

$$ \begin{bmatrix} \dot{X}_i \\ \dot{Y}_i \\ \dot{Z}_i \end{bmatrix} = [-\mathbf{I}_3 [\overrightarrow{CM_i}]_\times] \tau_C \qquad (9) $$

where $[\overrightarrow{CM_i}]_\times$ is the antisymmetric matrix related to vector $\overrightarrow{CM_i}$. By injecting (9) in (7), we immediately obtain

$$ \begin{bmatrix} \dot{x}_i \\ \dot{y}_i \\ \dot{z}_i \end{bmatrix} = \mathbf{L}_i \tau_C = \begin{bmatrix} \mathbf{f}_1^T (\overrightarrow{CM_i} \times \mathbf{f}_1)^T \\ \mathbf{f}_2^T (\overrightarrow{CM_i} \times \mathbf{f}_2)^T \end{bmatrix} \tau_C \qquad (10) $$

where $\mathbf{f}_1^T = \frac{1}{Z_i^2} [-Z_i \ 0 \ X_i]$ and $\mathbf{f}_2^T = \frac{1}{Z_i^2} [0 \ -Z_i \ Y_i]$. The two rows of the interaction matrix $\mathbf{L}_i$ can thus be interpreted as the Plücker representation of two straight lines passing through point $M_i$ and such that they are orthogonal to the plane whose normal is the line of sight $CM_i$ (note that $\overrightarrow{CM_i} = [X_i \ Y_i \ Z_i]$ is orthogonal to $\mathbf{f}_1$ and $\mathbf{f}_2$).

This result was extended in (Briot et al, 2017a) in which we proved that the basis $\xi_i$ of the rows of the matrix $\mathbf{L}_i$ can be written as:

$$ \xi_i = \begin{bmatrix} \xi_{i1} \\ \xi_{i2} \end{bmatrix} \qquad (11) $$

in which $\xi_{ij} = [\mathbf{f}_{ij}^T \mathbf{m}_{ij}(Q)^T] (j = 1, 2)$, where $\mathbf{f}_{ij}$ is any vector orthogonal to $\overrightarrow{CM_i}$, and $\mathbf{m}_{ij}(Q) = \overrightarrow{QM_i} \times \mathbf{f}_{ij}$, $Q$ being any arbitrary point. $\xi_{ij}$ is the Plücker representation of any line $\mathcal{L}_{ij}$ of direction $\mathbf{f}_{ij}$ contained in the plane $\mathcal{P}_i$, which is the plane passing through $M_i$ and whose normal is parallel to $\overrightarrow{CM_i}$ (see Fig. 2). A necessary and a sufficient condition for the basis $\xi_i$ to be of full rank is that $\mathbf{f}_{i1}$ should not be parallel to $\mathbf{f}_{i2}$.

The fact that the basis $\xi_i$ for the rows of the matrix $\mathbf{L}_i$ is spanned by a system of two lines allows mentioning the following remarks/results:

- Even if the matrix $\mathbf{L}_i$ at (10) is composed of a set of two Plücker lines whose expressions are given with respect to point C, degeneracy of any system of lines is independent of the location of the points where the lines are expressed (Merlet, 2006). However, when providing an analytical expression of the basis $\xi_i$, all lines expressions should be given at the same point.

- Degeneracy conditions of systems of lines are only dependent on the relative configuration of the lines between them. Therefore, the degeneracy conditions appearing in the observation of n points are independent of the relative orientation between the camera and object frames, as already shown in (Briot et al, 2017a).

<figure><img src="image_3.png" /><figcaption>Fig. 2: Representation of the planes $\mathcal{P}_1, \mathcal{P}_2, \mathcal{P}_3$ and $\mathcal{P}_{123}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3: Condition of singularity for the observation of three points: the four planes planes $\mathcal{P}_1, \mathcal{P}_2, \mathcal{P}_3$ and $\mathcal{P}_{123}$ intersect in a point O.</figcaption></figure>

- There are infinitely many lines $\mathcal{L}_{ij}$ spanning the basis $\xi_i$. However, we can select pairs of lines that simplify the final expression of the rows' basis of the interaction matrix, as it is shown below.

In (Briot et al, 2017a), we also showed that the singularity for the observation of three non-aligned points $M_1, M_2$ and $M_3$ occurs if and only if the planes $\mathcal{P}_1, \mathcal{P}_2, \mathcal{P}_3$ and the plane $\mathcal{P}_{123}$ containing the three points intersect in a point, that we denote by O (see Fig. 3). If the three points are aligned, $\mathcal{P}_{123}$ is not defined anymore, nonetheless, the planes $\mathcal{P}_1, \mathcal{P}_2$ and $\mathcal{P}_3$ still intersect in a point at infinity. In the general case, it
<figure><img src="image_5.png" /><figcaption>Fig. 4: Orthogonal projection over the plane $P_{123}$ of the four points, of the camera center $C$, and point $O$ at one possible configuration in which the camera is in the singularity cylinder of $M_1, M_2, M_3$.</figcaption></figure>

is interesting to note that point $O$ was proven to be on the circle $\mathcal{R}$ passing through $M_1, M_2$ and $M_3$. Moreover, it is contained on the diameter of $\mathcal{R}$ containing $C'$, the projection of the camera center $C$ on the plane $P_{123}$, on the extremity opposed to $C'$ (see Fig. 4). As a result, in order to simplify the expression of the basis $\xi_{(3)}$ of the rows of the interaction matrix $L_{(3)}$, as proposed in (Briot et al, 2017a), it is possible to take as a basis $\xi_i$ for the observation of the point $M_i$ ($i = 1, 2, 3$) the two lines $L_{ij}$ defined by:

$$ \mathbf{f}_{i1} = \mathbf{n}_{123} \times \overrightarrow{CM_i}, \mathbf{f}_{i2} = \overrightarrow{CM_i} \times \mathbf{f}_{i1}, \quad (12) $$

$$ \mathbf{m}_{i1}(O) = \overrightarrow{OM_i} \times \mathbf{f}_{i1}, \mathbf{m}_{i2}(O) = \overrightarrow{OM_i} \times \mathbf{f}_{i2}, \quad (13) $$

where $\mathbf{n}_{123}$ is the normal vector to plane $P_{123}$. By doing this, the line $L_{i1}$ contains both $M_i$ and $O$, ensuring $\mathbf{m}_{i1}(O) = \mathbf{0}$. As a result, in singularity, the basis $\xi_{(3)}$, when computing the lines expressions at point $O$, takes the following block-triangular form (after reordering the columns):

$$ \xi_{(3)} = \begin{bmatrix} \xi_{11} \\ \xi_{21} \\ \xi_{31} \\ \xi_{12} \\ \xi_{22} \\ \xi_{32} \end{bmatrix} = \begin{bmatrix} \mathbf{f}_{11}^T & \mathbf{0} \\ \mathbf{f}_{21}^T & \mathbf{0} \\ \mathbf{f}_{31}^T & \mathbf{0} \\ \mathbf{f}_{12}^T & \mathbf{m}_{12}^T \\ \mathbf{f}_{22}^T & \mathbf{m}_{22}^T \\ \mathbf{f}_{32}^T & \mathbf{m}_{32}^T \end{bmatrix} \quad (14) $$

Now, let us consider the interaction matrix $L_{(4)}$ associated with the observation of four points $M_1, M_2, M_3$ and $M_4$. As mentioned previously, a necessary condition for singularity is that all singularity cylinders associated with all subset of three points among the four observed ones must intersect in at least a single point. This means that the interaction matrix $L_{(3)}$ associated with the observation of the points $M_1, M_2$, and $M_3$ must be singular, and thus that point $O$ still exists. As a result, the row basis $\xi_{(3)}$ of $L_{(3)}$ still takes the form of Eq. (14). Moreover, the basis for the rows of the matrix $L_{(4)}$, defined by:

$$ \xi_{(4)} = \begin{bmatrix} \xi_{11}^T & \xi_{21}^T & \xi_{31}^T & \xi_{12}^T & \xi_{22}^T & \xi_{32}^T & \xi_{41}^T & \xi_{42}^T \end{bmatrix}^T \quad (15) $$

in which the expressions of $\xi_{i1}$ and $\xi_{i2}$ are provided at (12) and (13), takes the following form:

$$ \xi_{(4)} = \begin{bmatrix} \mathbf{f}_{11}^T & \mathbf{0} \\ \mathbf{f}_{21}^T & \mathbf{0} \\ \mathbf{f}_{31}^T & \mathbf{0} \\ \mathbf{f}_{12}^T & \mathbf{m}_{12}^T \\ \mathbf{f}_{22}^T & \mathbf{m}_{22}^T \\ \mathbf{f}_{32}^T & \mathbf{m}_{32}^T \\ \mathbf{f}_{41}^T & \mathbf{m}_{41}^T \\ \mathbf{f}_{42}^T & \mathbf{m}_{42}^T \end{bmatrix} \quad (16) $$

It should be mentioned that, when $\overrightarrow{CM_i} \perp P_{123}$ ($i = 1, \dots, 4$), then $\mathbf{f}_{i1} = \mathbf{0}$ and $\xi_{i1} = \mathbf{0}$ according to the choice we have made in (12). Thus, our representation of the basis degenerates. The way to deal with this issue is discussed below.

### 3 Singularities of the P4P problem

In this section, we first describe our parametrization of the system. Then, we analyze the singularities for P4P, i.e., the singularities of the matrix $\xi_{(4)}$.

#### 3.1 Parametrization of the configuration

Consider the fixed object frame $\mathcal{F}_P : (P, \mathbf{x}_P, \mathbf{y}_P, \mathbf{z}_P)$ where the vectors $\mathbf{x}_P, \mathbf{y}_P$ and $\mathbf{z}_P$ form an orthonormal right-handed basis. Without loss of generality, the first observed point $M_1$ is considered as the origin $P$, $M_2$ is placed on the $\mathbf{x}_P$ axis, $M_3$ in the $\mathbf{x}_P\mathbf{y}_P$ plane while $M_4$ is arbitrary. As long as $M_1, M_2$ and $M_3$ are not aligned, the plane $P_{123}$ is the same as the plane $\mathbf{z}_P = \mathbf{0}$. Consequently, in frame $\mathcal{F}_P$, the four observed points have the following coordinates (see Fig. 5):

$$ \begin{aligned} M_1 &= [0 \ 0 \ 0]^T, \\ M_2 &= [1 \ 0 \ 0]^T, \\ M_3 &= [d_1 \ d_2 \ 0]^T, \\ M_4 &= [d_3 \ d_4 \ d_5]^T, \end{aligned} \quad (17) $$

where $d_1, d_2, d_3, d_4, d_5$ are real values that parametrize points $M_3$ and $M_4$. Note that the distances are scaled so that the distance between $M_1$ and $M_2$ is equal to 1. However, the results obtained in the paper can be easily extended to distances without scaling as it will be shown in Section 3.2.3.

We will always assume that no two points coincide, since that would be a degenerate configuration in which we would have less than four points in practice. This means that the
following five cases will be discarded along the whole paper:

$$
\begin{align*}
d_1 &= d_2 = 0 (M_1 \equiv M_3), \\
d_1 - 1 &= d_2 = 0 (M_2 \equiv M_3), \\
d_3 &= d_4 = d_5 = 0 (M_1 \equiv M_4), \\
d_3 - 1 &= d_4 = d_5 = 0 (M_2 \equiv M_4), \\
d_3 - d_1 &= d_4 - d_2 = d_5 = 0 (M_3 \equiv M_4).
\end{align*}
$$

We will consider separately the configuration where the four
points are coplanar (in Section 3.2.2), that is, when

$$
either \ d_5 = 0 \text{ or } d_2 = 0, \tag{18}
$$

from the general configuration (in Section 3.2.1).

Another special case arises when

$$
\overrightarrow{CM_i} \perp \mathcal{P}_{123} \text{ for some } i = 1, \dots, 4. \tag{19}
$$

Note that, from the definition of $\mathcal{F}_P$, the normal of plane $\mathcal{P}_{123}$ is nothing but the $z_P$-axis. As a consequence, geometrically, this means that the camera is placed so that $C$ is vertically aligned with the point $M_i$ (in Fig. 5 it would appear right under $M_i$). In this case, as already mentioned at the end of Section 2.2, $\xi_{i1} = 0$, hence we should make a different choice of the line direction $\mathbf{f}_{i1}$ in order for $\xi_{(4)}$ to span the whole line space. We could do that, but instead, we choose here to study this case using the following observation: if we have a configuration in which $\overrightarrow{CM_1}$ is perpendicular to the plane $\mathcal{P}_{123}$, after renumbering the points, for example considering $M'_1 = M_1, M'_2 = M_2, M'_3 = M_4, M'_4 = M_3$, now $\overrightarrow{CM'_1}$ is not perpendicular to the new plane $\mathcal{P}_{1'2'3'}$ unless the four points are coplanar (remember that the case where all points are coplanar will be studied afterwards in Section 3.2.2). In fact, when any $\overrightarrow{CM'_i}$, $i = 1 \dots 4$ is perpendicular to $\mathcal{P}_{123}$, it will not be perpendicular to the three remaining planes $\mathcal{P}_{124}, \mathcal{P}_{134}$ or $\mathcal{P}_{234}$. Therefore, permuting the points such that the plane $\mathcal{P}_{123}$ is replaced with one of those planes makes $\overrightarrow{CM_i}$ no longer perpendicular to the latter. This fact can be exploited so that the matrix $\xi_{(4)}$ given in (16) is always well-defined for the singularity analysis with respect to four non coplanar points whatever the camera position.

Let us thus express the position of the camera with respect
to the object frame $\mathcal{F}_P$ by

$$
C = [XY Z]^T \tag{20}
$$

We will determine the coordinates $(X, Y, Z)$ of the cam-
era center in $\mathcal{F}_P$ that correspond to singular configurations.
Note that they are valid whatever the camera orientation. For
expressing the coordinates $(X_i, Y_i, Z_i)$ of points $M_i$ in the cam-
era frame $\mathcal{F}_C$ where the perspective equations (2) apply, the
following change of frame might be necessary:

$$
\begin{bmatrix} X_i \\ Y_i \\ Z_i \end{bmatrix} = {}^c\mathbf{R}_P \mathbf{M}_i - \begin{bmatrix} X \\ Y \\ Z \end{bmatrix} \quad (21)
$$

<figure><img src="image_6.png" /><figcaption>Fig. 5: Observation of four points. The coordinates of the points in $\mathcal{F}_P$ are also marked.</figcaption></figure>

where $^c\mathbf{R}_P$ is the rotation matrix between frames $\mathcal{F}_P$ and $\mathcal{F}_C$
and $\mathbf{M}_i$ are the coordinates of point $M_i$ in $\mathcal{F}_P$ given in (17).

For our study, it does not make sense to consider the cases
in which the camera position coincides with one observed
point, so the following cases will be discarded:

$$
\begin{align*}
X &= Y = Z = 0 (C \equiv M_1), \\
X - 1 &= Y = Z = 0 (C \equiv M_2), \\
X - d_1 &= Y - d_2 = Z = 0 (C \equiv M_3), \\
X - d_3 &= Y - d_4 = Z - d_5 = 0 (C \equiv M_4).
\end{align*}
$$

Let us recall now that, in order for the camera to be in
a singular position with respect to the observation of four
points, in particular it must be in a singular position with
respect to the observation of $M_1, M_2$ and $M_3$. For this reason,
we consider that the camera is in the singularity cylinder of
$M_1, M_2$ and $M_3$, assuming they are non-collinear (see (47)
in Appendix A and the supplementary Maple file for calcu-
lations). In such a case, as already mentioned, there exists a
point $O$ where the planes $\mathcal{P}_1, \mathcal{P}_2, \mathcal{P}_3$ and $\mathcal{P}_{123}$ intersect (Briot
et al, 2017a). Geometrically, this point $O$ is placed in the
circle containing $M_1, M_2$ and $M_3$, in the antipodal of the or-
thogonal projection $C'$ of the camera center $C$ over the plane
$\mathcal{P}_{123}$. Knowing the coordinates of points $M_1, M_2$ and $M_3$ in
the frame $\mathcal{F}_P$, the coordinates of point $O$ can be calculated in
the same frame. They are given by (see Fig. 4)

$$
O = \left[ 1 - X, \frac{d_1^2 + d_2^2 - d_1}{d_2} - Y, 0 \right]^T . \quad (22)
$$

In the next section, we determine the singularities of the
matrix $\xi_{(4)}$.
## 3.2 Singularities as intersection of four cylinders

In this section, we make the analysis that leads to a description of the singular configurations for the P4P, based on the basis $\xi_{(4)}$ provided in (16). In a first step, we consider the general case where the four points are not coplanar and we give a system of equations characterizing the positions of the camera leading to a singular configuration.

Later in Section 3.2.2, we will show how to treat the case where all points are coplanar.

### 3.2.1 General case (non coplanar points)

From (18), we consider here that

$$d_2 \neq 0 \text{ and } d_5 \neq 0. \qquad (23)$$

A fixed choice of the four observed points (a relative configuration between them) is determined by a choice of the values of the parameters $d_1, \dots, d_5$. We aim at giving equations for the singular positions $[XY Z]$ of the camera depending on the parameters $d_1, \dots, d_5$.

The most direct way of obtaining these equations is by computing the minors of the matrix $\xi_{(4)}$: whenever the 28 minors of size $(6 \times 6)$ of the matrix $\xi_{(4)}$ vanish simultaneously, the configuration is singular. However, it is not possible to deal with the polynomials arising from the computation of these minors by hand: we obtain a system of 28 polynomials in 3 variables $(X, Y, Z)$ with coefficients in 5 parameters $(d_1, d_2, d_3, d_4, d_5)$.

One preliminary observation is that, in order to have a singular configuration for the camera with respect to the four points, it is necessary that the camera is on the singularity cylinder generated by each choice of three points out of the four. Hence, as mentioned before, **the singularity locus for the observation of four points is contained in the intersection of four cylinders**, each of which is the singularity cylinder corresponding to three out of the four points.

If we choose three out of the four cylinders, an affine real version of Bézout theorem (Bézout (1779)) prescribes<fnref n="1" /> that there will be at most eight single real 3D points in this intersection. There is still one cylinder left to intersect, hence the intersection will be a (possibly empty) subset of these (at most eight) points. Let us highlight that the four cylinders have a special relation: in particular, three of them determine the fourth. However, it is unclear in principle if being in the intersection of the four cylinders is a sufficient condition for singularity. Indeed, the equations of the four cylinders correspond to four very specific minors of the matrix $\xi_{(4)}$. It is clear that in order for the 28 minors to vanish, in particular these four should vanish. This means that satisfying the

equations of the four cylinders is a necessary condition for singularity, although it is in principle not clear whether the vanishing of these four minors will automatically imply the vanishing of the 24 minors left.

To check that, let us simplify the equations of the minors. They can be factorized (see the supplementary Maple file), and all of them contain at least two of the following four factors:

$$X^2 + Y^2, X^2 + Y^2 - 2X + 1, \qquad (24)$$

$$d_1^2 - 2d_1X + d_2^2 - 2d_2Y + X^2 + Y^2,$$

$$d_3^2 - 2d_3X + d_4^2 - 2d_4Y + X^2 + Y^2.$$

In fact, all of these factors can be rewritten as $(X - X_{M_i})^2 + (Y - Y_{M_i})^2$, where $X_{M_i}$ and $Y_{M_i}$ are the X and Y coordinates of the point $M_i$, $i = 1, \dots, 4$, respectively. They only vanish when the camera is vertically aligned with $M_i$. So, we discard these factors from all the minors and call the resulting polynomials $g_1, \dots, g_{28}$. Consequently, the resulting polynomials corresponding to the cylinders are denoted by $C_{123}$, $C_{124}$, $C_{134}$ and $C_{234}$. They are displayed in (47)-(50) in Appendix A, where $C_{ijk}$ is the cylinder corresponding to points $M_i, M_j$ and $M_k$.

The polynomials $g_1, \dots, g_{28}$ generate an ideal (see Appendix B and Cox et al (2007) for further reading), denoted by $I_{28} = \langle g_1, \dots, g_{28} \rangle$. It is essentially the set of all (linear and nonlinear) combinations of the constituent polynomials using coefficients from the set $\mathbb{R}(d_1, d_2, d_3, d_4, d_5)[X, Y, Z]$ of all polynomials in the variables $X, Y, Z$ whose coefficients are rational fractions depending on $d_1, \dots, d_5$. Similarly, let $I_4 = \langle C_{123}, C_{124}, C_{134}, C_{234} \rangle$ be the ideal generated by the polynomials corresponding to the four simplified cylinder equations.

For each choice $(d_1, d_2, d_3, d_4, d_5) \in \mathbb{R}^5$ of the design parameters, the values of camera coordinates $(X, Y, Z) \in \mathbb{R}^3$ for which all polynomials in $I_{28}$ vanish lead us to the singular positions of the camera. This set is the real variety (see Appendix B) in $\mathbb{R}^3$ defined by the ideal $I_{28}$, denoted by $V(I_{28})$. Similarly, we write $V(I_4)$ for the set of all solutions in $\mathbb{R}^3$ to all polynomials in the ideal $I_4$.

In general, when we write $V(I_4)$ and $V(I_{28})$ we will mean varieties in $\mathbb{R}^3$ as just explained. However, for certain tasks we will consider all $d_1, \dots, d_5$ as variables also. This will then be made explicit, and at these moments $V(I_4)$, $V(I_{28})$ will consequently denote the corresponding varieties in $\mathbb{R}^8$.

Because the generating set of $I_4$ is contained in the one of $I_{28}$, $I_4$ is the smaller ideal, i.e. $I_4 \subset I_{28}$, and it turns out to be computationally easier to handle. However, the variety<fnref n="2" /> $V(I_4)$ might contain more points than $V(I_{28})$ ($V(I_{28}) \subseteq V(I_4) \subset \mathbb{R}^8$, see Appendix B). The following ideal $I_d$ computed in Maple (see Appendix B and the supplementary Maple file)

<footnote n="1">¹ Under a certain generic condition, which in this case is reduced to the cylinders being different from each other, and holds in any case for $d_5 \neq 0$.</footnote>

<footnote n="2">² We are considering varieties in $\mathbb{R}[d_1, d_2, d_3, d_4, d_5, X, Y, Z]$ in this discussion.</footnote>
has the property that $V(I_d)$ is the smallest affine variety containing all points which are in $V(I_4)$ but not in $V(I_{28})$:

$$I_d = \langle d_5, \rangle \qquad (25)$$

$$-d_1^2 d_4 - d_2^2 d_4 + d_2 d_3^2 + d_2 d_4^2 + d_1 d_4 - d_2 d_3, \qquad (26)$$

$$d_1^2 Y + d_2^2 Y - d_2 X^2 - d_2 Y^2 - d_1 Y + d_2 X, \qquad (27)$$

$$d_3^2 Y + d_4^2 Y - d_4 X^2 - d_4 Y^2 - d_3 Y + d_4 X. \qquad (28)$$

Note that the above ideal contains the polynomial $d_5$, which means that the equation $d_5 = 0$ must be satisfied by any point in $V(I_d)$. Therefore, if $d_5 \neq 0$ as assumed in the beginning of this section, $V(I_d) = \emptyset$, implying $V(I_4) = V(I_{28})$. Consequently, the real solutions of the minors' system $g_1 = 0, \dots, g_{28} = 0$ are exactly the intersection points of the four cylinders and henceforth only $I_4$ needs to be considered for further analysis.

These four polynomials are still not small, but computationally speaking they are small enough. In order to determine the variety$^3 V(I_4)$, a nicer set of generators of the ideal $I_4$ (that is, a nicer set of equations for the same geometric object) can be found by calculating a Gröbner basis (see Appendix B and refer Cox et al (2007)). The computer algebra system Maple is used to do so (see the supplementary Maple file). In particular, by choosing the lexicographical ordering$^4$ with $Z <_{lex} Y <_{lex} X$, we obtain a Gröbner basis consisting of three polynomials $F_1, F_2, F_3$ in the variables $X, Y, Z$ with parameters $d_1, \dots, d_5$ in the coefficients. These polynomials have the following form:

$$F_1(X) = c_6 X^6 + c_5 X^5 + c_4 X^4 + c_3 X^3 + c_2 X^2 + c_1 X + c_0, \qquad (29)$$

$$F_2(X, Y) = c_y Y + b_2(X), \qquad (30)$$

$$F_3(X, Z) = c_z Z + b_3(X). \qquad (31)$$

where $c_z, c_y, c_j, j = 1, \dots, 6$ are polynomial expressions in the parameters $d_1, \dots, d_5$ with real coefficients, while $b_2, b_3$ are polynomials in $X$, whose coefficients are themselves polynomial expressions in $d_1, \dots, d_5$ with real coefficients. Let us note that $c_6 = 16d_2^4 d_5^4 \neq 0$ (remember (23)).

Therefore, with the help of Gröbner basis, an initially indiscernible set of four quadratic equations in three variables transforms into a structured set of three equations of which the first one is univariate.

Observe that for parameter values at which the product $c_6 c_y c_z$ does not vanish, one obtains a triangular system which parameterizes the $(Y, Z)$-coordinates w.r.t. the roots of a univariate polynomial depending on $X$. A natural question

arises here: for those parameter values at which $c_y c_z \neq 0$ (remember $c_6 \neq 0$), do the specializations of the three polynomials $F_1, F_2, F_3$ describe the solutions to our initial equations where we have specialized the parameters?

To understand this better, one can start with an analogy with solving a linear system depending on parameters: when one has a linear system $AX = b$ where the entries of $A, b$ depend on parameters and $X$ is a vector of unknowns, one can express the solutions with respect to parameters using Cramer's formulas. Here the entries of $X$ will be obtained as rational fractions depending on the parameters where the denominator is simply the determinant of $A$; hence these formulas are valid after specialization as long as this determinant is not zero. In our context, things are similar since computing Gröbner bases requires only a finite amount of additions, subtractions, multiplications and divisions depending on the input coefficients. There exists a subset of the parameter space, defined by a polynomial equation $q = 0$ such that specializing the obtained Gröbner basis computation outside this is valid. To define $q$, one could take for instance the product of all denominators that appear during the Gröbner basis computation but it is demonstrated in Appendix C that this is an overkill and it is actually enough to take $q = c_y c_z$.

As a consequence, outside the subset $\{q = 0\} \subset \mathbb{R}^5$ of the parameter space, the set of equations in (29-31) determines the singularities for the observation of four points.

Besides, even when $q = 0$, we can still use (29) to solve for $X$. Therefore, all the results in the following regarding bounds on the number of singular positions and their geometrical interpretation are still valid. Once $X$ is determined, the system $I_4$ containing the four cylinder equations $C_{ijk}$ in Appendix A can be used instead of the Gröbner basis to solve for $Y$ and $Z$.

After specialization of the parameters, as long as $q \neq 0$, the singular positions for the camera will be those points $(X_S, Y_S, Z_S) \in \mathbb{R}^3$ such that

$$F_1(X_S) = F_2(X_S, Y_S) = F_3(X_S, Z_S) = 0. \qquad (32)$$

Indeed, given a fixed set of four points (or equivalently, a fixed set of values for the parameters $d_1, \dots, d_5$), satisfying the assumptions in this section, it is possible to obtain the singular positions for the camera, by just substituting the parameter values in the coefficients $c_y, c_z, c_6, \dots, c_0, b_2$ and $b_3$ and solving the resulting system in $X, Y, Z$. Moreover, we can evaluate the number of singular positions using the univariate polynomial in the Gröbner basis.

A first observation to be done from the univariate polynomial $F_1$ in (29) comes from its degree: there are at most 6 real values of $X$ for which the polynomial $F_1$ vanishes. Moreover, since $F_2$ and $F_3$ are linear polynomials in $Y$ and $Z$ resp., for any choice of $d_1, \dots, d_5$ such that $d_2 \neq 0, d_5 \neq 0$, there are at most 6 singular positions of the camera when observing four non coplanar points.

<footnote>$^3$ Through the rest of this section we are only regarding $X, Y, Z$ as variables, so we are considering varieties in $\mathbb{R}^3$.</footnote>

<footnote>$^4$ Gröbner basis with lexicographical ordering renders the system triangular and is analogous to the triangular form used by Gaussian elimination algorithm to solve linear equations (see Appendix B for a formal definition).</footnote>
One could wonder whether there is a choice of the pa-
rameter values for which we can guarantee that the above
system has less than 6 real solutions, maybe 4, 2 or even no
solutions. Observing the polynomial $F_1$ shown in the supple-
mentary Maple file, some of the coefficients $c_j, j = 1, ..., 6$
are dependent, which allows eliminating some of them as
follows:

$$
\begin{equation}
\begin{split}
F_1 = c_6 X^6 & - 3 c_6 X^5 + (c_2 + c_1 + 3 c_6) X^4 \\
& - (2 c_1 + 2 c_2 + c_6) X^3 + c_2 X^2 + c_1 X + c_0.
\end{split}
\tag{33}
\end{equation}
$$

Since $c_6 = 16d_2^4 d_5^4 \neq 0$ (according to (23)), we can normalize
$F_1$, obtaining

$$
\begin{equation}
\begin{split}
\bar{F}_1 ={}& X^6 - 3X^5 + \frac{(c_2 + c_1 + 3c_6)}{c_6} X^4 - \frac{(2c_1 + 2c_2 + c_6)}{c_6} X^3 \\
& + \frac{c_2}{c_6} X^2 + \frac{c_1}{c_6} X + \frac{c_0}{c_6}.
\end{split}
\tag{34}
\end{equation}
$$

Consider now the variable change (known as the linear Tschirnhausen transformation) $X = t - \frac{c_5}{6c_6} = t + \frac{1}{2}$. This change on $\overline{F}_1$ is supposed to remove only the fifth degree term, but thanks to the symmetry of the problem, it eliminates all odd powered terms giving rise to the following simpler polynomial in $t$:

$$
\overline{F}_1(\overline{X}) = G_1(t) = t^6 + C_4 t^4 + C_2 t^2 + C_0. \quad (35)
$$

The latter is a cubic equation in $t^2 = T$:

$$
\tilde{G}_1(T) = T^3 + C_4 T^2 + C_2 T + C_0. \quad (36)
$$

It has to be noted that $Z <_{lex} Y <_{lex} X$ is the only variable
order for which the term $\frac{c_5}{6c_6} = \frac{1}{2}$ in Tschirnhausen transfor-
mation is a constant. This helps us to geometrically interpret
the solutions to $F_1 = 0$.

Let $T_1$ be a real positive root of $\tilde{G}_1(T) = 0$. Then, $t = \pm \sqrt{T_1}$ are roots of $G_1(t)$ and hence two roots of the polynomial $F_1(X)$ are $X_1 = \sqrt{T_1} + \frac{1}{2}$ and $X_2 = -\sqrt{T_1} + \frac{1}{2}$. It leads to $X_2 = 1 - X_1$. Furthermore, substituting $X_1$ and $X_2$ in (30) and (31) shows that the solutions to the cylinder equations always appear in pairs. If $(X_1, Y_1, Z_1)$ is a singular point, then its pair $(X_2, Y_2, Z_2)$ is given by:

$$
X_2 = 2X^* - X_1, \quad Y_2 = 2Y^* - Y_1, \quad Z_2 = 2Z^* - Z_1
$$

$$
\begin{align*}
& \text{with } X^* = \frac{1}{2}, \quad Y^* = \frac{d_1^2 + d_2^2 - d_1}{2d_2}, \\
& Z^* = - \frac{d_4(d_1^2 + d_2^2) - d_2(d_3^2 + d_4^2 + d_5^2) - d_1d_4 + d_2d_3}{2d_2d_5},
\end{align*}
$$

where $(X^*, Y^*, Z^*)$ is the center of the circum-sphere of the
tetrahedron $M_1M_2M_3M_4$. Thus, it can be deduced that ev-
ery pair of singular points is symmetric about the center of
the sphere circumscribing the four observed points. In other
words, if we know one singular point of the camera, its pair

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6: A pair of singular points being symmetrical w.r.t. the center of the sphere circumscribing the four observed points (on top), and orthogonal projections of all the six singular points onto the plane $\mathcal{P}_{123}$ (bottom). +</figcaption></figure>

can be found on the sphere containing this point and with
center $(X^*, Y^*, Z^*)$. Moreover, the projection of these solu-
tions into the plane $\mathcal{P}_{123}$ shown in Fig. 4 reveals that for a
singular point $C(X, Y, Z)$, the projected pair of $C'(X, Y, 0)$ is
the antipodal point $O$. This is depicted in Fig. 6 for a config-
uration with $d_1 = -2$, $d_2 = -3$, $d_3 = 2$, $d_4 = -1$ and $d_5 = 7$,
where the pairs of $C$ are represented as $C_1$ and $C_2$ while their
projections on $\mathcal{P}_{123}$ as $C'_1$ and $C'_2 \equiv O_1$, $O_1$ being the antipo-
dal point of $C_1$. In this particular case, there exist six singular
positions for the camera. Their orthogonal projection onto
the plane $\mathcal{P}_{123}$ is also shown on this figure.

The only case in which $F_1$ will have real roots is when
$\tilde{G}_1(T)$ has at least one positive real root. It turns out that
for any values of $d_1, \dots, d_5$, $\tilde{G}_1(0) = -4p^2 \le 0$, where $p$ is
a common factor of polynomials $c_y$ and $c_z$ that appears in
$q$ (see the supplementary Maple file). So $\tilde{G}_1(0) < 0$ when
$q \ne 0$. Also, $\tilde{G}_1$ takes only positive values as $T$ grows to
infinity, so the intermediate value theorem tells that there is
at least one positive root $T$ for $\tilde{G}_1$, and therefore at least two
real roots for $F_1$.
Meanwhile, if $p = 0$ and hence $\tilde{G}_1(0) = 0$, the cubic in (36) can be factorized as $\tilde{G}_1(T) = T(T^2 + C_4T + C_2)$ implying that $T = 0$ is always a solution. This leads to two singular positions of the camera (due to symmetry), both with $X = \frac{1}{2}$. Geometrically, the line connecting the projections of these solutions in $\mathcal{P}_{123}$ is parallel to the Y-axis (see Fig. 6). However, since $p = 0$ implies $q = 0$, this case is particular because (30) and (31) cannot be used anymore to deduce Y and Z. This issue can easily be solved by substituting $X = \frac{1}{2}$ in the first cylinder equation (47), which results in the quadratic equation

$$Y^2d_2 - (d_1^2 + d_2^2 - d_1)Y - \frac{d_2}{4} = 0$$

whose discriminant is a sum of squares. So it yields two values for Y as expected. The solutions for X and Y can then be used to determine Z analytically from the other cylinder equations (48), (49) or (50).

As a result, $\tilde{G}_1(T)$ has at least one non-negative root resulting in at least two real roots for $F_1$. Therefore, there are at least 2 singular positions of the camera when observing four non coplanar points. The enclosed Maple file allows determining all the singular positions for the camera given any configuration for the four points.

Having established that there are at least two and at most six singular positions of the camera depending on the number of roots of the cubic (36), there can be exactly two or six singular positions as long as its discriminant does not vanish. On the contrary, if the discriminant vanishes, the cubic has a repeated solution and hence we expect four singular positions. As we will see just below, this is not entirely true. The discriminant of the cubic $\tilde{G}_1(T)$ is given by

$$\Delta = -16777216 d_2^6 d_5^4 (d_1 - d_3)^2 (d_1^2 d_4 + d_2^2 d_4 - d_2 d_3^2 - d_2 d_4^2 - d_1 d_4 + d_2 d_3)^2 (d_1^2 d_4^2 + d_1^2 d_5^2 + d_2^2 d_4^2 - d_2 d_3^2 d_4 - d_2 d_4^3 - d_2 d_4 d_5^2 - d_1 d_4^2 - d_1 d_5^2 + d_2 d_3 d_4)^2 (d_1^{12} d_4^6 + 3 d_1^{12} d_4^4 d_5^2 + \dots 1600 \text{ terms } + d_2^6 d_3^6 + d_2^6 d_3^4 d_5^2). \quad (37)$$

The last factor of $\Delta$ is difficult to interpret geometrically, while the rest of the algebraic factors also appear in the product $c_6c_yc_z$ (see (29)-(31) and Appendix C). As already mentioned, even when $q = c_yc_z = 0$, the current Gröbner basis can be used to determine bounds on the number of singular positions. However, it cannot be used to determine the exact number of solutions. For instance, when $d_1 = d_3$, $c_y = c_z = 0$ and (29) has four solutions for X but the number of solutions for Y and Z cannot be determined using (30) and (31). In fact, calculating a new Gröbner basis with $d_1 = d_3$ yields a quartic polynomial in X, followed by a quadratic polynomial in Y, linear in X and finally a linear polynomial in Z and X, which implies that there still could be six singular positions.

### 3.2.2 Case of four coplanar points

We recall that the four points are coplanar iff either $d_5 = 0$ or $d_2 = 0$ (see (17) and (18)). It is obvious when $d_5 = 0$, and it also occurs when $d_2 = 0$ implying that the points $M_1$, $M_2$ and $M_3$ are aligned. In both these cases, the discussion in Section 3.2.1 is no longer valid. Moreover, when $d_2 = 0$, we cannot use anymore the form of $\xi_{(4)}$ given by (16) since $d_2$ appears in the denominator of the Y-coordinate of point O (see (22)) with respect to which $m_{ij}$ is defined in (13). This can be overcome by simply re-defining the Plücker representation $\xi_{ij}$ in (11) with respect to a different point (say $M_1$ as shown below) although the simple structure in (16) has to be compromised.

A new basis is defined with each column $\xi_{ij}$, $i = 1, 2, 3, 4$, $j = 1, 2$ such that $f_{ij}$ is chosen according to (12) while the moment $m_{ij}$ of the line is chosen as follows:

$$\mathbf{m}_{i1}(M_1) = \overline{M_1 M_i} \times \mathbf{f}_{i1}, \quad \mathbf{m}_{i2}(M_1) = \overline{M_1 M_i} \times \mathbf{f}_{i2}. \quad (38)$$

Following the same procedure as in Section 3.2, vanishing of all $(6 \times 6)$ minors of the modified basis should yield the singular configurations.

By substituting $d_5 = 0$, the minors can be further decomposed into different polynomial systems, each of them defining one smaller component of the whole solution set. After eliminating the components that have no solutions, we obtain the following sets of polynomials, each of which defines a meaningful component of the solution set (see the supplementary Maple file):

$$h_{a_1} = [d_4, d_2], \quad (39)$$

$$h_{a_2} = [Z, \begin{aligned} &X^2d_4 + Y^2d_4 - Yd_3^2 - Yd_4^2 - Xd_4 + Yd_3, \\ &X^2d_2 + Y^2d_2 - Yd_1^2 - Yd_2^2 - Xd_2 + Yd_1, \\ &- d_1^2d_4 - d_2^2d_4 + d_2d_3^2 + d_2d_4^2 + d_1d_4 - d_2d_3. \end{aligned}]$$

The polynomials in $h_{a_1}$ vanish when the four points are aligned (that is $d_2 = d_4 = d_5 = 0$). Then, we obtain the well-known and intuitive result that the interaction matrix is always singular in that case.

The polynomials in $h_{a_2}$ vanish when the point $M_4$ and the camera center lie on the circle formed by points $M_1$, $M_2$ and $M_3$, which is the result already given in Briot et al (2017a)). Thus, a singular configuration occurs when the four points and the camera share the same circle.

Similarly, substituting $d_2 = 0$, the minors can be further decomposed as follows (see the supplementary Maple file):

$$h_{a_3} = [d_4, d_5], \quad (40)$$

$$h_{a_4} = [X - d_3, d_4(d_4 - Y) + d_5(d_5 - Z)].$$

The variety of $h_{a_3}$ with $d_2 = 0$ is the same as that of $h_{a_1}$ with $d_5 = 0$ in (39). In these cases where the four points are aligned, any configuration of the camera is singular.
The polynomials in $h_{a4}$ vanish when the camera lies on a line defined by

$$ \begin{cases} X = d_3 \\ d_4(d_4 - Y) + d_5(d_5 - Z) = 0. \end{cases} \tag{41} $$

Geometrically, this line is orthogonal to the plane containing the four points and passes through the fourth point $M_4$. Thus, **if the first three points are aligned but not the fourth one, then there is a line of singularities.**

Finally, remembering the discussion about (19), we have to consider the case where the four points are coplanar and the camera is vertically aligned with one of the points, let us say $M_i$. In that case, we recall that $\mathbf{f}_{i1} = 0$ inducing $\xi_{i1} = 0$, in which case the form of $\xi_{(4)}$ given by (16) cannot be used anymore. Yet again, a different basis $\xi_i$ should be chosen to analyse this case. From the definition of the basis in (11), since the camera is vertically aligned with $M_i$, $\overrightarrow{CM_i}$ is along $z_P$-axis in $\mathcal{F}_P$ and the direction vectors of the lines $\xi_{ij}$ can be simply chosen as parallel to $x_P$ and $y_P$ axes as follows:

$$ \mathbf{f}_{i1} = [1, 0, 0]^T, \quad \mathbf{f}_{i2} = \overrightarrow{CM_i} \times \mathbf{f}_{i1} = [0, 1, 0]^T \tag{42} $$

The moments $m_{ij}$ of the corresponding lines can be calculated using (38) to avoid missing the case $d_2 = 0$. Then, the minors of the modified basis $\xi_{(4)}$ consisting of the aforementioned direction and moment vectors are determined. For each $i = 1, 2, 3, 4$, the corresponding conditions such that $CM_i$ is along the $z_P$ axis are substituted in the minors and are further decomposed to obtain the following non-trivial sets of polynomials (see the supplementary Maple file):

$$ CM_1 \perp \mathcal{P}_{123} : \{X = Y = 0\} \qquad (43) $$

$$ CM_2 \perp \mathcal{P}_{123} : \{X - 1 = Y = 0\} $$

$$ CM_3 \perp \mathcal{P}_{123} : \{X - d_1 = Y - d_2 = 0\} $$

$$ CM_4 \perp \mathcal{P}_{123} : \{X - d_3 = Y - d_4 = 0\} $$

When the camera is vertically aligned with $M_1$ (that is $X = Y = 0$), we obtain $h_{b_1}$ whose element vanishes when the vectors $\overrightarrow{M_2M_3}$ and $\overrightarrow{M_2M_4}$ are parallel, implying the collinearity of points $M_2, M_3$ and $M_4$. Similarly, the components defined by $h_{b_2}, h_{b_3}$ and $h_{b_4}$ correspond to the vertical alignment of the camera with $M_2, M_3$ and $M_4$, respectively while the other three points are collinear. It is noteworthy that when the camera is vertically aligned with $M_4$ (that is $X = d_3, Y = d_4$ while $Z$ is arbitrary), we obtain $d_2 = d_5 = 0$ which is coherent with the variety of $h_{a4}$ in (40) with $d_2 = 0$. Thus, generalizing the result obtained from $h_{a4}$ in (40), **when the camera is vertically aligned with one of the points while the other three are collinear, the singularities are located**

<table><caption>Table 1: Different cases of singularities when four points are observed.</caption><thead><tr><th>Cases</th><th>Singularity configurations</th></tr></thead><tbody><tr><td>Four points non-coplanar (d<sub>2</sub> ≠ 0 and d<sub>5</sub> ≠ 0).</td><td>At least 2 and at most 6 real singularities</td></tr><tr><td>All four points are coplanar, but not concyclic, and no three of them collinear.</td><td>No singularities.</td></tr><tr><td>All four points are concyclic.</td><td>C concyclic with the points.</td></tr><tr><td>Only three points M<sub>i</sub>, M<sub>j</sub>, M<sub>k</sub> are collinear but not M<sub>l</sub>, (i, j, k, l all different ∈ 1, ..., 4)</td><td>C on the line through M<sub>l</sub> orthogonal to the plane containing the four points.</td></tr><tr><td>All four points are collinear.</td><td>Always singular.</td></tr></tbody></table>

**on the line passing through the former point and normal to the plane containing all the points.**

All cases corresponding to the different singular configurations are summarized in Table 1.

### 3.2.3 Homogenization of the cylinder equations

The cylinder equations (47)-(50) were derived by considering the point $M_2 = [100]^T$ as shown in (17), implying that the vector $\overrightarrow{M_1M_2}$ was normalized to be of unit length. For practical purposes, let us now consider $M_2, M_3, M_4$ and the camera center $C$ as $[D_000]^T$, $[D_1D_20]^T$, $[D_3D_4D_5]^T$ and $[X_dY_dZ_d]^T$, respectively with $D_0 ≠ 0$. Thus, the new cylinder equations are homogenized forms of the ones derived in this paper (they are listed in the supplementary Maple file). For instance, the first cylinder equation (47) becomes:

$$ -X_d^2D_2 + X_dD_2D_0 - Y_d^2D_2 + (D_1^2 + D_2^2 - D_0D_1)Y_d = 0. $$

By substituting $D_i = d_iD_0, i = 1, ..., 5, X_d = XD_0$ and $Y_d = YD_0$ in the above equation, leads to:

$$ D_0^3(-X^2d_2 + Xd_2 - Y^2d_2 + (d_1^2 + d_2^2 - d_1)d_3)Y = 0, $$

which is exactly (47) since $D_0 ≠ 0$. Consequently, when point $M_2$ is chosen as $[D_0, 0, 0]^T$, all points should be scaled such that $d_i = \frac{D_i}{D_0}$, $i = 0, ..., 5$ to use the equations proposed in this paper. Then, the coordinates of the obtained camera center must be multiplied with the same scaling factor $D_0$ to obtain the singularities corresponding to the chosen points.

## 4 Simulation results

This section illustrates the impact of the exposed singularities on the behavior of classical pose estimation methods and basic image-based visual servoing. We first introduce the considered study cases.
<table><caption>Table 2: Examples of points configuration and singular position of the camera. All units are in meters. In practice, only Ex. 1 is used in this paper as all configurations lead to similar behaviors.</caption><thead><tr><th>Ex.</th><th>d<sub>1</sub></th><th>d<sub>2</sub></th><th>d<sub>3</sub></th><th>d<sub>4</sub></th><th>d<sub>5</sub></th><th>Singular point (X<sub>s</sub>, Y<sub>s</sub>, Z<sub>s</sub>)</th></tr></thead><tbody><tr><td>1</td><td>0.6</td><td>0.4</td><td>0.3</td><td>0.9</td><td>0.3</td><td>(1.009, -0.082, 2.504)</td></tr><tr><td>2</td><td>0.7</td><td>0.3</td><td>0.3</td><td>0.8</td><td>0.3</td><td>(1.019, -0.341, 2.546)</td></tr><tr><td>3</td><td>0.9</td><td>0.3</td><td>0.3</td><td>0.9</td><td>0.4</td><td>(0.94, -0.237, 2.599)</td></tr><tr><td>4</td><td>0.6</td><td>0.4</td><td>0.9</td><td>0.9</td><td>0.3</td><td>(0, 0, 2.7)</td></tr><tr><td>5</td><td>0.5</td><td>0.5</td><td>0.9</td><td>0.9</td><td>0.3</td><td>(0, 0, 2.7)</td></tr><tr><td>6</td><td>0.9</td><td>0.5</td><td>0.3</td><td>0.9</td><td>0.3</td><td>(0.96, -0.092, 2.734)</td></tr><tr><td>7</td><td>0.5</td><td>0.4</td><td>0.9</td><td>0.9</td><td>0.3</td><td>(-0.012, -0.120, 3.015)</td></tr><tr><td>8</td><td>0.8</td><td>0.3</td><td>0.3</td><td>0.9</td><td>0.3</td><td>(0.988, -0.275, 3.177)</td></tr><tr><td>9</td><td>0.5</td><td>0.3</td><td>0.9</td><td>0.9</td><td>0.3</td><td>(-0.061, -0.339, 3.447)</td></tr><tr><td>10</td><td>0.8</td><td>0.3</td><td>0.3</td><td>0.8</td><td>0.3</td><td>(0.977, -0.306, 2.648)</td></tr></tbody></table>

## 4.1 Singular configurations

From the results obtained in Section 3.2.1, ten singularity configurations have been selected according to two physically relevant constraints:

- the points should fit in a box of dimension $1 \text{ m} \times 1 \text{ m} \times 1 \text{ m}$;

- the position of the camera should be at a distance between 2.4 and 4 m from the object.

The parameters describing the configurations of the points and the corresponding position of the camera are listed in Table 2. In practice, all configurations lead to a similar behavior and only Ex. 1 will be illustrated in this section<fnref n="5" />.

As already stated, the singular configurations are valid whatever the camera orientation. However, in practice, the points have to be in front of the camera to be observable. In addition, they should be in the camera field of view so that they are visible in the image plane. That is why a simple algorithm has been designed for determining a suitable orientation for the camera ensuring the following constraints:

$$ \forall i \in [1, 4] : |x_i| < 0.4, |y_i| < 0.4, Z_i > 0 $$

where $x_i$ and $y_i$ are the image point coordinates obtained from (2) (the value 0.4 corresponds to a typical field of view around $30^\circ$), and $Z_i$ is the depth of the point expressed in the camera frame thanks to (21).

## 4.2 Pose estimation from four points

Four classical pose estimation algorithms have been considered. The first two are ancient methods while the last two represent state-of-the-art solvers:

- the non-planar DeMenthon's method (DeMenthon and Davis, 1995);

- the Efficient PnP (Lepetit et al, 2009);

- the Universal PnP (Kneip et al, 2014); and

- a P4P approach based on lambda-twist P3P (Persson and Nordberg, 2018). All four P3P's are computed and the solution with the least reprojection error for the four points is kept.

In addition, the reprojection error (1) of the output solutions can be locally optimized either from:

- the Virtual Visual Servoing (VVS) method (Marchand and Chaumette, 2002). This approach performs first-order local optimization and relies directly on the interaction matrix;

- Levenberg-Marquardt (LM) method (Lowe, 1991) that uses second-order optimization.

Even if EPnP and UPnP have been designed to reduce the algorithmic complexity for a large number of points, we have chosen to discuss their behavior for the minimal configuration $n = 4$ where it can be applied.

All simulations presented in the following have been performed using the ViSP software library (Marchand et al, 2005), in which DeMenthon, VVS and Levenberg-Marquardt methods are implemented. Lambda-twist (P4P) implementation relies on the open-source code provided by the authors (Persson and Nordberg, 2018). Finally, EPnP and UPnP implementations are from OpenGV (Kneip and Furgale, 2014).

### 4.2.1 Pose estimation along growing spheres

This section highlights the impact of the singularity when performing pose computation in its neighborhood. All methods have been tested along spheres of growing radii (from 1 cm to 30 cm), centered on the singularity given by Ex. 1 in Table 2. For each radius, 1000 poses have been considered along the sphere.

As a singular configuration should be very sensitive to noise, a Gaussian noise of standard deviation $\sigma = 10^{-5}$ on the normalized image point coordinates is considered. This small noise level was chosen to illustrate the numerical sensitivity in the considered workspace.

Two metrics are extracted from the simulations:

- the translation error, as the Euclidean distance between the ground truth and the estimated translation; and

- the rotation error, as the absolute value of the error angle

$$ \theta_e = \left| \arccos \left( \text{tr}(\mathbf{R}_P^c \widehat{\mathbf{R}_P}^T) - 1 \right) / 2 \right| \quad \text{where } \mathbf{R}_P \text{ and } \widehat{\mathbf{R}_P} \text{ are respectively the true and estimated rotation matrices.} $$

The results are displayed in Fig. 7 and show that all methods fail to compute an accurate estimate (less than 10 mm and 0.5 dg error) near the singular point. As the sphere grows, the percentage of coarse estimations decreases. Raw estimators show very different behaviors:

<footnote n="5">⁵ Interested reader can perform all simulations using the code available on https://github.com/oKermorgant/singularity_4points</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 7: Percentage of estimation error greater than 10 mm (left) and 0.5 dg (right) along growing spheres for Ex. 1, with an image noise of $\sigma = 10^{-5}$. When near to the singularity, the error is larger than 10 mm and 0.5 dg on the whole sphere for all methods. As the poses are computed farther to the singular point, almost all methods lead to better estimations. Some raw estimators (DeMenthon, P4P, EPnP) still show a significant percentage larger than 10 mm or 0.5 dg when far from the singularity. On the opposite, UPnP and both local optimization variants (VVS, LM) lead to very accurate estimations as soon as the pose is computed far enough from the singular point, i.e., around 10 cm.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8: (a) 3D configuration for Ex. 1 in Table 2. A small rose pattern (blue) is performed around the singularity (red), at around 2.5 m from the four observed points (blue/orange). (b) Pose estimation from Dementhon and its refinements. (c) Estimation from UPnP and its refinements. While the initial estimation methods (blue) lead to non-perfect measurement (systematic for DeMenthon, only around singularity for UPnP), local optimization can usually reduce the error to 0. Yet, first order optimization (orange) may still yield imperfect estimation, especially around the singular point. Second-order optimization (green) is almost always perfect, except for one "petal" of the rose when starting from DeMenthon's result, where a part of the ground truth is still visible in red.</figcaption></figure>

- DeMenthon (blue) and EPnP (orange) are sensitive to the noise even far from the singularity.

- P4P (green) being computed from one particular P3P triplet, it is particularly sensitive to noise. The impact of the singularity is clearly visible, as the sensitivity is much higher when close to the singular point.

- UPnP (red) shows the best results of the four, with very good accuracy as soon as the pose is computed far enough from the singularity.

Local optimization, whether it is first-order (VVS) or second-order (LM), allows reducing the pose estimation error as soon as it is done far enough from the singularity. The only

exception is when starting from the EPnP estimate, where local optimization still fails in around 15 % of the points. As it will be shown in Section 4.3.4, this is due to the singularity point creating a local minimum for optimization-based methods. As the EPnP estimate is quite far from the real value, a significant percentage of local optimization are trapped into such local minimum. These results show that even with a very small noise level, the singularity point has a strong impact on the accuracy of pose estimation around its neighborhood.

We now illustrate the local behavior of the pose computation around the singularity with perfect image measurements.
### 4.2.2 Pose estimation along a rose pattern trajectory

An open-loop trajectory was simulated around the singular point given by Ex 1. This trajectory is a rose pattern with seven petals, of equation:

$$ \begin{cases} X_c(t) = X_s + 0.1 \cos(7t) \cos(t) \\ Y_c(t) = Y_s + 0.1 \cos(7t) \sin(t) \\ Z_c(t) = Z_s \end{cases} \quad (44) $$

where $(X_c, Y_c, Z_c)$ are the coordinates of the camera center while $(X_s, Y_s, Z_s)$ are the coordinates of the singularity point. The sampling time for this simulation is $\delta t = 0.01$ s. The corresponding trajectory and the observed points are displayed in Fig. 8a.

Note that no noise has been introduced, either in the coordinates of the 3D points, or in their image coordinates. Any inaccuracy in the results presented afterwards is thus only due to the proximity to a singular configuration.

The estimated position of the camera along the rose is depicted in Fig. 8.b for DeMenthon and Fig. 8.c for UPnP. EPnP is not shown as the errors are so large that the figure is unreadable. As we will see, P4P leads to perfect estimations in this noiseless simulation, and the corresponding figure is also omitted. DeMenthon's method estimates the whole rose with some offset (Fig. 8b). Local refinement allows going back to the rose, except when near the singularity for VVS, and except for a whole “petal” for both VVS and LM. On the opposite, the initial estimation from UPnP is much less continuous (Fig. 8c). This time, local refinement easily recovers the rose pattern except when really near to the singularity for VVS.

The corresponding translation errors are depicted in Fig. 9 for DeMenthon, UPnP and P4P, together with their local refinement. Again, EPnP is not shown due to very bad results. Even without any noise, DeMenthon's method shows a systematic offset in the pose estimation. On the opposite, UPnP has an error only around the singularity point. P4P leads to perfect estimations and only numerical noise (less than $10^{-6}$ m) appears when passing through the exact singular point. Local optimization highlights a key-impact of the singularity even in the case of perfect measurements, that is, the presence of local minima in the neighborhood of the singularity. The initial estimate from DeMenthon's algorithm being quite inaccurate, the wrong “petal” appears around iteration 250 on the left plot. This occurs both for first and second order optimizations. The vicinity of the singular point makes first-order optimization more sensitive to local minima, as it can be seen on the orange curves for DeMenthon and UPnP's initial guesses. On the opposite, second-order optimization (green curves) can usually reduce the error to 0 even at the singular point.

As a conclusion, even with perfect measurements the singularity has an impact in estimating the pose with raw estimators, except for P4P. This method is indeed perfectly

suited for noiseless measurements and is surprisingly stable in the vicinity of the singularity. Local refinement on the other hand may lead to better estimation, but still suffers from potential local minima.

Finally, the impact of the distance to the singularity on the condition number of the interaction matrix is depicted in Fig. 10, and illustrates a numerical validation of Section 3.2.

The next section details the local behavior of visual servoing around a singular configuration, and explains some of the observations made for local optimization-based approaches.

## 4.3 Visual servoing near a singularity

Visual servoing consists of controlling the motion of a robot in closed loop from visual data. At each iteration, the control output is computed from the error between the current value of a set of visual features $\mathbf{s}(t)$ and their desired value $\mathbf{s}^*$. The most basic controller tries to ensure a decoupled exponential decrease of the error $\mathbf{s} - \mathbf{s}^*$:

$$ \dot{\mathbf{s}} = -\lambda(\mathbf{s} - \mathbf{s}^*) \quad (45) $$

where $\lambda$ is a positive scalar. The resulting control law has the following form (Chaumette and Hutchinson, 2008):

$$ \tau_c = -\lambda \hat{\mathbf{L}}_s^+ (\mathbf{s} - \mathbf{s}^*) \quad (46) $$

where $\tau_c$ is the camera velocity sent to the robot controller, $\hat{\mathbf{L}}_s$ an estimation of the interaction matrix of $\mathbf{s}$, and $\hat{\mathbf{L}}_s^+$ its Moore-Penrose pseudo-inverse. In our case, the feature vector $\mathbf{s}$ is composed of the Cartesian coordinates of four image points, so $\mathbf{L}_s = \mathbf{L}_{(4)}$ given in (5). The depths $Z_i$ of the points appear in $\mathbf{L}_{(4)}$. When performing actual visual servoing, these depths are not known and have to be estimated in order to apply (46). On the opposite, in Virtual Visual Servoing, the depths are available from the current estimation of the camera pose. In this section, we focus on the local behavior of the servoing or VVS optimization, hence noiseless measurements are considered and the true value of the Z-depths are used.

### 4.3.1 Simulation setup

From Ex. 1 in Table 2, six servoings have been simulated starting from various initial positions, while the desired position was always the same. They are defined relative to the coordinates of the singularity point, and are listed in Table 3.

For the desired and initial poses, the camera orientation is computed automatically from their 3D position in order to ensure a positive depth for all points and to center them in the image, following the strategy explained in Section 4.1. For the control law (46), $\hat{\mathbf{L}}_s$ is computed from the current position $(x_i, y_i)$ of the points in the image and from the true value of the depths $Z_i$.
<figure><img src="image_11.png" /><figcaption>Fig. 9: Translation error along the rose pattern for Ex. 1. Vertical lines show when the camera passes at the singular point. DeMenthon (left), UPnP (middle) and P4P (right). Without noise, P4P only suffers from negligible numerical instability at the exact singular configuration. UPnP shows a large error around the singular point. Local optimization is most of the time able to reduce this error to 0, with a clear advantage for second-order method (green). Yet, a local minimum is still present for one petal of the rose when the optimization is initialized with the poor DeMenthon's output (left plot at iteration 250).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10: Distance to the singularity and inverse condition number of $L_{(4)}$ while following the rose pattern. Ill-condition clearly appears when passing near the singularity point.</figcaption></figure>

<table><caption>Table 3: Desired and starting positions for the visual servoing simulations. Positions are relative to the singularity position. All units are in meters.</caption><thead><tr><th>Position</th><th>ΔX</th><th>ΔY</th><th>ΔZ</th><th>Note on starting position</th></tr></thead><tbody><tr><td>Desired</td><td>-0.35</td><td>-0.35</td><td>0</td><td>-</td></tr><tr><td>Start 1</td><td>0.35</td><td>0.35</td><td>0</td><td>Exactly opposed to desired</td></tr><tr><td>Start 2</td><td>0.02</td><td>0.02</td><td>0</td><td>Almost at singularity</td></tr><tr><td>Start 3</td><td>-0.2</td><td>0.6</td><td>0</td><td>Away from singularity</td></tr><tr><td>Start 4</td><td>-0.7</td><td>0.3</td><td>0</td><td>Away from singularity</td></tr><tr><td>Start 5</td><td>0.5</td><td>-0.05</td><td>0</td><td>Almost opposed to desired</td></tr><tr><td>Start 6</td><td>0.3</td><td>0.4</td><td>0</td><td>Very near to Start 1</td></tr></tbody></table>

Figure 11 shows the trajectories achieved by the camera depending on the initial poses. As it can be seen, very different behaviors are obtained. We now highlight the key aspects of these results.

### 4.3.2 Classical behavior from Starts 1, 4 & 5

The most classical behavior is obtained when starting from Start 4 (red), where the camera converges directly to its desired pose. Starts 1 (blue) and 5 (purple), which are opposed to the desired position with regards to the singularity point, also converge to the desired pose. In these two cases, the

camera does not follow a straight line at all but instead avoids the singularity point. It suggests that the singularity point is a local repeller for the controller (46). This explains the imperfect performances of pose estimation by VVS near the singularity point, as seen in Section 4.2.2.

### 4.3.3 Unstable control from Start 2

Start 2 (orange) is very near to the singularity point, and the control scheme still converges to the desired pose. However, the interaction matrix is poorly conditioned as shown in Fig. 12 (left). This implies very high velocity commands around the singularity, as shown on the same figure. In practice, this velocity would have a very bad impact on a real robot controller. For VVS, it implies an erratic update on the estimated pose.

### 4.3.4 Reaching a local minimum from Starts 3 & 6

Starts 3 (green) and 6 (brown) lead the camera to converge to a local minimum, displayed as a dashed camera in Fig. 11. The evolution of the visual error norm $||s - s^*||$ for all starting poses is shown in Fig. 13. While this error follows the expected exponential decrease (45) for Starts 1, 2, 4, and 5, it converges to around $10^{-3}$ for Starts 3 & 6. This confirms that the final pose reached for these cases is not another solution to $s = s^*$ but an actual local minimum of Eq. (1). This means that the camera has reached a particular pose where $s - s^*$ is in the null space of $\hat{L}_s^T$. This phenomenon explains why pose refinement through VVS and LM sometimes converges to a very bad estimation of the pose depending on its initial guess, as was shown in the results given in Fig. 7.

### 4.3.5 Energy landscape around the singular point

In order to highlight the local behavior of optimization approaches, Fig. 14 shows the resulting energy landscape $||s(X, Y) - s^*||$ for the horizontal slice at the singularity point ($Z = Z_s$),
<figure><img src="image_13.png" /><figcaption>Fig. 11: Visual servoing simulations from 6 different starting poses (colored cameras). Top view is centered on the singularity point from Ex. 1 in Table 2. Starting from Pose 3 & 6, the system converges to a local minimum, displayed as a dashed camera. Bottom view shows the trajectories of the four points in the image plane. As the automatic camera orientation strategy tries to align all image points, their positions in the image are quite close even if they correspond to large camera displacements. Local minima appear in the zoomed window, where the black diamond is the desired 2D position while the red dot is the final one for Starts 3 & 6.</figcaption></figure>

with the same orientation as the desired one. Yellow color indicates a high value, while dark blue means a small one. While these results are very partial (only 2D while the search space of the optimization is SE(3)), a local minimum is clearly visible next to the singularity. It explains the potential local minima that arise when refining a pose in this region. As the local minimum seems to be only on one side of the singularity, the optimization being able to converge to the true pose depends on the initial guess. As an example, most of DeMenthon's initial guess along the rose pattern in Sec-

<figure><img src="image_14.png" /><figcaption>Fig. 12: Inverse condition number of $L_{(4)}$ (left) during the first iterations for Start 2 (orange). Induced camera velocity during the servoing (right). The poor conditioning of $L_{(4)}$ leads to very high velocities when passing near the singularity point.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13: Evolution of the visual error norm $||s - s^*||$ (in log-scale). Starts 3 & 6 lead to converge to a local minimum, while the other initial poses provide a perfect exponential decrease according to (45).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14: Energy landscape on a 1 m horizontal square centered on the desired position with high values in yellow to low ones in dark blue. In the left picture the desired position is the singular point and reveals a local repeller next to the singular point. On the opposite, the right picture shows the same energy landscape when the desired position is far from any singularity. In this case, the landscape shows a single global minimum.</figcaption></figure>

tion 4.2.2 is on the good side of the global minimum, except for a petal of the rose.

## 5 Conclusions

Our work dealt with the complete determination of the singularity configurations when observing four points with a
camera. We were able to identify a basis of the rows of the corresponding Jacobian matrix, expressed as the interaction matrix (5). The determination of the singular configurations is then done from the 28 (6 × 6)-minors of this basis. This turns out to be a better approach in terms of analytical computations compared to considering directly the interaction matrix, and allowed us to obtain the equations for determining the singular configurations for each choice of 4 points. Using this approach, we were able to prove that, for any relative orientation between the object and camera frames, and for any generic choice of four points, singular positions for the camera center do exist, which is a new result. These singular positions are just a finite number of isolated points in 3D space, which can be computed as the intersection of four cylinders. More precisely, there are at least two and at most six singular positions. They have been obtained through Gröbner basis computations.

To verify our theoretical results, some realistic example configurations of the four points were considered. The simulations showed that pose estimation algorithms may fail, sometimes drastically, when the camera approaches a singularity. As for image-based visual servoing, when the camera moves near a singularity, the control scheme exhibits high velocities, which was expected. It is also sometimes unable to reach its desired location, but instead converges to a local minimum.

Since any generic configuration of four points is shown to have at least 2 and at most 6 singular locations of the camera, we expect that adding a fifth point in a generic configuration should have no singularities. However, there might be special locations of the fifth point leading to singularities. Thus, future work will focus on studying the potential singular configurations when observing five (or more) points.

# Appendices

## A. Cylinders Equations

Equations of the four cylinders $C_{ijk}$ containing points $M_i, M_j, M_k$:

$$C_{123} = -X^2 d_2 + X d_2 - Y^2 d_2 + (d_1^2 + d_2^2 - d_1) Y, \quad (47)$$

$$C_{124} = (-d_4^2 - d_5^2) X^2 + (d_4^2 + d_5^2) X - Y^2 d_4^2 - 2 Y Z d_4 d_5 + d_4 (d_3^2 + d_4^2 + d_5^2 - d_3) Y - Z^2 d_5^2 + d_5 (d_3^2 + d_4^2 + d_5^2 - d_3) Z, \quad (48)$$

$$C_{134} = (-d_1^2 d_4^2 - d_1^2 d_5^2 + 2 d_1 d_2 d_3 d_4 - d_2^2 d_3^2) X^2 - 2 d_1 d_2 d_5^2 XY + 2 d_2 d_5 (d_1 d_4 - d_2 d_3) ZX + (d_1^3 d_4^2 + d_1^3 d_5^2 - d_1^2 d_2 d_3 d_4 + d_1 d_2^2 d_4^2 + d_1 d_2^2 d_5^2 - d_1 d_2 d_3^2 d_4 - d_1 d_2 d_4^3 - d_1 d_2 d_4 d_5^2 - d_2^3 d_3 d_4 + d_2^2 d_3^3 + d_2^2 d_3 d_4^2 + d_2^2 d_3 d_5^2) X + (-d_1^2 d_4^2 + 2 d_1 d_2 d_3 d_4 - d_2^2 d_3^2 - d_2^2 d_5^2) Y^2 - 2 d_1 d_5 (d_1 d_4 - d_2 d_3) ZY + (-d_1^3 d_3 d_4 + d_1^2 d_2 d_3^2 + d_1^2 d_2 d_5^2 + d_1^2 d_3^2 d_4 + d_1^2 d_4^3 + d_1^2 d_4 d_5^2 - d_1 d_2^2 d_3 d_4 - d_1 d_2 d_3^3 - d_1 d_2 d_3 d_4^2 - d_1 d_2 d_3 d_5^2 + d_2^3 d_3^2 + d_2^3 d_5^2) Y - d_5^2 (d_1^2 + d_2^2) Z^2 - d_5 (d_1^2 + d_2^2) (d_1 d_3 + d_2 d_4 - d_3^2 - d_4^2 - d_5^2) Z, \quad (49)$$

$$C_{234} = -d_5 (d_1^3 d_3 + d_1^2 d_2 d_4 - d_1^2 d_3^2 - d_1^2 d_4^2 - d_1^2 d_5^2 + d_1 d_2^2 d_3 + d_2^3 d_4 - d_2^2 d_3^2 - d_2^2 d_4^2 - d_2^2 d_5^2 - d_1^3 - d_1^2 d_3 - d_1 d_2^2 + 2 d_1 d_3^2 + 2 d_1 d_4^2 + 2 d_1 d_5^2 - d_2^2 d_3 + 2 d_1^2 - d_1 d_3 + 2 d_2^2 - d_2 d_4 - d_3^2 - d_4^2 - d_5^2 - d_1 + d_3) Z + (-d_1^3 d_3 d_4 + d_1^2 d_2 d_3^2 + d_1^2 d_2 d_5^2 + d_1^2 d_3^2 d_4 + d_1^2 d_4^3 + d_1^2 d_4 d_5^2 - d_1 d_2^2 d_3 d_4 - d_1 d_2 d_3^3 - d_1 d_2 d_3 d_4^2 - d_1 d_2 d_3 d_5^2 + d_2^3 d_3^2 + d_2^3 d_5^2 + d_1^3 d_4 - 2 d_1^2 d_2 d_3 + d_1^2 d_3 d_4 + d_1 d_2^2 d_4 + d_1 d_2 d_5^2 + d_1 d_2^2 d_3^2 + d_1 d_2 d_3 d_4 + d_1 d_2 d_5^2 - 2 d_1 d_3^2 d_4 - 2 d_1 d_4^3 - 2 d_1 d_4 d_5^2 - 2 d_2^3 d_3 + d_2^2 d_3 d_4 + d_2 d_3^3 + d_2 d_3 d_4^2 + d_2 d_3 d_5^2 + d_1^2 d_2 - 2 d_1^2 d_4 + d_1 d_2 d_3 + d_1 d_3 d_4 + d_2^3 - d_2^2 d_4 - 2 d_2 d_3^2 - d_2 d_4^2 - 2 d_2 d_5^2 + d_2^3 d_4 + d_4^3 + d_4 d_5^2 - d_1 d_2 + d_1 d_4 + d_2 d_3 - d_3 d_4) Y - d_5^2 (d_1^2 + d_2^2 - 2 d_1 + 1) Z^2 + (d_1^3 d_4^2 + d_1^3 d_5^2 - d_1^2 d_2 d_3 d_4 + d_1 d_2^2 d_4^2 + d_1 d_2^2 d_5^2 - d_1 d_2 d_3^2 d_4 - d_1 d_2 d_4^3 - d_1 d_2 d_4 d_5^2 - d_2^3 d_3 d_4 + d_2^2 d_3^3 + d_2^2 d_3 d_4^2 + d_2^2 d_3 d_5^2 + d_1^2 d_2 d_4 - d_1^2 d_4^2 - d_1^2 d_5^2 + d_2^3 d_4 - d_2^2 d_3^2 - 2 d_2^2 d_4^2 - 2 d_2^2 d_5^2 + d_2 d_3^2 d_4 + d_2 d_4^3 + d_2 d_4 d_5^2 + d_1 d_2 d_4 - d_1 d_4^2 - d_1 d_5^2 - d_2^2 d_3 + d_2 d_3 d_4 + d_2^2 - 2 d_2 d_4 + d_4^2 + d_5^2) X + (-d_1^2 d_4^2 + 2 d_1 d_2 d_3 d_4 - d_2^2 d_3^2 - d_2^2 d_5^2 - 2 d_1 d_2 d_4 + 2 d_1 d_4^2 + 2 d_2^2 d_3 - 2 d_2 d_3 d_4 - d_2^2 + 2 d_2 d_4 - d_4^2) Y^2 - d_2 d_3^2 d_4 - d_2 d_4 d_5^2 + d_2 d_3 d_4 + (-d_1^2 d_4^2 - d_1^2 d_5^2 + 2 d_1 d_2 d_3 d_4 - d_2^2 d_3^2 - 2 d_1 d_2 d_4 + 2 d_1 d_4^2 + 2 d_1 d_5^2 + 2 d_2^2 d_3 - 2 d_2 d_3 d_4 - d_2^2 + 2 d_2 d_4 - d_4^2 - d_5^2) X^2 - 2 d_2 d_5^2 (d_1 - 1) XY + 2 d_2 d_5 (d_1 d_4 - d_2 d_3 + d_2 - d_4) ZX - 2 d_5 (d_1 - 1) (d_1 d_4 - d_2 d_3 + d_2 - d_4) ZY - d_2 d_4^3 - d_1 d_4^2 - d_1 d_5^2 - d_2^2 d_3 - d_2^3 d_4 + 2 d_2^2 d_4^2 - d_1^2 d_2 d_4 + d_1 d_2 d_4 - d_1 d_2^2 d_5^2 + d_2^3 d_3 d_4 + d_1 d_2 d_4^3 - d_1 d_2^2 d_4^2 - d_2^2 d_3 d_4^2 - d_2^2 d_3 d_5^2 - 2 d_1 d_2 d_3 d_4 + d_1^2 d_2 d_3 d_4 + d_1 d_2 d_3^2 d_4 + d_1 d_2 d_4 d_5^2 + 2 d_2^2 d_5^2 - d_1^3 d_4^2 - d_1^3 d_5^2 - d_2^2 d_3^3 + 2 d_1^2 d_4^2 + 2 d_1^2 d_5^2 + 2 d_2^2 d_3^2 \quad (50)$$

## B. Algebraic geometry tools

Most of the following definitions and examples are taken from Cox et al (2007).
## Ideal

A subset $\mathcal{I}$ of the ring of polynomials in the variables $x_1, \dots, x_n$ with coefficients in a field $k$, $k[x_1, \dots, x_n]$ is an **ideal** if it satisfies:

$$ (i) \quad 0 \in \mathcal{I} \\ (ii) \quad \text{If } f, g \in \mathcal{I} \text{ then } f + g \in \mathcal{I} \qquad (51) $$

$$ (iii) \quad \text{If } f \in \mathcal{I} \text{ and } h \in k[x_1, \dots, x_n], \text{ then } hf \in \mathcal{I}. $$

Let $f_1, \dots, f_s$ be polynomials in $k[x_1, \dots, x_n]$, then

$$ \langle f_1, \dots, f_s \rangle = \left\{ \sum_{i=1}^{s} h_i f_i : h_1, \dots, h_s \in k[x_1, \dots, x_n] \right\} $$

is an ideal generated by $f_1, \dots, f_s$. This set is interesting because, for any field $K$ containing $k$, all polynomials $f \in \langle f_1, \dots, f_s \rangle$ give polynomial equations $f = 0$ that hold at all points of $K^n$ at which $f_1 = 0, \dots, f_s = 0$ hold simultaneously, plus it has a structure that allows many useful algebraic operations.

## Variety

If $f_1, \dots, f_s$ are polynomials in the ring $k[x_1, \dots, x_n]$ and if

$$ V(f_1, \dots, f_s) = \{(a_1, \dots, a_n) \in k^n : f_i(a_1, \dots, a_n) = 0, \text{ for all } 1 \le i \le s\} $$

then $V(f_1, \dots, f_s)$ is called an **affine variety** defined by the polynomials $f_i$ (note that in many other text books, *affine varieties* are defined the same way but when $k$ is algebraically closed, e.g. $k = \mathbb{C}$). It is the set of all solutions with coordinates in $k$ of the system of equations $f_1(x_1, \dots, x_n) = \dots = f_s(x_1, \dots, x_n) = 0$.

There is a dual relation between certain varieties and certain ideals (see Cox et al (2007) for the precise conditions on them):

$$ \text{varieties} \leftrightarrow \text{ideals} \\ V = V(I) \leftrightarrow I = I(V), $$

where $V(I)$ is the set of all points in $k^n$ at which all polynomials $f \in \mathbb{I}$ vanish, and $I(V)$ is the set of all polynomials in $k[x_1, \dots, x_n]$ vanishing at all points in $V$.

If $I$ and $J$ are two ideals in $k[x_1, \dots, x_n]$ and $I \subset J$ (that is, $J$ contains more polynomials than $I$), then $V(J) \subset V(I)$. This is because the points where all polynomials in $J$ vanish make also all the polynomials in $I$ vanish, but the additional polynomials in $J$ impose additional conditions, so not all points in $V(I)$ are necessarily in $V(J)$ (more equations mean less solutions).

### Difference of varieties duality

*Difference of varieties* $\xleftarrow{\text{duality}}$ *quotient of ideals*

The **set difference** of two affine varieties is generally not an affine variety. For instance in $\mathbb{C}^2$, the difference $V(\langle x \rangle) \setminus V(\langle y \rangle)$ is the set resulting from removing the $x$-axis from the $y$-axis. This leaves the $y$-axis minus the point $\{0, 0\}$, which is not a variety but an open subset of a variety. It cannot be written as the set of solutions of a system of polynomial equations (it is not an affine variety). The smallest affine variety which contains it, is called the **Zariski closure** of the difference, denoted with an overline. Loosely speaking, taking the Zariski closure amounts to patching up the holes in the open set. Therefore, in this example, the Zariski closure of the difference is the $y$-axis itself:

$$ \overline{V(\langle x \rangle) \setminus V(\langle y \rangle)} = V(\langle x \rangle). $$

Using the correspondence between ideals and varieties, an ideal defining $\overline{V} \setminus \overline{W}$ where $V$ and $W$ are affine varieties is obtained as the saturation of an ideal $I$ defining $V$ with an ideal $J$ defining $W$. It is

denoted by $I: J^\infty$. The colon symbolizes the **quotient** of one ideal w.r.t. the other, and it removes factors from the polynomials in the first ideal which appear as polynomials themselves in the second ideal. The infinity symbol corresponds to considering all products of arbitrarily many polynomials in the second ideal. A couple of examples:

$$ \langle x \rangle : \langle y \rangle^\infty = \langle x \rangle, \\ \langle xy^2 \rangle : \langle y \rangle = \langle xy \rangle \subset \langle xy^2 \rangle : \langle y^2 \rangle = \langle xy^2 \rangle : \langle y \rangle^\infty = \langle x \rangle. $$

Then,

$$ \langle x \rangle : \langle y \rangle^\infty = \langle x \rangle. \leftrightarrow \overline{V(\langle x \rangle) \setminus V(\langle y \rangle)} = V(\langle x \rangle), \\ \langle xy^2 \rangle : \langle y \rangle^\infty = \langle x \rangle. \leftrightarrow \overline{V(\langle xy \rangle) \setminus V(\langle y \rangle)} = V(\langle x \rangle). $$

In the context of this paper, to find out whether actually $V(I_4) = V(I_{28})$ or the first one strictly contains the second one, we are interested in computing the Zariski closure of the difference $\overline{V(I_4)} \setminus \overline{V(I_{28})}$. Dually, determining this set amounts to computing the saturation of the ideal $I_4$ w.r.t. the ideal $I_{28}$, denoted by $I_4: I_{28}^\infty$. This is equivalent to removing from the polynomials in $I_4$ all the factors which are themselves polynomials in $I_{28}$. It can be done in a computational environment such as Maple or Singular and it results in (25), where

$$ I_d = I_4 : I_{28}^\infty, \text{ and } V(I_d) = \overline{V(I_4)} \setminus \overline{V(I_{28})}. $$

## Gröbner Basis

We can find different sets of polynomial equations having the same solution set. An especially useful one is the **Gröbner basis**. A Gröbner Basis of an ideal $I$ can be defined as a finite set generating the ideal and such that the leading terms of the elements in this set (according to some prescribed ordering of the monomials$^6$) generate the leading terms of all elements in the ideal $I$. One of its many uses is to solve polynomial equations. To do so, it requires an ordering of the terms in the polynomials. The ordering used in this paper is **lexicographic ordering** defined formally as follows:

Let $\alpha = (\alpha_1, \dots, \alpha_n)$ and $\beta = (\beta_1, \dots, \beta_n) \in \mathbb{Z}_{\ge 0}^n$. It is said that $\alpha >_{lex} \beta$ if, in the vector difference $\alpha - \beta \in \mathbb{N}$, the leftmost nonzero entry is positive. The monomials $x^\alpha >_{lex} x^\beta$ if $\alpha >_{lex} \beta$.

## C. Validity of Gröbner basis of the cylinder polynomials

The fact that the set of polynomials $\{F_1, F_2, F_3\}$ in (29-31) is a Gröbner basis of the ideal $\langle C_{123}, C_{124}, C_{134}, C_{234} \rangle$ implies that the polynomial system $\{F_1 = F_2 = F_3 = 0\}$ has the same solutions as the polynomial system $\{C_{123} = C_{124} = C_{134} = C_{234} = 0\}$.

When one has, as in this case, two sets of generators for an ideal, then one can write the polynomials in the first set as a linear combination of the polynomials in the second set with coefficients in the real numbers and the parameters, and vice versa. On one side, the fact that the set $\{F_1, F_2, F_3\}$ can be expressed as a combination of the set $\{C_{123}, C_{124}, C_{134}, C_{234}\}$ means that for $i = 1, 2, 3$:

$$ F_i = \frac{\alpha_i}{\alpha'_i} C_{123} + \frac{\beta_i}{\beta'_i} C_{124} + \frac{\gamma_i}{\gamma'_i} C_{134} + \frac{\delta_i}{\delta'_i} C_{234}, \qquad (52) $$

where $\alpha_i, \beta_i, \gamma_i, \delta_i, \alpha'_i, \beta'_i, \gamma'_i, \delta'_i$ are some polynomials in the variables $d_1, \dots, d_5$ with real coefficients. However, we have computed a Gröbner basis with parameters in the coefficients, and this means that, for some choices of the parameter values, the denominators $\alpha'_i, \beta'_i, \gamma'_i, \delta'_i$ could

<footnote>$^6$ For a polynomial of degree $m$ in one variable, $f = a_0x^m + a_1x^{m-1} + \dots + a_m$, $a_0x^m$ is the leading term. The reader is directed to Cox et al (2007) for further details on Gröbner basis and definition of the leading term in case of multi-variate polynomials.</footnote>
be zero. For such values of the parameters, the Gröbner basis we have
obtained will not be valid. Similarly, for $i, j, k \in \{1, 2, 3, 4\}$,

$$
C_{ijk} = \frac{A_{ijk}}{A'_{ijk}} F_1 + \frac{B_{ijk}}{B'_{ijk}} F_2 + \frac{C_{ijk}}{C'_{ijk}} F_3 \quad (53)
$$

with $A_{ijk}$, $B_{ijk}$, $C_{ijk}$, $A'_{ijk}$, $B'_{ijk}$, $C'_{ijk} \in \mathbb{R}[d_1, \dots, d_5]$. For the same reason as before, whenever the choice of the parameter values is such that $A'_{ijk}$, $B'_{ijk}$ or $C'_{ijk}$ vanish for some $i, j, k$, then $C_{ijk}$ cannot be written as a linear combination of $F_1, F_2, F_3$. Computing the coefficients of these combinations, we come to the conclusion that the Gröbner basis is valid as long as the following polynomial given by the least common multiple of the denominators in (52) and (53) does not vanish:

$$
q = \operatorname{lcm}(\alpha'_1, \beta'_i, \gamma'_i, \delta'_i, A'_{ijk}, B'_{ijk}, C'_{ijk}), \quad i, j, k \in \{1, 2, 3, 4\}. \quad (54)
$$

In fact, the above calculations can be avoided due to the fact that the product $c_6 c_y c_z$ (see 29-31) contains all the factors of $q$ (except the trivial factors corresponding to the coincidence of some $M_i$ and $M_j$, $i \neq j = 1, 2, 3, 4$) and hence for any parameter choice for which $c_6 c_y c_z = 0$, the calculated Gröbner basis is not valid (see the supplementary Maple file). We refer to Kalkbrener (1997) for more involved results on the stability of Gröbner bases under specialization.

References

Ben-Horin P, Shoham P (2006) Singularity analysis of a class of par-
allel robots based on Grassmann-Cayley algebra. Mechanism and
Machine Theory 41(8):958-970

Bézout E (1779) Théorie générale des équations algébriques. (Paris)

Briot S, Martinet P (2013) Minimal representation for the control of Gough-Stewart platforms via leg observation considering a hidden robot model. In: Proceedings of the 2013 IEEE International Conference on Robotics and Automation (ICRA 2013), Karlsruhe, Germany

Briot S, Martinet P, Rosenzveig V (2015) The hidden robot: an efficient concept contributing to the analysis of the controllability of parallel robots in advanced visual servoing techniques. IEEE Transactions on Robotics 31(6):1337–1352

Briot S, Rosenzveig V, Martinet P, Özgür E, Bouton N (2016) Minimal representation for the control of parallel robots via leg observation considering a hidden robot model. Mechanism and Machine Theory 106:115–147

Briot S, Chaumette F, Martinet P (2017a) Revisiting the determination of the singularity cases in the visual servoing of image points through the concept of “hidden robot”. IEEE Transactions on Robotics 33(2):536–546

Briot S, Martinet P, Chaumette F (2017b) Singularity cases in the visual servoing of three image lines. IEEE Robotics and Automation Letters 2(2):412–419

Burschka D, Geiman J, Hager G (2003) Optimal landmark configura-
tion for vision-based control of mobile robots. Proceedings - IEEE
International Conference on Robotics and Automation 3:3917–3922

Chaumette F, Hutchinson S (2008) Visual Servoing and Visual Tracking, chapter 24 of Handbook of Robotics. Springer

Corke P (2010) Spherical image-based visual servo and structure estima-
tion. In: Proceedings of the 2010 IEEE International Conference on
Robotics and Automation (ICRA 2010), Anchorage, pp 5550–5555

Cox DA, Little J, O'Shea D (2007) Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, 3/e (Undergraduate Texts in Mathematics). Springer-Verlag, Berlin, Heidelberg

DeMenthon D, Davis L (1995) Model-based object pose in 25 lines of codes. International Journal of Computer Vision 15:123–141

Gao XS, Hou XR, Tang J, Cheng HF (2003) Complete solution classification for the perspective-three-point problem. IEEE Transactions on Pattern Analysis and Machine Intelligence 25(8):930–943

Hamel T, Mahony R (2002) Visual servoing of an under-actuated dy-
namic rigid-body system: An image based approach. IEEE Transac-
tions on Robotics and Automation 18:187–198

Hamel T, Samson C (2017) Riccati observers for the non-stationary PnP problem. IEEE Transactions on Automatic Control 63(3):726–741

Horaud R, Conio B, Leboulleux O, Lacolle B (1989) An analytic solution for the perspective 4-point problem. Computer Vision, Graphics, and Image Processing 47(1):33–44

Hutchinson S, Hager G, Corke P (1996) A tutorial on visual servo control. IEEE Transactions on Robotics and Automation 12:651–670

Iwatsuki M, Okiyama N (2005) A new formulation for visual servo-
ing based on cylindrical coordinate system. IEEE Transactions on
Robotics 21(2):266–273

Kalkbrener M (1997) On the stability of gröbner bases under specializa-
tions. Journal of Symbolic Computation 24(1):51–58

Kneip L, Furgale P (2014) Opengv: A unified and generalized approach to real-time calibrated geometric vision. In: 2014 IEEE International Conference on Robotics and Automation (ICRA), IEEE, pp 1–8

Kneip L, Li H, Seo Y (2014) Upnp: An optimal o (n) solution to the absolute pose problem with universal applicability. In: European Conference on Computer Vision, Springer, pp 127–142

Lepetit V, Moreno-Noguer F, Fua P (2009) Epnp: An accurate o(n) so-
lution to the pnp problem. International Journal of Computer Vision
81(2):155

Liu M, Pradalier C, Siegwart R (2010) A bearing-only 2D/3D-homing method under a visual servoing framework. In: Proceedings of the 2010 IEEE International Conference on Robotics and Automation (ICRA 2010), Anchorage, pp 4062–4067

Longuet-Higgins HC, Prazdny K (1980) The interpretation of a moving retinal image. Proceedings of the Royal Society of London, Section B 208:385–397

Lowe DG (1991) Fitting parameterized three-dimensional models to im-
ages. IEEE Transactions on Pattern Analysis & Machine Intelligence
13(5):441–450

Lu C, Hager G, Mjolness E (2010) Fast and globally convergent pose es-
timation from video images. IEEE Transactions on Pattern Analysis
and Machine Intelligence 22(6):610–622

Marchand E, Chaumette F (2002) Virtual visual servoing: A frame-
work for real-time augmented reality. Computer Graphics Forum
21(3):289–298

Marchand E, Spindler F, Chaumette F (2005) Visp for visual servoing: a generic software platform with a wide class of robot control skills. IEEE Robotics and Automation Magazine 12(4):40–52

Marchand E, Uchiyama H, Spindler F (2016) Pose estimation for aug-
mented reality: a hands-on survey. IEEE Trans on Visualization and
Computer Graphics 22(12):2633–2651

Merlet JP (2006) Parallel Robots, 2nd edn. Springer

Michel H, Rives P (1993) Singularities in the determination of the situation of a robot effector from the perspective view of 3 points. Tech. rep., INRIA

Papanikolopoulos N (1995) Selection of features and evaluation of visual measurements during robotics visual servoing tasks. Journal of Intelligent and Robotics Systems 13(3):279–304

Persson M, Nordberg K (2018) Lambda twist: An accurate fast robust perspective three point (p3p) solver. In: Proceedings of the European Conference on Computer Vision (ECCV), pp 318–332

Rosenzveig V, Briot S, Martinet P, Özgür E, Bouton N (2014) A method for simplifying the analysis of leg-based visual servoing of parallel robots. In: Proc. 2014 IEEE Int. Conf. on Robotics and Automation (ICRA 2014), Hong Kong, China

Tahri O, Chaumette F (2005) Point-based and region-based image mo-
ments for visual servoing of planar objects. IEEE Transactions on
Robotics 21(6):1116–1127

Tatsambon Fomena R, Tahri O, Chaumette F (2011) Distance-based and orientation-based visual servoing from three points. IEEE Transactions on Robotics 27(2):256–267
Thompson E (1966) Space resection: failure cases. The Photogrammet-
ric Record 5(27):201-204