Beatriz Pascual-Escudero
email: beatriz.pascual.escudero@gmail.com

Abhilash Nayak
email: abhilash.nayak@ls2n.fr

Sébastien Briot
email: sebastien.briot@ls2n.fr

Olivier Kermorgant
email: olivier.kermorgant@ls2n.fr

Philippe Martinet
email: philippe.martinet@inria.fr

• Mohab

Din • Franc ¸ois Safey El

Chaumette

M Safey

El Din

F Chaumette
email: francois.chaumette@inria.fr

Complete Singularity Analysis for the Perspective-Four-Point Problem

Keywords: Pose estimation, Visual servoing, Singularity, PnP

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,

Introduction

Pose estimation is a classical problem in computer vision with many applications in augmented reality [START_REF] Marchand | Pose estimation for augmented reality: a hands-on survey[END_REF] and robotics, especially in visual servoing [START_REF] Hutchinson | A tutorial on visual servo control[END_REF]. 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) [START_REF] Horaud | An analytic solution for the perspective 4-point problem[END_REF][START_REF] Gao | Complete solution classification for the perspective-three-point problem[END_REF]. While many efficient methods are available for a long time [START_REF] Lowe | Fitting parameterized three-dimensional models to images[END_REF][START_REF] Dementhon | Model-based object pose in 25 lines of codes[END_REF][START_REF] Lu | Fast and globally convergent pose estimation from video images[END_REF][START_REF] Lepetit | Epnp: An accurate o(n) solution to the pnp problem[END_REF]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 [START_REF] Lowe | Fitting parameterized three-dimensional models to images[END_REF]:

p = arg min p n i=1 (x i -Π c M o (p) X i ) 2
(1)

where p are the six pose parameters representing the change of frame c M o between the object frame and the camera frame, x i are the coordinates of the n image points measured, and X i are the known coordinates of their corresponding 3D-points expressed in the object frame, which are, after the change of frames c 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 ∂x i ∂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 S E 3 , this Jacobian matrix relates the image point velocity to the relative camera-object velocity, leading to the basic motion field equations [START_REF] Longuet-Higgins | The interpretation of a moving retinal image[END_REF]. In that case, the Jacobian matrix is equivalent to the interaction matrix at the basis of any image-based visual servoing scheme [START_REF] Chaumette | Visual Servoing and Visual Tracking, chapter 24 of Handbook of Robotics[END_REF]. Using this form allows simple derivations and avoids any pose representation ambiguity [START_REF] Marchand | Virtual visual servoing: A framework for real-time augmented reality[END_REF].

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 [START_REF] Dementhon | Model-based object pose in 25 lines of codes[END_REF][START_REF] Lepetit | Epnp: An accurate o(n) solution to the pnp problem[END_REF]Kneip et al, 2014), [START_REF] Persson | Lambda twist: An accurate fast robust perspective three point (p3p) solver[END_REF]). 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 [START_REF] Michel | Singularities in the determination of the situation of a robot effector from the perspective view of 3 points[END_REF], 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 [START_REF] Papanikolopoulos | Selection of features and evaluation of visual measurements during robotics visual servoing tasks[END_REF] and more recently in [START_REF] Fomena | Distance-based and orientation-based visual servoing from three points[END_REF].

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 ≥ 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 [START_REF] Iwatsuki | A new formulation for visual servoing based on cylindrical coordinate system[END_REF], the redundant spherical coordinates [START_REF] Hamel | Visual servoing of an under-actuated dynamic rigid-body system: An image based approach[END_REF] C (camera center)

Image plane

A: Cylinder of singularities

M 1 M 2 M 3 m 1 m 2 m 3
Fig. 1: Cylinder of singularities when observing 3 points. [START_REF] Burschka | Optimal landmark configuration for vision-based control of mobile robots[END_REF][START_REF] Corke | Spherical image-based visual servo and structure estimation[END_REF], or other angular representations [START_REF] Liu | A bearing-only 2D/3D-homing method under a visual servoing framework[END_REF], does not modify the problem. It is also the case for non linear combinations of the image points, such as moments [START_REF] Tahri | Point-based and region-based image moments for visual servoing of planar objects[END_REF]. 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 [START_REF] Briot | Minimal representation for the control of Gough-Stewart platforms via leg observation considering a hidden robot model[END_REF], in which it was proven to be efficient to determine:

the singularity cases of vision-based controllers dedicated to parallel robots [START_REF] Rosenzveig | A method for simplifying the analysis of leg-based visual servoing of parallel robots[END_REF][START_REF] Briot | The hidden robot: an efficient concept contributing to the analysis of the controllability of parallel robots in advanced visual servoing techniques[END_REF][START_REF] Briot | Minimal representation for the control of parallel robots via leg observation considering a hidden robot model[END_REF]. 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 [START_REF] Ben-Horin | Singularity analysis of a class of parallel robots based on Grassmann-Cayley algebra[END_REF] and/or the Grassmann geometry [START_REF] Merlet | Parallel Robots[END_REF], 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 [START_REF] Hamel | Riccati observers for the non-stationary PnP problem[END_REF]. 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.

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).

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 [START_REF] Michel | Singularities in the determination of the situation of a robot effector from the perspective view of 3 points[END_REF]. A 3-D point M i with coordinates [X i Y i Z i ] T in the camera frame 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 = X i Z i , y i = Y i Z i (2)
The relation between the image point velocities ẋi and ẏi and the relative camera-object velocity τ c expressed in the camera frame is obtained by differentiation of (2) and is

given by [START_REF] Longuet-Higgins | The interpretation of a moving retinal image[END_REF][START_REF] Chaumette | Visual Servoing and Visual Tracking, chapter 24 of Handbook of Robotics[END_REF] ẋi ẏi = L i τ c

where

L i = -1 Z i 0 x i Z i x i y i -(1 + x 2 i ) y i 0 -1 Z i y i Z i 1 + y 2 i -x i y i -x i (3)
Then, considering the observation of n points (n ≥ 3) M 1 , M 2 , ...., M n , the interaction matrix linking the velocities of the points m i (i = 1, 2, .., n) grouped in the vector ṡ = [ ẋ1 ẏ1 ẋ2 ẏ2 ..., ẋn ẏn ] T to the camera velocity τ c by the relation

ṡ = L (n) τ c (4)
is given by

L (n) =                   L 1 L 2 . . . L n                   =                                     -1 Z 1 0 x 1 Z 1 x 1 y 1 -(1 + x 2 1 ) y 1 0 -1 Z 1 y 1 Z 1 1 + y 2 1 -x 1 y 1 -x 1 -1 Z 2 0 x 2 Z 2 x 2 y 2 -(1 + x 2 2 ) y 2 0 -1 Z 2 y 2 Z 2 1 + y 2 2 -x 2 y 2 -x 2 . . . . . . . . . . . . . . . . . . -1 Z n 0 x n Z n x n y n -(1 + x 2 n ) y n 0 -1 Z n y n Z n 1 + y 2 n -x n y n -x n                                     (5)
We recall that this matrix is equivalent to the Jacobian matrix ∂s ∂p involved in the minimization of (1).

Singularities appear when the matrix L (n) is rank deficient, that is, when rank(L (n) ) < 6. It amounts to the vanishing of all (6 × 6) minors of 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.

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 ≥ 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 L passing through a point M takes the following form:

L : f m(Q) =       f --→ QM × f       (6)
in which f is the direction of the line, and m(Q) = --→ QM × f is the moment of the line with respect to point Q. Now, let us analyze the expression of the interaction matrix L i given in (3). It is obtained by first differentiating the perspective projection equations (2):

ẋi ẏi = ( Ẋi Z i -X i Żi )/Z 2 ( Ẏi Z i -Y i Żi )/Z 2 = 1 Z 2 i Z i 0 -X i 0 Z i -Y i           Ẋi Ẏi Żi           , (7)
and then by using the kinematics equation

          Ẋi Ẏi Żi           = -v C -ω C × ---→ CM i (8)
where v C and ω C are respectively the translational and rotational components of the camera velocity

τ C = [v T C ω T C ]
T . This relation can be rewritten under matrix form:

          Ẋi Ẏi Żi           = -I 3 [ ---→ CM i ] × τ C (9)
where [ ---→ CM i ] × is the antisymmetric matrix related to vector ---→ CM i . By injecting ( 9) in ( 7), we immediately obtain

ẋi ẏi = L i τ C =             f T 1 ---→ CM i × f 1 T f T 2 ---→ CM i × f 2 T             τ C ( 10
)
where

f T 1 = 1 Z 2 i [-Z i 0 X i ] and f T 2 = 1 Z 2 i [0 -Z i Y i ].
The two rows of the interaction matrix 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

---→ CM i = [X i Y i Z i ] is orthogonal to f 1 and f 2 ).
This result was extended in (Briot et al, 2017a) in which we proved that the basis ξ i of the rows of the matrix L i can be written as:

ξ i = ξ i1 ξ i2 (11) in which ξ i j = [f T i j m i j (Q) T ] ( j = 1, 2)
, where f i j is any vector orthogonal to ---→ CM i , and m i j (Q) = ---→ QM i × f i j , Q being any arbitrary point. ξ i j is the Plücker representation of any line L i j of direction f i j contained in the plane P i , which is the plane passing through M i and whose normal is parallel to ---→ CM i (see Fig. 2). A necessary and a sufficient condition for the basis ξ i to be of full rank is that f i1 should not be parallel to f i2 .

The fact that the basis ξ i for the rows of the matrix L i is spanned by a system of two lines allows mentioning the following remarks/results:

-Even if the matrix 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 [START_REF] Merlet | Parallel Robots[END_REF]. However, when providing an analytical expression of the basis ξ 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).

C M 1 M 2 M 3 D 1 D 3 D 2 P 1 P 2 P 3 P 123
Fig. 2: Representation of the planes P 1 , P 2 , P 3 and P 123 .

C M 1 M 2 M 3 D 1 P 123 D 2 D 3 P 1 P 3 P 2 O R
Fig. 3: Condition of singularity for the observation of three points: the four planes planes P 1 , P 2 , P 3 and P 123 intersect in a point O.

-There are infinitely many lines L i j spanning the basis ξ 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 P 1 , P 2 , P 3 and the plane P 123 containing the three points intersect in a point, that we denote by O (see Fig. 3). If the three points are aligned, P 123 is not defined anymore, nonetheless, the planes P 1 , P 2 and P 3 still intersect in a point at infinity. In the general case, it

C' (X, Y, 0) M 3 M 2 M 1 I 1 I 3 I 2 M' 4 1-X, -Y, 0 ( ) d 2 d 1 2 +d 2 2 -d 1 O
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 .

is interesting to note that point O was proven to be on the circle R passing through M 1 , M 2 and M 3 . Moreover, it is contained on the diameter of 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 ξ (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 ξ i for the observation of the point M i (i = 1, 2, 3) the two lines L i j defined by:

f i1 = n 123 × ---→ CM i , f i2 = ---→ CM i × f i1 , (12)
m i1 (O) = ---→ OM i × f i1 , m i2 (O) = ---→ OM i × f i2 , ( 13
)
where n 123 is the normal vector to plane P 123 . By doing this, the line L i1 contains both M i and O, ensuring m i1 (O) = 0. As a result, in singularity, the basis ξ (3) , when computing the lines expressions at point O, takes the following blocktriangular form (after reordering the columns):

ξ (3) ≡                           ξ 11 ξ 21 ξ 31 ξ 12 ξ 22 ξ 32                           =                           f T 11 0 f T 21 0 f T 31 0 f T 12 m T 12 f T 22 m T 22 f 32 m T 32                           (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 ξ (3) of L (3) still takes the form of Eq. ( 14). Moreover, the basis for the rows of the matrix L (4) , defined by:

ξ (4) = ξ T 11 ξ T 21 ξ T 31 ξ T 12 ξ T 22 ξ T 32 ξ T 41 ξ T 42 T (15)
in which the expressions of ξ i1 and ξ i2 are provided at (12) and ( 13), takes the following form:

ξ (4) ≡                                      f T 11 0 f T 21 0 f T 31 0 f T 12 m T 12 f T 22 m T 22 f 32 m T 32 f T 41 m T 41 f T 42 m T 42                                      (16)
It should be mentioned that, when ---→ CM i ⊥ P 123 (i = 1, . . . , 4), then f i1 = 0 and ξ i1 = 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.

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 ξ (4) .

Parametrization of the configuration

Consider the fixed object frame F P : (P, x P , y P , z P ) where the vectors x P , y P and 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 x P axis, M 3 in the x P 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 z P = 0. Consequently, in frame F P , the four observed points have the following coordinates (see Fig. 5):

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 , (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:

d 1 = d 2 = 0 (M 1 ≡ M 3 ), d 1 -1 = d 2 = 0 (M 2 ≡ M 3 ), d 3 = d 4 = d 5 = 0 (M 1 ≡ M 4 ), d 3 -1 = d 4 = d 5 = 0 (M 2 ≡ M 4 ), d 3 -d 1 = d 4 -d 2 = d 5 = 0 (M 3 ≡ M 4 ).
We will consider separately the configuration where the four points are coplanar (in Section 3.2.2), that is, when

either d 5 = 0 or d 2 = 0, ( 18
)
from the general configuration (in Section 3.2.1).

Another special case arises when ---→ CM i ⊥ P 123 for some i = 1, . . . , 4.

(

) 19
Note that, from the definition of F P , the normal of plane 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, ξ i1 = 0, hence we should make a different choice of the line direction f i1 in order for ξ (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 ---→ CM 1 is perpendicular to the plane 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 ---→
CM 1 is not perpendicular to the new plane 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 ---→ CM i , i = 1 . . . 4 is perpendicular to P 123 , it will not be perpendicular to the three remaining planes P 124 , P 134 or P 234 . Therefore, permuting the points such that the plane P 123 is replaced with one of those planes makes ---→ CM i no longer perpendicular to the latter. This fact can be exploited so that the matrix ξ (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

F P by C = [X Y Z] T (20)
We will determine the coordinates (X, Y, Z) of the camera center in 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 camera frame F C where the perspective equations (2) apply, the following change of frame might be necessary: where c R P is the rotation matrix between frames F P and F C and M i are the coordinates of point M i in F P given in ( 17).

          X i Y i Z i           = c R P M i -           X Y Z           (21) M 1 (0, 0, 0) D 4 D 3 D 2 D 1 x P z P M 2 (1, 0, 0) M 3 (d 1 , d 2 , 0) M 4 (d 3 , d 4 , d 5 ) C (X, Y, Z)
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:

X = Y = Z = 0 (C ≡ M 1 ), X -1 = Y = Z = 0 (C ≡ M 2 ), X -d 1 = Y -d 2 = Z = 0 (C ≡ M 3 ), X -d 3 = Y -d 4 = Z -d 5 = 0 (C ≡ M 4 ).
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 calculations). In such a case, as already mentioned, there exists a point O where the planes P 1 , P 2 , P 3 and 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 orthogonal projection C of the camera center C over the plane P 123 . Knowing the coordinates of points M 1 , M 2 and M 3 in the frame F P , the coordinates of point O can be calculated in the same frame. They are given by (see Fig. 4)

O =       1 -X, d 2 1 + d 2 2 -d 1 d 2 -Y, 0       T . (22)
In the next section, we determine the singularities of the matrix ξ (4) .

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 ξ (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.

General case (non coplanar points)

From ( 18), we consider here that

d 2 = 0 and d 5 = 0. ( 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 , . . . , d 5 . We aim at giving equations for the singular positions [X Y Z] of the camera depending on the parameters d 1 , . . . , d 5 .

The most direct way of obtaining these equations is by computing the minors of the matrix ξ (4) : whenever the 28 minors of size (6×6) of the matrix ξ (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)) prescribes1 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 ξ (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, ( 24
)
d 2 1 -2d 1 X + d 2 2 -2d 2 Y + X 2 + Y 2 , d 2 3 -2d 3 X + d 2 4 -2d 4 Y + X 2 + Y 2 .
In fact, all of these factors can be rewritten as (X - For each choice (d 1 , d 2 , d 3 , d 4 , d 5 ) ∈ R5 of the design parameters, the values of camera coordinates (X, Y, Z) ∈ R3 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 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 R 3 to all polynomials in the ideal I 4 .

X M i ) 2 + (Y -Y M i ) 2 ,
In general, when we write V(I 4 ) and V(I 28 ) we will mean varieties in R 3 as just explained. However, for certain tasks we will consider all d 1 , ..., 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 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 ⊂ I 28 , and it turns out to be computationally easier to handle. However, the variety 2 V(I 4 ) might contain more points than V(I 28 ) ( V(I 28 ) ⊆ V(I 4 ) ⊂ R 8 , see Appendix B). The following ideal I d computed in Maple (see Appendix B and the supplementary Maple file) 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 = d 5 , (25)
-d 2 1 d 4 -d 2 2 d 4 + d 2 d 2 3 + d 2 d 2 4 + d 1 d 4 -d 2 d 3 , (26) d 2 1 Y + d 2 2 Y -d 2 X 2 -d 2 Y 2 -d 1 Y + d 2 X, (27)
d 2 3 Y + d 2 4 Y -d 4 X 2 -d 4 Y 2 -d 3 Y + d 4 X . ( 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 = 0 as assumed in the beginning of this section,

V(I d ) = ∅, implying V(I 4 ) = V(I 28 ).
Consequently, the real solutions of the minors' system g 1 = 0, . . . , 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 variety3 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 [START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, 3/e[END_REF]). The computer algebra system Maple is used to do so (see the supplementary Maple file). In particular, by choosing the lexicographical ordering4 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 , . . . , 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 , (29)
F 2 (X, Y) = c y Y + b 2 (X), (30)
F 3 (X, Z) = c z Z + b 3 (X). ( 31
)
where c z , c y , c j , j = 1, ..., 6 are polynomial expressions in the parameters d 1 , . . . , d 5 with real coefficients, while b 2 , b 3 are polynomials in X, whose coefficients are themselves polynomial expressions in d 1 , . . . , d 5 with real coefficients. Let us note that c 6 = 16d 4 2 d 4 5 = 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 = 0 (remember c 6 = 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} ⊂ 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 i jk 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 = 0, the singular positions for the camera will be those points (X S , Y S , Z S ) ∈ R 3 such that

F 1 (X S ) = F 2 (X S , Y S ) = F 3 (X S , Z S ) = 0. ( 32
)
Indeed, given a fixed set of four points (or equivalently, a fixed set of values for the parameters d 1 , . . . , 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 , . . . , 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 , . . . , d 5 such that d 2 = 0, d 5 = 0, there are at most 6 singular positions of the camera when observing four non coplanar points.

One could wonder whether there is a choice of the parameter 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 supplementary Maple file, some of the coefficients c j , j = 1, ..., 6 are dependent, which allows eliminating some of them as follows:

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 . ( 33
)
Since c 6 = 16d 4 2 d 4 5 = 0 (according to ( 23)), we can normalize F 1 , obtaining

F 1 = X 6 -3 X 5 + (c 2 + c 1 + 3 c 6 ) c 6 X 4 - (2c 1 + 2c 2 + c 6 ) c 6 X 3 + c 2 c 6 X 2 + c 1 c 6 X + c 0 c 6 . ( 34
)
Consider now the variable change (known as the linear Tschi-

rnhausen transformation) X = t - c 5 6c 6 = t + 1 2
. This change on 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:

F 1 (X) = G 1 (t) = t 6 + C 4 t 4 + C 2 t 2 + C 0 . ( 35
)
The latter is a cubic equation in t 2 = T :

G1 (T ) = T 3 + C 4 T 2 + C 2 T + C 0 . ( 36
)
It has to be noted that Z < lex Y < lex X is the only variable order for which the term c 5 6c 6 = 1 2 in Tschirnhausen transformation is a constant. This helps us to geometrically interpret the solutions to F 1 = 0. Let T 1 be a real positive root of G1 (T ) = 0. Then, t = ± √ T 1 are roots of G 1 (t) and hence two roots of the polyno-

mial F 1 (X) are X 1 = √ T 1 + 1 2 and X 2 = - √ T 1 + 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 , Y 2 = 2Y * -Y 1 , Z 2 = 2Z * -Z 1 with X * = 1 2 , Y * = d 2 1 + d 2 2 -d 1 2d 2 , Z * = - d 4 (d 2 1 + d 2 2 ) -d 2 (d 2 3 + d 2 4 + d 2 5 ) -d 1 d 4 + d 2 d 3 2d 2 d 5 ,
where (X * , Y * , Z * ) is the center of the circum-sphere of the tetrahedron M 1 M 2 M 3 M 4 . Thus, it can be deduced that every 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

6 6 6 2 4 -2 -2 -2 -2 -2 -2 -2 -2 -4 -4 -4 -4 -4 -4 -4 x x x y -6 -6 -8 -4 -4 -4 -6 -6 2 2 2 2 2 -2 z M 3 (-2, -3, 0) M 1 M 2 C 1 (0.65, 0.05, -7.47) C' 1 M 4 (2, -1, -7)
C 2 (0.35, -5.05, 0.75)

C' 2 O 1 C' 1 (0.65, 0.05) M 3 M 2 M 1 C' 4 (1.34, -0.09) C' 2 O 1 (0.35, -5.05) C' 4 O 2 (-0.34, -4.91) C' 6 O 5 (-1.62, -3.92) C' 5 (2.62, -1.08) X Y P 123 P 123
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 P 123 (bottom).

+ can be found on the sphere containing this point and with center (X * , Y * , Z * ). Moreover, the projection of these solutions into the plane 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 configuration 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 P 123 as C 1 and C 2 ≡ O 1 , O 1 being the antipodal point of C 1 . In this particular case, there exist six singular positions for the camera. Their orthogonal projection onto the plane P 123 is also shown on this figure.

The only case in which F 1 will have real roots is when G1 (T ) has at least one positive real root. It turns out that for any values of d 1 , . . . , d 5 , G1 (0) = -4p 2 ≤ 0, where p is a common factor of polynomials c y and c z that appears in q (see the supplementary Maple file). So G1 (0) < 0 when q = 0. Also, G1 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 G1 , and therefore at least two real roots for F 1 .

Meanwhile, if p = 0 and hence G1 (0) = 0, the cubic in (36) can be factorized as G1 (T ) = T (T 2 + C 4 T + 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 = 1 2 . Geometrically, the line connecting the projections of these solutions in 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 = 1 2 in the first cylinder equation ( 47), which results in the quadratic equation

Y 2 d 2 -(d 2 1 + d 2 2 -d 1 )Y - 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 [START_REF] Chaumette | Algebraic geometry tools Most of the following definitions and examples are taken from Cox et al[END_REF].

As a result, G1 (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 G1 (T ) is given by

∆ = -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 + ... 1600 terms ... + d 2 6 d 3 6 + d 2 6 d 3 4 d 5 2 ). ( 37
)
The last factor of ∆ is difficult to interpret geometrically, while the rest of the algebraic factors also appear in the product c 6 c y c z (see ( 29)-( 31) and Appendix C). As already mentioned, even when q = c y c 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.

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 ξ (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 i j is defined in ( 13). This can be overcome by simply re-defining the Plücker representation ξ i j 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 ξ i j , i = 1, 2, 3, 4, j = 1, 2 such that f i j is chosen according to (12) while the moment m i j of the line is chosen as follows:

m i1 (M 1 ) = ----→ M 1 M i × f i1 , m i2 (M 1 ) = ----→ M 1 M i × f i2 . ( 38
)
Following the same procedure as in Section 3.2, vanishing of all (6 × 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 ], ( 39
)
h a 2 = [Z, X 2 d 4 + Y 2 d 4 -Yd 3 2 -Yd 4 2 -Xd 4 + Yd 3 , X 2 d 2 + Y 2 d 2 -Yd 1 2 -Yd 2 2 -Xd 2 + Yd 1 , -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 ].
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 ], ( 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 a 4 vanish when the camera lies on a line defined by

X = d 3 d 4 (d 4 -Y) + d 5 (d 5 -Z) = 0. (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 f i1 = 0 inducing ξ i1 = 0, in which case the form of ξ (4) given by ( 16) cannot be used anymore. Yet again, a different basis ξ 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 , ---→ CM i is along z P -axis in F P and the direction vectors of the lines ξ i j can be simply chosen as parallel to x P and y P axes as follows:

f i1 = [1, 0, 0] T , f i2 = ---→ CM i × f i1 = [0, 1, 0] T (42)
The moments m i j of the corresponding lines can be calculated using (38) to avoid missing the case d 2 = 0. Then, the minors of the modified basis ξ (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 ⊥ P 123 : {X = Y = 0} (43) h b 1 = [d 5 , -d 4 (d 1 -1) + d 2 (d 3 -1)], CM 2 ⊥ P 123 : {X -1 = Y = 0} h b 2 = [d 5 , d 1 d 4 -d 2 d 3 ], CM 3 ⊥ P 123 : {X -d 1 = Y -d 2 = 0} h b 3 = [d 5 , d 4 ], CM 4 ⊥ P 123 : {X -d 3 = Y -d 4 = 0} h b 4 = [d 5 , d 2 ].
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 -----→ M 2 M 3 and -----→ M 2 M 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 a 4 in (40) with d 2 = 0. Thus, generalizing the result obtained from h a 4 in (40), when the camera is vertically aligned with one of the points while the other three are collinear, the singularities are located

Cases Singularity configurations

Four points non-coplanar (d 2 = 0 and d 5 = 0).

At least 2 and at most 6 real singularities All four points are coplanar, but not concyclic, and no three of them collinear.

No singularities.

All four points are concyclic.

C concyclic with the points. Only three points M i , M j , M k are collinear but not M l , (i, j, k, l all different ∈ 1, . . . , 4) C on the line through M l orthogonal to the plane containing the four points.

All four points are collinear.

Always singular.

Table 1: Different cases of singularities when four points are observed.

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.

Homogenization of the cylinder equations

The cylinder equations ( 47)-( 50) were derived by considering the point M 2 = [1 0 0] T as shown in ( 17), implying that the vector -----→ M 1 M 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 0 0 0] T , [D 1 D 2 0] T , [D 3 D 4 D 5 ] T and [X d Y d Z 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 2 D 2 + X d D 2 D 0 -Y d 2 D 2 + D 1 2 + D 2 2 -D 0 D 1 Y d = 0.
By substituting D i = d i D 0 , i = 1, . . . , 5, X d = XD 0 and Y d = Y D 0 in the above equation, leads to:

D 3 0 -X 2 d 2 + Xd 2 -Y 2 d 2 + d 1 2 + d 2 2 -d 1 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 = 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.

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. 0.9 0.3 0.3 0.9 0.4 (0.94, -0.237, 2.599) 4 0.6 0.4 0.9 0.9 0.3 (0, 0, 2.7) 5 0.5 0.5 0.9 0.9 0.3 (0, 0, 2.7) 6 0.9 0.5 0.3 0.9 0.3 (0.96, -0.092, 2.734) 7 0.5 0.4 0.9 0.9 0.3 (-0.012, -0.120, 3.015) 8 0.8 0.3 0.3 0.9 0.3 (0.988, -0.275, 3.177) 9 0.5 0.3 0.9 0.9 0.3 (-0.061, -0.339, 3.447) 10 0.8 0.3 0.3 0.8 0.3 (0.977, -0.306, 2.648) 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.

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 m × 1 m × 1m; the position of the camera should be at a distance betweeen 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 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:

∀i ∈ [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 o ), and Z i is the depth of the point expressed in the camera frame thanks to (21).

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:

5 Interested reader can perform all simulations using the code available on https://github.com/oKermorgant/singularity 4points the non-planar DeMenthon's method [START_REF] Dementhon | Model-based object pose in 25 lines of codes[END_REF]; the Efficient PnP [START_REF] Lepetit | Epnp: An accurate o(n) solution to the pnp problem[END_REF]; the Universal PnP (Kneip et al, 2014); and a P4P approach based on lambda-twist P3P [START_REF] Persson | Lambda twist: An accurate fast robust perspective three point (p3p) solver[END_REF]. 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 [START_REF] Marchand | Virtual visual servoing: A framework for real-time augmented reality[END_REF]. This approach performs first-order local optimization and relies directly on the interaction matrix; -Levenberg-Marquardt (LM) method [START_REF] Lowe | Fitting parameterized three-dimensional models to images[END_REF] 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 [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF], in which DeMenthon, VVS and Levenberg-Marquardt methods are implemented. Lambda-twist (P4P) implementation relies on the open-source code provided by the authors [START_REF] Persson | Lambda twist: An accurate fast robust perspective three point (p3p) solver[END_REF]. Finally, EPnP and UPnP implementations are from OpenGV [START_REF] Kneip | Opengv: A unified and generalized approach to real-time calibrated geometric vision[END_REF].

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 σ = 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 θ e = arccos tr( c R P c R P T ) -1)/2 where c R P and c R P 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: . 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. 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.

-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 secondorder (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.

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:

           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 (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 δ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 neighboorhood 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 firstorder 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.

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 s(t) and their desired value s * . The most basic controller tries to ensure a decoupled exponential decrease of the error ss * :

ṡ = -λ(s -s * ) ( 45
)
where λ is a positive scalar. The resulting control law has the following form [START_REF] Chaumette | Visual Servoing and Visual Tracking, chapter 24 of Handbook of Robotics[END_REF]:

τ c = -λ L s + (s -s * ) (46)
where τ c is the camera velocity sent to the robot controller, L s an estimation of the interaction matrix of s, and L s + its Moore-Penrose pseudo-inverse. In our case, the feature vector s is composed of the Cartesian coordinates of four image points, so L s = L (4) given in (5). The depths Z i of the points appear in 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.

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), 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 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.

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.

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.

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 ss * 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 ss * is in the null space of 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.

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 ), tion 4.2.2 is on the good side of the global minimum, except for a petal of the rose.

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.

Ideal

A subset I of the ring of polynomials in the variables x 1 , . . . , x n with coefficients in a field k, k[x 1 , ..., x n ] is an ideal if it satisfies:

(i) 0 ∈ I (ii) If f, g ∈ I then f + g ∈ I (iii) If f ∈ I and h ∈ k[x 1 , ..., x n ], then h f ∈ I. (51) Let f 1 , ..., f s be polynomials in k[x 1 , ..., x n ], then f 1 , ..., f s = s i=1 h i f i : h 1 , ..., h s ∈ k[x 1 , ..., x n ]
is an ideal generated by f 1 , ..., f s . This set is interesting because, for any field K containing k, all polynomials f ∈ f 1 , ..., f s give polynomial equations f = 0 that hold at all points of K n at which f 1 = 0, . . . , f s = 0 hold simultaneously, plus it has a structure that allows many useful algebraic operations.

Variety

If f 1 , ..., f s are polynomials in the ring k[x 1 , ..., x n ] and if V( f 1 , ..., f s ) = {(a 1 , ..., a n ) ∈ k n : f i (a 1 , ..., a n ) = 0, for all 1 ≤ i ≤ s} then V( f 1 , ..., 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 = C). It is the set of all solutions with coordinates in k of the system of equations f 1 (x 1 , ..., x n ) = ... = f s (x 1 , ..., x n ) = 0.

There is a dual relation between certain varieties and certain ideals (see [START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, 3/e[END_REF] for the precise conditions on them): varieties ←→ ideals V = V(I) ←→ I = I(V), where V(I) is the set of all points in k n at which all polynomials f ∈ I vanish, and I(V) is the set of all polynomials in k[x 1 , . . . , x n ] vanishing at all points in V.

If I and J are two ideals in k[x 1 , . . . , x n ] and I ⊂ J (that is, J contains more polynomials than I), then V(J) ⊂ 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 ← ----→ quotient of ideals

The set difference of two affine varieties is generally not an affine variety. For instance in C 2 , the difference V( x ) \ V( y ) 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:

V( x ) \ V( y ) = V( x ).
Using the correspondence between ideals and varieties, an ideal defining V \ 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 ∞ . 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: 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 V(I 4 ) \ 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 . 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

x : y ∞ = x ,
I d = I 4 : I ∞
28 , and V(I d ) = V(I 4 ) \ 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 monomials6 ) 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 α = (α 1 , ..., α n ) and β = (β 1 , ..., β n ) ∈ Z n ≥0 . It is said that α > lex β if, in the vector difference αβ ∈ n, the leftmost nonzero entry is positive. The monomials x α > lex x β if α > lex β.

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 C 123 , C 124 , C 134 , C 234 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 = α i α i C 123 + β i β i C 124 + γ i γ i C 134 + δ i δ i C 234 , (52)
where α i , β i , γ i , δ i , α i , β i , γ i , δ i are some polynomials in the variables d 1 , . . . , 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 α i , β i , γ i , δ i could be zero. For such values of the parameters, the Gröbner basis we have obtained will not be valid. Similarly, for i, j, k ∈ {1, 2, 3, 4},

C i jk = A i jk A i jk F 1 + B i jk B i jk F 2 + C i jk C i jk F 3 (53)
with A i jk , B i jk , C i jk , A i jk , B i jk , C i jk ∈ R[d 1 , . . . , d 5 ]. For the same reason as before, whenever the choice of the parameter values is such that A i jk , B i jk or C i jk vanish for some i, j, k, then C i jk 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 = lcm α 1 , β i , γ i , δ i , A i jk , B i jk , C i jk , i, j, k ∈ {1, 2, 3, 4}.

(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 = 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 [START_REF] Kalkbrener | On the stability of gröbner bases under specializations[END_REF] for more involved results on the stability of Gröbner bases under specialization.

Fig. 5 :

 5 Fig. 5: Observation of four points. The coordinates of the points in F P are also marked.

Fig. 7 :

 7 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 σ = 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.

  Fig. 8: (a) 3D configuration for Ex. 1 in Table2. 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.

Fig. 9 :

 9 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).

Fig. 10 :

 10 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.

Fig. 11 :Fig. 12 :Fig. 13 :

 111213 Fig. 11: Visual servoing simulations from 6 different starting poses (colored cameras). Top view is centered on the singularity point from Ex. 1 in Table2. 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.

Fig. 14 :

 14 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.

  xy 2 : y = xy ⊂ xy 2 : y 2 = xy 2 : y ∞ = x . Then, x : y ∞ = x . ↔ V( x ) \ V( y ) = V( x ), xy 2 : y ∞ = x . ↔ V( xy ) \ V( y ) = V( x ).

  where X M i and Y M i are the X and Y coordinates of the point M i , i = 1, . . . , 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 , . . . , 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 i jk is the cylinder corresponding to points M i , M j and M k .The polynomials g 1 , . . . , g 28 generate an ideal (see Appendix B and[START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, 3/e[END_REF] for further reading), denoted by I 28 = g 1 , . . . , g 28 . It is essentially the set of all (linear and nonlinear) combinations of the constituent polynomials using coefficients from the set 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 , . . . , d 5 . Similarly, let I 4 = C 123 , C 124 , C 134 , C 234 be the ideal generated by the polynomials corresponding to the four simplified cylinder equations.

  Ex. d 1 d 2 d 3 d 4 d 5 Singular point (X s , Y s , Z s )

	1	0.6 0.4 0.3 0.9 0.3 (1.009, -0.082, 2.504)
	2	0.7 0.3 0.3 0.8 0.3 (1.019, -0.341, 2.546)
	3

Table 3 :

 3 Desired and starting positions for the visual servoing simulations. Positions are relative to the singularity position. All units are in meters.

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 = 0.

We are considering varieties in R[d 1 , d 2 , d

, d

, d

, X, Y, Z] in this discussion.

Through the rest of this section we are only regarding X, Y, Z as variables, so we are considering varieties in R 3 .

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).

For a polynomial of degree m in one variable, f = a 0 x m + a 1 x m-1 + . . . + a m , a 0 x m is the leading term. The reader is directed to[START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra, 3/e[END_REF] for further details on Gröbner basis and definition of the leading term in case of multi-variate polynomials.

Appendices

A. Cylinders Equations

Equations of the four cylinders C i jk containing points M i , M j , M k :

Thompson E (1966) Space resection: failure cases. The Photogrammetric Record 5(27):201-204