Riccardo Spica
email: riccardo.spica@irisa.fr

Paolo Robuffo Giordano

Franc ¸ois Chaumette
email: francois.chaumette@irisa.fr

F Chaumette

Active Structure from Motion: Application to Point, Sphere and Cylinder

In this paper, we illustrate the application of a nonlinear active structure estimation from motion (SfM) strategy to three problems, namely, 3D structure estimation for (i) a point, (ii) a sphere and, (iii) a cylinder. In all three cases, an appropriate parametrization reduces the problem to the estimation of a single quantity. Knowledge of this estimated quantity and of the available measurements allows for then retrieving the full 3D structure of the observed objects. Furthermore, in the point feature case, two different parametrizations based on either a planar or a spherical projection model are critically compared: indeed, the two models yield, somehow unexpectedly, to different convergence properties for the SfM estimation task. The reported simulative and experimental results fully support the theoretical analysis and clearly show the benefits of the proposed active estimation strategy, which is in particular able to impose a desired transient response to the estimation error equivalent to that of a reference linear second-order system with assigned poles.

I. INTRODUCTION

T HE problem of Structure from Motion (SfM), i.e., how to recover the missing structure of the observed scene from images taken by a moving camera, is a very classical and well-studied topic in computer and robot vision. One solution is to rely on prior knowledge of the scene as, e.g., known size of a tracked object. Alternatively, one can exploit the possibility of observing the same scene from different points of view, and fuse together the acquired images with the known camera displacement among them. When processing consecutive images over time, a possibility is to treat SfM as a recursive/filtering task: images and camera motion can be elaborated online for obtaining an incremental estimation of the scene structure. Other approaches (e.g., bundle adjustment) rely, instead, on global optimization methods meant to solve SfM problems by processing altogether information acquired over an extended time period. A recent discussion about the pros/cons of both approaches in the context of Visual SLAM can be found in [START_REF] Strasdat | Visual SLAM: Why filter?[END_REF].

Within the first class of (recursive) methods, a vast literature exists for addressing SfM: for instance, as a non-exhaustive list, Extended Kalman Filter (EKF)-based solutions have been proposed in [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF]- [START_REF] Civera | 1-Point RANSAC for EKF Filtering: Application to Real-Time Structure from Motion and Visual Odometry[END_REF], and, along similar lines, an Unscented Kalman Filter was exploited in [START_REF] Omari | Metric Visual-Inertial Navigation System Using Single Optical Flow Feature[END_REF]. All these strategies have the considerable advantage of being 'aware', to some extent, of the measurement and process noise (when modeled as Gaussian distributions). On the other hand, they require a certain level of approximation of the system dynamics which may affect the estimation performance. Other approaches exploiting tools from (deterministic) nonlinear observation can instead be found in [START_REF] Dixon | Range Identification for Perspective Vision Systems[END_REF]- [START_REF] Martinelli | Vision and IMU data fusion: Closed-form solutions for attitude, speed, absolute scale, and bias determination[END_REF] and references therein, while [START_REF] Dani | Single Camera Structure and Motion[END_REF] has tackled the more challenging problem of structure and motion estimation, i.e., how to simultaneously recover the missing structure along with the (partially) unknown camera velocity. This class of methods does not typically involve any linearization of the system dynamics and allows for some formal characterization of the estimation error convergence. However, presence of noise is not explicitly taken into account, with the filter design being developed in a fully deterministic setting. A recent experimental comparison of a EKF solution versus a deterministic nonlinear filter in the context of SfM for a quadrotor UAV can also be found in [START_REF] Grabe | A Comparison of Scale Estimation Schemes for a Quadrotor UAV based on Optical Flow and IMU Measurements[END_REF].

While all these works study the general issue of structure estimation from motion in different contexts, much less attention has been devoted to the problem of actively imposing a desired (e.g., optimized) convergence behavior to a SfM estimation task by acting on the motion imposed to the camera and on the employed estimation gains. For instance, in [START_REF] Chaumette | Structure from controlled motion[END_REF] an active strategy for minimizing the effects of image noise and discretization errors was proposed and experimentally tested, but without the aim of also imposing a desired estimation transient response. In [START_REF] Davison | Mobile Robot Localisation Using Active Vision[END_REF], the problem of actively selecting which features to track for improving the indoor localization of a wheeled mobile robot is successfully addressed; however, no attempt is made to actively shape the robot motion so as to optimize the SfM convergence (the robot navigates in an 'uninformed' way w.r.t. the estimation task). In [START_REF] Achtelik | Path Planning for Motion Dependent State Estimation on Micro Aerial Vehicles[END_REF] an EKFbased SfM estimation scheme for a UAV is integrated with a path planning strategy aiming at minimizing the covariance matrix of the estimated states at the end of the motion. Nevertheless, one needs to assume full pre-knowledge of the surrounding environment (e.g., obstacles) so as to numerically propagate the EKF filter along all the edges of a randomly constructed roadmap (the method is, thus, only amenable for an off-line/planning use).

With respect to this state-of-the-art, this paper then tackles the problem of designing an online and active algorithm for structure from known and controlled motion, i.e., assuming that the camera velocity can be measured and actively modified by the robot actuators (as it is often the case in robotic applications). The active component of the scheme makes it possible to impose an estimation error transient response equivalent to that of a reference linear second-order system with desired poles. The developments build upon the theoretical framework presented in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF]: in a nonlinear context, the observability properties of the states under consideration are not (in general) time-invariant but may depend on the current state and on the current inputs applied to the system. It is then possible to simultaneously act on the estimation gains and system inputs (i.e., the camera velocity for SfM) in order to optimize the observation process.

The methodology proposed in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] can be applied to all those systems in which an invertible function of the unknown states can appear linearly in the system dynamics, as it is indeed the case for SfM problems. We then exploit this fact and propose three concrete active SfM applications: 3D structure estimation for (i) a point feature, (ii) a spherical target and, (iii) a cylindrical target. The estimation of the depth of a point feature has already been well studied in the past literature (see, e.g., [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF], [START_REF] Dixon | Range Identification for Perspective Vision Systems[END_REF], [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], [START_REF] Morbidi | Observer design via immersion and invariance for vision-based leader-follower formation control[END_REF], [START_REF] Mahony | Adaptive depth estimation in image based visual servo control of dynamic systems[END_REF], [START_REF] Corke | Spherical Image-Based Visual Servo and Structure Estimation[END_REF]) although by never considering the active perspective taken in this work. On the other hand, the machinery proposed for the spherical and cylindrical objects represents a novel contribution also in terms of the chosen parameterization. Indeed, we show that a suitable transformation of the observed features allows to express the 3D sphere/cylinder structures in terms of image measurements and of only one unknown constant parameter (the sphere/cylinder radius) rather than the classical (and timevarying) scaled orientation of the limb surface in the camera frame. This, of course, significantly simplifies the SfM task.

We conclude by highlighting that the ability of both characterizing and optimizing the transient response of the estimation error brings several added values compared to more classical inactive estimation strategies: for instance, it allows obtaining the 'best' estimation error convergence when subject to realworld constraints such as limited camera velocity or upper bounds on the estimation gains due to noise, discretization, or other typical non-idealities. Furthermore, from a more theoretical perspective, the proposed methodology can also be used to get insights into the optimal camera trajectories needed to estimate the scene structure for particular classes of SfM problems (e.g., when dealing with point features or specific 3D geometrical primitives). Finally, we note the many similarities between the SfM approach adopted in this work and the notion of "sensor-based" or "ego-centric" Visual SLAM, see, e.g., [START_REF] Guerreiro | Globally Asymptotically Stable Sensor-Based Simultaneous Localization and Mapping[END_REF] for a recent overview. In both cases, a robot/camera builds a 3D model of the environment in its own body/sensor frame via a filtering technique: an EKF in [START_REF] Guerreiro | Globally Asymptotically Stable Sensor-Based Simultaneous Localization and Mapping[END_REF] and similar works, and the deterministic filter (but with a fully characterized and actively optimizable transient response) derived from [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] in this paper.

The rest of the paper is organized as follows: Sect. II reviews the SfM problem in the context of nonlinear state observation and briefly summarizes the methodology developed in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] for actively imposing a desired transient behavior to the estimation error. Section III then focuses on the three SfM problems considered in this work. Subsequently, Sect. IV reports the corresponding simulative and experimental results obtained with a manipulator equipped with an eye-in-hand camera. Finally, Sect. V concludes the paper and discusses some future directions.

II. PRELIMINARIES

In this section, we briefly summarize the active estimation framework originally proposed in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF]. This is then applied to SfM case studies discussed in the next Sect. III.

A. A nonlinear observation scheme

Let (s , χ) ∈ R m+p be the state of a dynamical system in the form

ṡ = f m (s, u, t) + Ω T (t)χ χ = f u (s, χ, u, t) (1)
where s ∈ R m and χ ∈ R p represent, respectively, a measurable and unmeasurable component of x, and u ∈ R v is the system input vector. In formulation (1) vector χ is required to appear linearly in the dynamics of s (first equation). Furthermore, matrix Ω(t) ∈ R p×m and vectors f m (•) ∈ R m and f u (•) ∈ R p are assumed to be generic but known and sufficiently smooth functions w.r.t. their arguments which are all available apart from the unknown value of χ in f u (•). SfM problems can be recast to formulation (1) by taking s as a set of visual features measured in the image, u = (v, ω) as the camera linear/angular velocity in camera frame, and χ as a suitable (and locally invertible) function of the unknown structure of the scene to be estimated. For instance, in the point feature case, χ can be taken as the inverse of the feature depth [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], and, for image moments of planar scenes, χ can be taken as the normal vector of the observed plane scaled by its distance from the camera optical center [START_REF] Giordano | 3D structure identification from image moments[END_REF]. Furthermore, in SfM one has Ω(t) = Ω(s(t), v(t)) with, in particular, Ω(s, 0) ≡ 0: the camera linear velocity v(t) plays a key role for the resolution of SfM problems1 .

For a system in form (1), a possible estimation scheme can be devised as follows [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF]: let (ŝ , χ) ∈ R m+p be the estimated state, ξ = s -ŝ, z = χ -χ, e = (ξ, z) be the total error vector, and consider the following observer

ṡ = f m (s, u, t) + Ω T (t) χ + Hξ χ = f u (s, χ, u, t) + ΛΩ(t)Qξ (2)
where H > 0, Λ = Λ T > 0 and Q = Q T > 0 are positive definite gain matrices. Note that observer (2) is function of only measured/known quantities, with in particular a feedback action on the measurable error component ξ. The corresponding estimation error dynamics is then given by

   ξ = -Hξ + Ω T (t) z ż = -ΛΩ (t) Qξ + (f u (s, χ, u) -f u (s, χ, u)) = -ΛΩ (t) Qξ + g(e, t) (3
) with g(e, t) being a 'perturbation term' vanishing w.r.t. the error vector e, i.e., such that g(0, t) = 0, ∀t. The origin of (3) can be proven to be locally exponentially stable if and only if (iff) the following Persistency of Excitation (PE) condition holds

t+T t Ω (τ ) Ω T (τ ) dτ ≥ γI p > 0, ∀t ≥ t 0 , (4)
for some T > 0 and γ > 0, with I n representing the n × n identity matrix2 .

Remark II.1. We note that the local stability properties of the error dynamics (3) are due to the perturbation term g(e, t) which affects an otherwise globally exponentially stable error system. Indeed, in the special case χ = 0 (unknown but constant parameters), one has g(e, t) ≡ 0 and global exponential convergence for the error system (3). This is, for instance, the case of the structure estimation problems for spherical and cylindrical objects considered in Sects. III-B and III-C. We stress, however, that the estimation scheme [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF] is not restricted to this particular situation but can be applied (with, in this case, only local convergence guarantees) to the more general case of state observation problems in which the unknown χ is subject to a non-negligible dynamics as in [START_REF] Strasdat | Visual SLAM: Why filter?[END_REF]. The depth estimation for a point feature discussed in the following Sect. III-A falls in this class.

The PE condition (4) plays the role of an observability criterium: convergence of the estimation error e(t) → 0 is possible iff the square matrix Ω (t) Ω T (t) ∈ R p×p remains full rank in the integral sense of (4). We note that if m ≥ p, that is, if the number of independent measurements s is larger or equal to the number of estimated quantities χ, then it is possible to instantaneously satisfy (4) by enforcing

Ω (t) Ω T (t) ≥ γ T I p , ∀t. (5)
In the rest of this work we will only consider this (more restrictive) observability condition.

B. An active estimation strategy

As clear from (4) and (5) (some measure of) the norm of matrix ΩΩ T determines the convergence properties of the error system (3). Furthermore, since in the SfM case it is Ω(t) = Ω(s(t), v(t)), it is meaningful to study how to optimize the camera linear velocity v in order to affect matrix ΩΩ T and, as a consequence, to shape the transient response of the error vector e(t). The active strategy developed in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] and summarized hereafter shows how to achieve this goal.

Let U ΣV T = Ω be the singular value decomposition of matrix Ω, where Σ = [S 0], S = diag(σ i ) ∈ R p×p , and 0 ≤ σ 1 ≤ . . . ≤ σ p are the p singular values of Ω. Let also Q = αI m and Λ = βI p , with α > 0, β > 0 (scalar gain matrices). By designing the gain matrix H in (2) as

H = V D 1 0 0 D 2 V T (6)
with

D 1 ∈ R p×p > 0, D 2 ∈ R (m-p)×(m-p) > 0, it is possible to show that, under the change of coordinates η = 1 √ αβ S -1 U T z (7)
and in the approximation S -1 U T ≈ const, the behavior of vector η (and hence of the estimation error z = χ -χ) is governed by the following linear (and almost diagonal)

dynamics η = (Π -D 1 ) η -αβS 2 η. (8)
System (8) can be interpreted as a (unit-)mass-spring-damper system with diagonal stiffness matrix αβS 2 and damping matrix D 1 , together with an additional 'perturbing' term Π whose full expression can be found in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF].

The convergence rate of (8) is then related to its slowest mode dictated by the 'stiffness value' αβσ 2 1 , with σ 2 1 being the smallest eigenvalue of the square matrix ΩΩ T . Therefore, for the sake of imposing a desired transient response to vector η(t) (i.e., to the estimation error z(t) = χ(t) -χ(t)), one can 'place the poles' of ( 8) by (i) shaping the damping factor D 1 in (6) (a free parameter), (ii) regulating the value of the smallest eigenvalue σ 2 1 by acting upon vector v, and (iii) suitably choosing the gain αβ (a free parameter).

For what concerns the design of matrix D 1 , we first note that, as explained in [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF], matrix Π in (8) can be regarded as a second-order perturbation term affecting the dissipative action induced by D 1 . Therefore, neglecting the effects of matrix Π and choosing D 1 = diag(c i ), c i > 0, allows obtaining a completely decoupled transient behavior for [START_REF] Metni | Visual Tracking Control of Aerial Robotic Systems with Adaptive Depth Estimation[END_REF] ηi

+ c i ηi + αβσ 2 i η i = 0, i = 1 . . . p. (9)
One can then take, for instance, c i = c * i = 2 √ αβσ i in order to impose a critically damped evolution to the estimation error (coincident eigenvalues for (9)).

As for the regulation of σ 1 (t), being Ω = Ω(s, v), it is

(σ 2 i ) = J v,i v + J s,i ṡ (10)
where the Jacobian matrices J v,i ∈ R 1×v and J s,i ∈ R 1×n can be computed in closed form, see [START_REF] Spica | A Framework for Active Estimation: Application to Structure from Motion[END_REF] for all the details. By inverting the differential mapping [START_REF] Giordano | 3D structure identification from image moments[END_REF], vector v can then be exploited so as to, e.g., asymptotically enforce σ 2 1 (t) → σ 2 1,des for some desired value σ 2 1,des > 0. We note that ensuring σ 2 1 (t) → σ 2 1,des > 0 also automatically satisfies the observability condition [START_REF] Civera | 1-Point RANSAC for EKF Filtering: Application to Real-Time Structure from Motion and Visual Odometry[END_REF].

Finally, the following considerations hold for the choice of gain αβ in [START_REF] Metni | Visual Tracking Control of Aerial Robotic Systems with Adaptive Depth Estimation[END_REF]. In the SfM context, the norm of matrix ΩΩ T is strongly related to the norm of the camera linear velocity v. Roughly speaking, the 'faster' the motion (∼ larger v ), the 'larger' the value of σ 2 1 (∼ larger ΩΩ T ). Therefore, in order to maximize the estimation convergence speed of (8) (dictated by αβσ 2 1 ), one can equivalently (i) travel at a larger speed v for a given gain αβ, or (ii) increase the gain αβ for a given v . While increasing the gain αβ may always appear more convenient in terms of reduced control effort, practical issues such as noise, discretization or quantization errors, may impose an upper limit on the possible value of αβ, thus necessarily requiring a larger v for obtaining the desired convergence speed. Furthermore, as in all SfM problems, a v = 0 is also mandatorily required for guaranteeing σ 2 1 > 0 (a non-translating camera cannot estimate the scene structure).

Remark II.2. We note that the proposed strategy is an active one since, in the general case, inversion of [START_REF] Giordano | 3D structure identification from image moments[END_REF] will result in a camera linear velocity v optimized as a function of the system measured state s in order to attain a desired σ 2 1,des over time. We also highlight the role played by the gain matrix H weighting the feedback term in observer [START_REF] Matthies | Kalman filter-based algorithms for estimating depth from image sequences[END_REF]: the proposed machinery in fact relies on a suitable state-dependent design of the damping matrix H via the choice of

D 1 = diag(c * i ) with c * i = 2 √ αβσ i .
For the interested reader, this statedependent design is conceptually equivalent to the shaping of the damping factor adopted in the context of impedance control for robot manipulator arms, see, e.g., [START_REF] Albu-Schaffer | Cartesian impedance control of redundant robots: recent results with the DLRlight-weight-arms[END_REF].

Remark II.3. We also note that, in general, it is not possible to fully compensate for the term J s,i ṡ when inverting [START_REF] Giordano | 3D structure identification from image moments[END_REF]. Indeed, the expression in ( 1) implies a direct dependence of ṡ from the unmeasurable χ, so that an exact evaluation of ṡ is not obtainable in practice. A possible solution could be to use an approximation ŝ of ṡ obtained by evaluating f u (•) on the current estimate χ. Another simple workaround is, however, to just enforce ṡ ≃ 0 by imposing a constraint on the camera motion. A combination of both strategies is, of course, also possible. The next sections will present some examples in this sense.

Remark II.4. It can be finally shown that, in the special situation p = 1 (only one quantity to be estimated), if σ 1 (t) ≡ const then S -1 U T ≡ const in [START_REF] Dixon | Range Identification for Perspective Vision Systems[END_REF] and matrix Π has no disturbing effects on [START_REF] Metni | Visual Tracking Control of Aerial Robotic Systems with Adaptive Depth Estimation[END_REF]. Therefore, in this case it is always possible to exactly enforce the ideal estimation error dynamics ( 9) by just keeping Ω(t) 2 = σ 2 1 (t) = const during the camera motion. This situation will apply to all the case studies discussed in the rest of the paper.

III. APPLICATIONS TO STRUCTURE FROM MOTION

In this section we illustrate the application of the proposed active estimation framework to three concrete SfM problems: (i) estimation of the 3D coordinates of a point feature, (ii) estimation of the 3D position and radius of a spherical target, and (iii) estimation of the 3D position and radius of a cylindrical target.

In the point feature case, the effects of the adopted projection model on the estimation convergence are also explicitly considered by discussing the differences between the two popular choices of planar and spherical projection models. For the spherical and cylindrical targets, we instead propose two novel minimal parameterizations that allow to express the sphere/cylinder 3D structures in terms of measured visual features and of a single unknown parameter (the sphere/cylinder radius). This allows, in all three cases, to reduce the SfM task to the estimation of a single unknown quantity (point feature depth or sphere/cylinder radius), thus satisfying the requirements of Remark II.4 for exactly imposing the ideal dynamics [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF] to the estimation error.

A. Depth estimation for a point feature 1) Planar projection model: Let p = (x, y, 1) = (X/Z, Y /Z, 1) ∈ R 3 be the perspective projection of a 3D point P = (X, Y, Z) onto the image plane of a calibrated pinhole camera. As it is well know [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF], the differential relationship between the image motion of a point feature and the camera linear/angular velocity u = (v, ω) ∈ R 6 expressed in camera frame is

ẋ ẏ =    - 1 Z 0 x Z xy -1 + x 2 y 0 - 1 Z y Z 1 + y 2 -xy -x    u (11) where Z is the depth of the feature point. The dynamics of Z is Ż = 0 0 -1 -yZ xZ 0 u.
The expression in [START_REF] Morbidi | Observer design via immersion and invariance for vision-based leader-follower formation control[END_REF] is not linear in Z but it is linear in 1/Z. Therefore, by defining s = (x, y) ∈ R 2 and χ = 1/Z, with then m = 2 and p = 1, we obtain for (1)

         f m (s, u, t) = xy -1 + x 2 y 1 + y 2 -xy -x ω Ω (s, v) = xv z -v x yv z -v y f u (s, χ, u, t) = v z χ 2 + (yω x -xω y ) χ , (12)
with the perturbation term g(e, t) in (3) taking the expression

g(e, t) = v z χ 2 -χ2 + (yω x -xω y ) z, (13)
so that g(0, t) = 0 as expected. Note that, once χ has been estimated, one can obviously retrieve the 3D position of the point feature as P = p/χ.

In the point feature case matrix ΩΩ T reduces to its single eigenvalue which, for a planar projection model, takes the expression

σ 2 1 = Ω 2 = (xv z -v x ) 2 + (yv z -v y ) 2 . ( 14
)
Furthermore, using [START_REF] Martinelli | Vision and IMU data fusion: Closed-form solutions for attitude, speed, absolute scale, and bias determination[END_REF], the Jacobian J v,1 in ( 10) is given by

J v,1 = 2   v x -xv z v y -yv z (xv z -v x ) x + (yv z -v y ) y   T . (15)
Since σ 2 1 does not depend on ω, it is then possible to freely exploit the camera angular velocity for fulfilling additional goals of interest without interfering with the regulation of σ 2 1 (t) (only affected by v). For instance, as in [START_REF] Chaumette | Structure from controlled motion[END_REF], one can use ω for keeping s ≃ const so as to make the effects of ṡ negligible when inverting (10) w.r.t. v, see Remark II.3.

We now note that σ 2 1 in ( 14) depends on both the camera linear velocity v and on the location p of the feature point on the image plane. Since the value of σ 2 1 directly affects the convergence speed of the estimation error, it is interesting to study what conditions on p and v result in the largest possible σ 2 1 (i.e., the fastest possible convergence for a given gain αβ). Letting e 3 = (0, 0, 1) being the camera optical axis, it is (by inspection)

Ω T 0 = [e 3 ] × [p] × v
where [v 1 ] × is the skew-symmetric matrix representing the cross product operator for 3D vectors (i.e.,

[v 1 ] × v 2 = v 1 × v 2 )
. Therefore,

σ 2 1 = Ω 0 Ω T 0 = [e 3 ] × [p] × v 2 = p 2 v 2 sin 2 (θ p,v ) sin 2 θ e3,[p] × v
where θ p,v and θ e3,[p] × v represent the angles between vectors (p, v) and vectors (e 3 , [p] × v), respectively. The maximum attainable value for σ 2 1 is then

σ 2 max = max p, v σ 2 1 = p 2 v 2 . ( 16
)
This maximum is obtained when the camera linear velocity v is such that p ⊥ v and e 3 ⊥ [p] × v, i.e., rearranging in matrix form

p T e T 3 [p] × v = x y 1 -y x 0 v = 0. ( 17
)
If p = e 3 (point feature not at the center of the image plane), system ( 17) has (full) rank 2 and admits the unique solution (up to a scalar factor)

v = δ [p] 2 × e 3 , δ ∈ R.
This requires the linear velocity v to be orthogonal to p and to lie on the plane defined by vectors p and e 3 (i.e., v must belong to a straight line as shown in Fig. 1a).

If p = e 3 (point feature at the center of the image plane), system [START_REF] Chaumette | Structure from controlled motion[END_REF] loses rank and any v ⊥ e 3 is a valid solution, see Fig. 1b.

It is then possible to draw the following conclusions: for a given norm of the linear velocity v (i.e., the amount of 'control effort'), system [START_REF] Chaumette | Structure from controlled motion[END_REF] determines the direction of v resulting in σ 2 1 = σ 2 max (maximization of σ 2 1 ). These conditions are summarized in Figs. 1a and1b. The value of σ 2 max is, however, also a function of the feature point location p which can be arbitrarily positioned on the image plane. In particular, σ 2 max = v 2 for p = e 3 and σ

2 max = p 2 v 2 > v 2 ∀ p = e 3 , with lim p →∞ σ 2 max (p) = ∞.
The value of p (distance of the point feature from the image center) thus acts as an amplification factor for σ 2 max . Therefore, 1) the smallest σ 2 max (i.e., the slowest 'optimal' convergence for the depth estimation error) is obtained for the smallest value of p , i.e., when p = e 3 =⇒ p = 1 (feature point at the center of the image plane). It is worth noting that in this case v z = 0 (from the condition v ⊥ p) and

σ 2 max = v 2 = v 2
x +v 2 y : the camera moves on the surface of a sphere with a constant radius (depth) pointing at the feature point. Also, being in this case χ = Ż/Z 2 = 0, one has g(e, t) ≡ 0 and global convergence for the estimation error (see Remark II.1); 2) the largest σ 2 max (i.e., the fastest 'optimal' convergence for the depth estimation error) is obtained for the largest possible value of p . In the usual case of a rectangular image plane centered at the origin, this translates into keeping the feature point positioned at one of the four image corners. However, compared with the previous case, this results in a g(e, t) = 0 and only local convergence for the estimation error. 2) Spherical projection model: We now develop the depth estimation machinery for the spherical projection model. In this case, the following quantity is taken as visual feature measured on the image plane

s = p p = P P ∈ S 2 ,
where S 2 represents the unit sphere and, as well-known [START_REF] Hamel | Visual servoing of an under-actuated rigid body system: An image based approach[END_REF],

ṡ = 1 P ss T -I 3 [s] × u, and
d dt 1 P = - 1 P 2 d P dt = - s T Ṗ P 2 = s T v P 2 .
Hence by taking χ = 1/ P one obtains for (1)

     f m (s, u) = [s] × ω Ω (s, v) = -v T I 3 -ss T f u (s, χ, u) = χ 2 s T v (18)
with m = 3, p = 1, and g(e, t) = (χ 2 -χ2 )s T v for the perturbation term in (3). We note that, although in this case m = 3, vector s is subject to the constraint s = 1, thus resulting in only two independent measurements (as in the previous case of planar projection). Moreover, from the estimated χ one can easily retrieve P = s/χ.

For the spherical projection model, the eigenvalue determining the convergence of the estimation error is

σ 2 1 = ΩΩ T = v T v -(s T v) 2 , with thus J v,1 = 2v T (I 3 -ss T ). (19)
As before, σ 2 1 does not depend on ω which can then be exploited to fulfil any additional task of interest (e.g., keeping s ≃ const during motion).

As for the conditions on s and v that yield maximization of σ 2 1 , one clearly has

σ 2 1 = σ 2 max = max s, v σ 2 1 = v 2 (20)
iff s T v = 0 (linear velocity orthogonal to the projection ray passing through P ). We also note that, in this case, one has χ = 0 and g(e, t) ≡ 0 (constant unknown state and global convergence for the estimation error) regardless of the location of s on the image plane.

3) Comparison between planar and spherical projection models: for a spherical projection model, maximization of the eigenvalue σ 2 1 imposes only one condition for the linear velocity v (s T v = 0). When this condition is met, one has σ 2 1 = σ 2 max = v 2 and global convergence for the estimation error whatever the location of the feature point s. This is equivalent to what was obtained for the planar projection case in the special situation p = e 3 (indeed the two projection models coincide for p = s = e 3 ). However, with a spherical projection model one also loses the possibility to increase the estimation convergence rate by suitably positioning the point feature s on the image plane (since in this case σ 2 max does not depend on s).

It is then worth noting the complementarity of the two cases: for a given v , and provided the optimal condition p T v = 0 is satisfied, the planar projection allows obtaining a faster error convergence at the price of local stability (increase of the perturbation g) by suitably positioning s = (x, y) (the larger s the faster the convergence). The spherical projection guarantees global error convergence for any location of the feature point but at the price of being always subject to the same convergence rate only function of the control effort v .

B. Structure estimation for a spherical target

We now detail the application of the proposed estimation machinery to the case of a spherical target. Consider a sphere O s of radius R and let P 0 = (X 0 , Y 0 , Z 0 ) be the coordinates of its center in the camera frame. Let also

L : n T X + d = 0
represent the planar limb surface associated to the sphere in the camera frame, where X ∈ R 3 is any 3D point on the plane, n ∈ S 2 is the plane unit normal vector and d ∈ R the plane distance to the camera center [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]. Figure 2 shows the quantities of interest.

The depth Z of any point X lying on L can be expressed in terms of its normalized image coordinates p = (x, y, 1) as

1 Z = X 0 K x + Y 0 K y + Z 0 K = χ T p, (21)
where K = P T 0 P 0 -R 2 and χ = P 0 /K = -n/d ∈ R 3 represent unmeasurable quantities (analogously to Z for the point feature case), see [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] for all the details. The interaction matrix of a generic (i, j)-th order moment m ij evaluated on the image of O s depends linearly on χ, see [START_REF] Giordano | 3D structure identification from image moments[END_REF], [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]. Therefore, a first possibility to retrieve the sphere 3D parameters (P 0 , R) would be to implement the estimation scheme (2) with s being a suitable collection of image moments (e.g., area and barycenter). It is in fact possible to show that (see Appendix A)

χ = - v K -[ω] × χ + 2χχ T v
and that K can be expressed in terms of image moments and of vector χ itself, so that, having estimated χ, one can consequently retrieve P 0 = χK and R = P T 0 P 0 -K. Although conceptually valid, this solution requires the concurrent estimation of three time-varying quantities (vector χ(t)). On the other hand, inspired by [START_REF] Fomena | Improvements on visual servoing from spherical targets using a spherical projection model[END_REF], we now describe a novel representation of the sphere projection on the image plane that allows to reformulate the structure estimation task in terms of a single unknown constant parameter, i.e., the sphere radius R.

To this end, define vector s = (s x , s y , s z ) ∈ R 3 as

                 s x = x g s z a 2 1 s y = y g s z a 2 1 s z = 1 + a 2 1 a 2 1 , (22)
where (x g , y g , n 20 , n 11 , n 02 ) represent the barycenter and normalized centered moments of order 2 measured from the elliptical projection of the sphere O s on the image plane, and a 1 is the minor axis of the observed ellipse with [27]

a 2 1 = 2 n 20 + n 02 -(n 20 -n 02 ) 2 + 4n 11 . (23)
We thus note that vector s can be directly evaluated in terms of measured image quantities.

From [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF], [START_REF] Fomena | Improvements on visual servoing from spherical targets using a spherical projection model[END_REF] one also has

x g = X 0 Z 0 Z 2 0 -R 2 , y g = Y 0 Z 0 Z 2 0 -R 2 , a 2 1 = R 2 Z 2 0 -R 2 (24
) which, when plugged in [START_REF] Corke | Spherical Image-Based Visual Servo and Structure Estimation[END_REF] and [START_REF] Guerreiro | Globally Asymptotically Stable Sensor-Based Simultaneous Localization and Mapping[END_REF], result in the equivalent expression s = P 0 /R. Since vector s can be computed from image measurements as in [START_REF] Corke | Spherical Image-Based Visual Servo and Structure Estimation[END_REF], estimation of the (unknown) sphere radius R allows to recover the 3D sphere center as P 0 = sR.

Exploiting now the results of [START_REF] Fomena | Improvements on visual servoing from spherical targets using a spherical projection model[END_REF], it is possible to show that

ṡ = - 1 R I 3 [s] × u. (25)
Since ( 25) is linear in 1/R, we can define χ = 1/R, with then m = 3 and p = 1, and obtain for (1) and ( 3)

           f m (s, u) = [s] × ω Ω (s, u) = -v T f u (s, χ, u) = 0 g(e, t) = 0 . ( 26
)
We note that in this case it is always possible to obtain global convergence for the estimation error since χ = 0 and therefore g(e, t) = 0 by construction (see Remark II.1). Furthermore, matrix ΩΩ T reduces again to its single eigenvalue σ 2 1 = v 2 and, if σ 2 1 (t) ≡ const > 0, the 'ideal' estimation error dynamics (9) can be exactly obtained. One also has Ω = Ω(v) and J v,1 = 2v T .

We finally note the following facts: first, contrarily to the previous cases, here ṡ has no effect on the regulation of σ 2 1 which is only function of the camera linear velocity v. It is then of course still possible to freely exploit the camera angular velocity ω for, e.g., keeping the sphere at the center of the image by regulating (s x , s y ) to zero. Second, we note the strong similarities with the previous optimal results obtained for a point feature under a spherical projection model (σ 2 max in (20)): in both cases the maximum estimation convergence rate for a given v does not depend on the position of the observed object on the image plane.

C. Structure estimation for a cylindrical target

We now finally consider the case of SfM for a 3D cylindrical object. A cylinder O c can be described by its radius R > 0 and by its main axis a ∈ S 2 passing through a 3D point P 0 = (X 0 , Y 0 , Z 0 ), with a = 1 and, w.l.o.g., a T P 0 = 0 (P 0 can be chosen as the closest point on a to the origin of the camera frame [START_REF] Chaumette | Structure from controlled motion[END_REF]). Moreover, analogously to the sphere, a cylinder is also associated with a planar limb surface L such that [START_REF] Mahony | Adaptive depth estimation in image based visual servo control of dynamic systems[END_REF] holds for any point on L with projection p = (x, y, 1). Therefore, a possibility is to estimate the three unknown parameters of the limb plane L (vector χ) by exploiting (at least) three image measurements, see [START_REF] Chaumette | Structure from controlled motion[END_REF] and Appendix B for some details in this sense. However, following the previous developments, we now propose a novel representation of the cylinder projection on the image plane which, again, allows to obtain the cylinder parameters (P 0 , a, R) in terms of image measurements and of the unknown but constant cylinder radius R which, therefore, represents the only quantity to be estimated.

Let (ρ 1 , θ 1 ) and (ρ 2 , θ 2 ) be the (measured) distance/angle parameters of the two straight lines resulting from the projection of the cylinder on the image plane, and

n 1 = (cos θ 1 , sin θ 1 , -ρ 1 ), n 2 = (cos θ 2 , sin θ 2 , -ρ 2 )
(27) be the normal vectors to the two planes passing through the origin of the camera frame and the two above-mentioned projected lines 3 . Figure 3 gives a graphical representation of the quantities of interest. Note that vectors n 1 and n 2 3 The two planes are therefore tangent to the surface of the cylinder. can be directly evaluated from image measurements (the line parameters). We then define vector s ∈ R 3 as

s = ∆ ∆ 2 (28)
with

∆ = 1 2 n 1 n 1 + n 2 n 2 . ( 29
)
Vector s is, thus, also directly obtainable in terms of image quantities.

We now note that, from [START_REF] Chaumette | Visual servoing using image features defined upon geometrical primitives[END_REF], an equivalent expression for vectors n 1 , n 2 in terms of the cylinder 3D geometry can be obtained as

n 1 = 1 N 1          R X 0 √ K -α R Y 0 √ K -β R Z 0 √ K -γ          , n 2 = 1 N 2        R X0 √ K + α R Y 0 √ K + β R Z 0 √ K + γ        (30)
with

                     K = P T 0 P 0 -R 2 (α, β, γ) = [P 0 ] × a N 1 = R X0 √ K -α 2 + R Y0 √ K -β 2 N 2 = R X0 √ K + α 2 + R Y0 √ K -β 2 , (31) thus yielding
                 n 1 n 1 = 1 P T 0 P 0   RX 0 -α √ K RY 0 -β √ K RZ 0 -γ √ K   n 2 n 2 = 1 P T 0 P 0   RX 0 + α √ K RY 0 + β √ K RZ 0 + γ √ K   . (32)
Plugging ( 32) in ( 29) results in the equivalent expression

∆ = R 2
P T 0 P 0 s which, using [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], finally yields the following relationship between image quantities and cylinder 3D structure

s = ∆ ∆ 2 = P 0 R . (33)
As for the cylinder axis a, exploiting (30) one has

[n 2 ] × n 1 = 2R N 1 N 2 √ K   Z 0 β -Y 0 γ X 0 γ -Z 0 α Y 0 α -X 0 β   = 2R N 1 N 2 √ K   α β γ   × P 0 = 2R N 1 N 2 √ K [P 0 ] × a × P 0 = 2RP T 0 P 0 N 1 N 2 √ K a (34)
where in the last step the property a T P 0 = 0 was used. Since a = 1, from (34) it is

a = [n 2 ] × n 1 [n 2 ] × n 1 . ( 35
)
The cylinder axis a can then be directly obtained in terms of only measured quantities. We now note that, as in the sphere case, the only unknown left is the cylinder radius R: once known, the cylinder 3D structure can be fully recovered from image measurements as P 0 = Rs from (33) and a from (35).

An estimation scheme for R can be obtained exploiting the following differential relationship whose derivation is given in Appendix

C ṡ = - 1 R I 3 -aa T [s] × u. ( 36
)
Note the similarity of (36) with [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] for the sphere case. Being (36) linear in 1/R, one can then apply observer (2) by choosing χ = 1/R with m = 3 and p = 1, and obtaining

           f m (s, u) = [s] × ω Ω (s, u) = -v T I 3 -aa T f u (s, χ, u) = 0 g(e, t) = 0 . ( 37
)
Note how, again, being χ = 0 it is g(e, t) = 0 (global convergence for the error system (3) as in the sphere case). Matrix ΩΩ T reduces to its single eigenvalue

σ 2 1 = ΩΩ T = v 2 -(a T v) 2 . ( 38
)
It is worth comparing (38) with the result obtained for the sphere (σ 2 1 = v 2 ). In the cylinder case, the convergence rate of the estimation error is affected by both the norm and the direction of the linear velocity v. In particular, for a given v = const, the maximum value for σ 2 1 is obtained when v has a null component along the cylinder axis a (a T v = 0) with, in this case, σ 2 1 = σ 2 max = v 2 . Intuitively, any camera motion along the cylinder axis does not provide any useful information to the estimation task. Furthermore, as in all previous cases, keeping a σ 2 1 (t) = const allows to exactly enforce the ideal estimation error dynamics [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], see Remark II.4.

Finally, from (38) one has

(σ 2 1 ) = J v,1 v + J a,1 ȧ = J v,1 v + J a,1 [a] × ω (39)
with J v,1 = 2v T I 3aa T and J a,1 = 2v T av T . Although (39) also depends on the angular velocity ω, it is possible to fully compensate for the effects of J a,1 [a] × ω (a known quantity) when inverting (39) w.r.t. v as discussed in Sect. IV-D. Therefore, one can act on v to regulate the value of σ 2 1 (t) and, at the same time and in a decoupled way, exploit the camera angular velocity ω for implementing additional tasks of interest such as keeping the cylinder axis a at the center of the image plane by enforcing (s x , s y ) = 0.

IV. SIMULATION AND EXPERIMENTAL RESULTS

In this section we show some experimental and simulation results meant to validate the theoretical developments of the previous sections. The experiments were run by employing a greyscale camera with a resolution of 640 × 480 px and a framerate of 30 fps. The camera was mounted on the endeffector of a 6-dofs Gantry robot commanded in velocity at a frequency of 100 Hz. All the image processing and feature tracking were implemented via the open-source ViSP library [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. Some snapshots of the three experiments are shown in Fig. 4 where the result of the image processing is highlighted in red.

A video of the reported experiments (including the three cases of point, sphere and cylinder) is also attached to the paper for the reader's convenience.

A. Comparison of planar and spherical projection models

We start by comparing via simulation results the effects of adopting a planar and spherical projection model for the depth estimation of a point feature as extensively discussed in Sect. III-A1 and Sect. III-A2. We considered three cases differing for the location on the image plane at which the point feature was (purposely) kept exploiting the camera angular velocity ω:

1) in case I the point feature was kept at the center of the image plane (red line in the following plots); 2) in case II the point feature was kept at one of the corners of an image plane with the same size of the camera used in the experiments (green line in the following plots); 3) in case III the point feature was kept at one of the corners of an image plane with a size five times larger than case II (blue line in the following plots). In all cases, a constant camera velocity v(t) ≡ v(t 0 ) = const was kept during motion, with the initial condition v(t 0 ) chosen so as to comply with the optimality conditions discussed in Sects. III-A1 and III-A2 for letting σ 2 1 = σ 2 max (e.g., with v(t 0 ) being a solution of [START_REF] Chaumette | Structure from controlled motion[END_REF] in the planar projection case).

Figure 5a shows the behavior of z(t) for the three cases when using a planar projection model. We can then note how the convergence rate of the estimation error increases from case I (slowest convergence) to case III (fastest convergence) as predicted by the theory (for the same v a larger p results in a larger σ 2 max ). Similarly, Fig. 5b reports the behavior of σ 2 1 (t) for the three cases: as expected, σ 2 1 (t) results largest for case III. Note also how σ 2 1 (t) for case II (green line) is only slightly larger than case I (red line). This is due to relatively small size of the image plane of case II whose dimensions were set as those of the real camera used for the experiments. Finally, Fig. 5c shows the behavior of the perturbation term g(e, t) in the three cases: here, one can verify how g = 0 for case I, with then an increasing |g| for cases II and III. Indeed, as discussed in Sect. III-A1, the 'amplification' effect on σ 2 max obtained by increasing p comes at the price of an increased magnitude of the perturbation g. This is also evident in Fig. 5a where the ideal response of ( 9) is plotted with dashed lines for the three considered cases. We can thus note how z(t) in case I presents a perfect match with its corresponding ideal response, with then an increasing (albeit very limited) mismatch in the other two cases due to the increased effect of the perturbation g.

As for the spherical projection model, Fig. 5d reports the behavior of the estimation error z(t) for the three cases under consideration, together with the ideal response [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF]. Here, the symbol z s (t) is used to denote the estimation error in the spherical projection case in order to distinguish it from the error obtained with the planar projection model. All the plots result perfectly superimposed as expected from the analysis of Sect. III-A2. Indeed, in the spherical projection case, σ 2 max = v 2 regardless of the location of p and g(t) ≡ 0. However, absence of perturbation terms is obtained at the expense of the convergence rate of z s (t), which indeed results slower or equal to that of z(t) in the planar projection case. This is shown in Fig. 5e where the behavior of z(t) -z s (t) is reported for the three cases. We can then note how z(t) -z s (t) = 0 only in case I, as the planar and spherical models coincide when the feature point is at the center of the image plane.

These results then fully confirm the validity of the theoretical analysis reported in Sects. III-A1 and III-A2. However, we also note the marginal effects of the two projection models on the estimation performance when applied to an image plane of size comparable to that of the real camera used in our experimental setup. Therefore, in the following experimental results we will only consider the case of planar projection model.

B. Depth estimation for a point feature

We here report some experimental results for the depth estimation of a point feature under a planar projection model (Sect. III-A1). The following experiments are meant to demonstrate how the proposed active estimation framework can be exploited to select online the 'best' camera motion. As visual target, we made use of a circular white dot of 5 mm radius painted on a planar black surface and sufficiently far from the camera in order to safely consider it as a 'point feature'.

Figure 6a shows the evolution of the estimation error z(t) = 1/Z(t) -1/ Ẑ(t) for two experiments 4 in which v(t) = v 0 but with its direction being either optimized in order to maximize the estimation convergence rate (case I, red line) or kept constant so that v(t) = v 0 = const (case II, blue line). This effect was obtained by using the following control law5

v = v v 2 k 1 (κ des -κ) + k 2 I 3 - vv T v 2 J T v,1 (40)
with [START_REF] Strasdat | Visual SLAM: Why filter?[END_REF] given by [START_REF] Dani | Single Camera Structure and Motion[END_REF]. In fact, the first term in (40) enforces the constraint v(t) = v 0 (same control effort in both cases), while the second term allows to implement either case I (k 2 > 0) or case II (k 2 = 0) (maximization of σ 2 1 ) within the null-space of the first constraint. In both cases, the angular velocity ω was exploited for keeping the point feature at the center of the image plane (x, y) → (0, 0). We note that, as discussed in Sect. III-A1, when (x, y) = (0, 0) one has σ 2 max = v 2 x +v 2 y from ( 16) and σ 2 1 = σ 2 max iff v z = 0 (circular motion around the point feature). The experiments were run with the following parameters: αβ = 10 3 for gains Q and Λ, c 1 = c * 1 for D 1 in (6), v(t 0 ) = v 0 = (0.03, 0, -0.04) m/s, k 1 = 5 and k 2 = 10 4 , thus resulting in the maximum value σ 2 max = 0.0025 for the eigenvalue σ 2 1 . As clear from Fig. 6b, while in case II the camera gets closer to the point feature, the use of the active strategy of case I results in a null component of v along the projection

k 1 > 0, k 2 ≥ 0, κ = 1 2 v T v, κ des = 1 2 v T 0 v 0 , and J v,
Fig. 5: Simulation results comparing the planar and spherical projection models for the depth estimation of a point feature. The following color coding is adopted for the three considered cases: red-case I, green-case II, blue-case III. (a) behavior of the estimation error z(t) in the planar projection case (solid lines) with superimposed the corresponding ideal response (9) (dashed lines). (b) behavior of σ 2 1 (t) for the three cases with, again, the largest σ 2 1 (t) in case III. (c) behavior of the perturbation term g(e, t) for the three cases. (d) behavior of the estimation error zs(t) for the spherical projection model in the three cases. (e) behavior of z(t) -zs(t). ray of the point feature (i.e., v z = 0) and in an associated circular trajectory centered on the tracked point as predicted by the theoretical analysis of Sect. IV-A. This then allows to move faster in the 'useful' directions (while keeping the same constant v ), and, thus, to increase the value of σ 2 1 towards its theoretical maximum σ 2 max = 0.0025 (Fig. 6c), resulting in an overall faster convergence for the estimation error (Fig. 6a). Furthermore, Fig. 6a also reports the ideal response of [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF] with desired poles at σ 2 max (dashed black line). We can then note the almost perfect match with case I (solid red line): indeed, as explained in Remark II.4, imposing a σ 2 1 (t) = const allows to exactly obtain the ideal behavior governed by [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF]. It is finally worth noting the accuracy of the reconstructed depth: Fig. 6a reports two vertical dashed lines indicating, for the two cases under consideration, the times T 1 = 4.95 s and T 2 = 9.85 s at which the estimation error z(t) becomes smaller than 5 mm. We then obtained a standard deviation of approx. 0.8 and 0.3 mm evaluated on a time window of 1 s after the times T 1 and T 2 , respectively. These results then also confirm the robustness of the proposed estimation approach despite the unavoidable presence of noise and discretization in the image acquisition. Note also that, as expected, the estimation error in the (active) case I reaches 'convergence' (i.e., drops below the threshold of 5 mm) significantly faster than case II (T 1 < T 2 ).

C. Structure Estimation for a Spherical Target

We now discuss some experimental results concerning the estimation of the radius of a spherical target: indeed, as explained in Sect. III-B, estimation of R allows to fully recover the sphere 3D position P 0 = sR where vector s is directly obtainable from image measurements, see [START_REF] Corke | Spherical Image-Based Visual Servo and Structure Estimation[END_REF]. As object to be tracked, we made use of a white table tennis ball placed on a black table and with a radius of 1.9 cm. As explained in Sect. III-B, the convergence rate of the estimation error for the sphere case only depends on the norm of the linear velocity v and not on its direction. This fact is proven by the first experiment where the estimation task is run twice starting from two different positions and imposing two different camera velocities but with same norm. These values were used during the experiments: [START_REF] Omari | Metric Visual-Inertial Navigation System Using Single Optical Flow Feature[END_REF], and v = (-0.05, 0, 0) m/s for case I and v = (0, 0.045, 0.02) m/s for case II, with v = 0.05 m/s in both cases. The camera angular velocity ω was exploited to keep (s x , s y ) ≃ (0, 0) (centered sphere).

αβ = 2 • 10 3 for gains Q and Λ, c 1 = c * 1 = 2 √ αβσ 1 for D 1 in
Figure 7a shows the behavior of the estimation errors (solid blue and red lines): note how the error transient response for the two cases is essentially coincident, and also equivalent to that of the reference second order system (9) with the desired poles, i.e., by setting σ 2 1 = v 2 = const and c 1 = c * 1 in (9) (dashed black line). The higher noise level in case II (red line) is due to the larger distance between the camera and the spherical target (see Fig. 7b) which negatively affects the estimation task. The standard deviation of the radius estimation error, computed on a time window of 1 s after z(t) has become smaller than 1 mm (vertical dashed lines in the plot), is 0.3 mm for case I and 0.2 mm for case II: we can note, again, the very satisfactory results obtained with the proposed estimation scheme in terms of accuracy of the reconstructed sphere radius. Note also how, in the two cases, the estimation error z(t) drops below the threshold of 1 mm at essentially the same time, as expected (same error transient respose).

Since the direction of the velocity does not play any role in this case, no optimization of σ 2 1 can be performed under the constraint v = const. On the other hand, the analysis of Sect. II-B clearly indicates the importance of choosing a proper value of c 1 for the damping matrix D 1 in (6). To show this fact, we report here three experiments characterized by the same camera trajectory of the previous case I, but by employing three different values for c 1 , that is, c * 1 , 2c * 1 and 0.5c * 1 . These correspond to a critically damped, overdamped and underdamped response for the ideal system [START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF], respectively. The experimental results reported in Fig. 8 show that the behavior of the estimation error z (solid lines) has an excellent match with that of (9) (represented by dashed lines), thus fully confirming (i) the validity of the proposed theoretical analysis,

D. Structure Estimation for a Cylindrical Target

In this final section we report some experimental results concerning the active estimation of the radius of a cylindrical object. Indeed, as in the sphere case, knowledge of R allows to fully recover the 3D point P 0 = Rs, with vector s from (33) and the cylinder axis a in (35) being directly obtainable from image measurements. For these experiments we used a white cardboard cylinder placed on a black table. The radius of the cylinder was approximately 4.2 cm.

In the cylinder case, the convergence rate of the estimation error depends both on the norm of the camera linear velocity v and on its direction w.r.t. the cylinder axis a, see (38). It is then interesting to optimize the direction of v under the constraint v = const for maximizing the eigenvalue σ 2 1 (i.e., so as to obtain the fastest convergence rate for a given 'control effort' v ). From (39), maximization of σ 2 1 (t) w.r.t. vector v can be obtained by choosing

v = J T v,1 -J † v,1 J a,1 [a] × ω, (41)
with A † being the pseudoinverse of a matrix A, i.e., by following the gradient of σ 2 1 w.r.t. v and by compensating for the (known) effects of input ω. In order to additionally enforce the constraint v = const during the eigenvalue maximization, eq. ( 41) can be modified as

v = v v 2 k 1 (κ des -κ) + k 2 I 3 - vv T v 2 (J T v,1 -J † v,1 J a,1 [a] × ω), (42)
with k 1 > 0 and k 2 > 0. Analogously to the point feature case, the first term in (42) asymptotically guarantees v(t) = v 0 while the second term projects (41) onto the null-space of the constraint v(t) = const. As for the angular velocity ω, we exploited it for keeping the axis of the cylinder at the center of the image plane by regulating (s x , s y ) to (0, 0).

We now present three experimental results structured as follows: in the first experiment (case I), the update rule (42) is fully implemented (k 1 > 0, k 2 > 0) for actively optimizing the direction of v. In the second experiment (case II), the camera starts from the same initial pose and velocity as in case I, but (42) is implemented with k 1 > 0 and k 2 = 0, i.e., without performing any optimization of σ 2 1 . Finally, in the third experiment (case III), the camera starts from a different initial pose and with a different velocity direction (but same norm) w.r.t. the previous two cases, and (42) is again fully implemented. This last case is meant to show how the convergence properties of the estimator are not affected by the direction of the camera linear velocity as long as it stays orthogonal to the cylinder axis a.

The experiments were run with the following conditions: αβ = 500 for gains Q and Λ, c 1 = c * 1 for D 1 in (6), k 1 = 10, k 2 = 1 for cases I and III, and k 2 = 0 for case II. As for the linear velocity, we set v(t 0 ) = v 0 = (-0.01, 0.05, 0.05) m/s for cases I and II, and v(t 0 ) = v 0 = (-0.05, 0.05, 0.01) m/s for case III (note how v 0 2 = 5.1 × 10 -3 m 2 /s 2 in all three cases).

The behavior of σ 2 1 (t) is shown in Fig. 9a: as explained at the end of Sect. III-C, under the constraint v = const, one has max v σ 2 1 = v 2 as the largest possible value for σ 2 1 (obtained when v T a = 0). It is then possible to verify that, indeed, σ 2 1 (t) → v 0 2 in cases I and III despite the different initial conditions of the experiments (different camera pose and direction of v). The optimization in (42) results in a null component of v along a, thus allowing to move faster in the 'useful' directions (while keeping a constant v ), and to increase the value of σ 2 1 to its maximum possible value. Also, note how the value of σ 2 1 (t) for case II results smaller than in the other two cases (as expected) since no optimization is present in this case.

The behavior of the estimation error z(t) is shown in Fig. 9b: again, we can note that the transient response for cases I and III results essentially coincident and in almost perfect agreement with that of the reference system (9) with desired poles (dashed black line). As expected, the response for case II (red line) is slower than in cases I and III. As in the point feature case, Fig. 9b reports, for the three cases under consideration, the times T 1 = 2.74 s, T 2 = 4.78 s and T 3 = 2.66 s at which the estimation error drops below the threshold 2 mm (vertical dashed lines). The standard deviation of the error evaluated on a time window of 1 s after convergence has been 'reached' resulted in the values of 0.4, 0.6 and 0.7 mm, respectively. We can then appreciate, again, the high accuracy of the proposed approach in estimating the cylinder radius R while also optimizing online for the camera motion. The higher estimation error in case III can be ascribed to the larger distance between the camera and the observed target, which increases the effect of discretization errors. Note also how T 1 ≈ T 3 < T 2 thanks to the active optimization of the error convergence rate. Finally, Fig. 9c depicts the camera trajectories for the three experiments with arrows indicating the direction of the optical axis. In case II the camera simply travels along a straight line (v(t) ≡ v 0 ), while in cases I and III the direction of v is suitably modified resulting in a trajectory lying on a plane orthogonal to a.

V. CONCLUSIONS

In this paper we have addressed the problem of active SfM for recovering the 3D structure of a point feature and of spherical and cylindrical objects by exploiting a novel active estimation strategy tailored to the three cases under consideration. For the depth estimation of a point feature, two possibilities differing in the adopted projection model (planar or spherical) were proposed and critically compared. The results showed the complementarity of the two models in terms of attainable convergence rates and basin of attraction for the estimation error. In the spherical and cylindrical cases, we instead showed how an adequate choice of the measured visual features allows to reduce the SfM task to the estimation of a single unknown constant quantity (the sphere/cylinder radius R) in place of the classical (and time-varying) three parameters (scaled normal vector of the planar limb surface). Availability of this quantity allows to then retrieve the full 3D structure of the observed targets. The reported experimental results fully confirmed the validity of the theoretical analysis and, in particular, the ability of the proposed active estimation strategy to impose, in all three cases, a desired transient response to the estimation error equivalent to that of a reference linear second-order system with desired poles.

We are currently investigating the use of similar active strategies for dealing with more complex 3D scenes. A possibility in this sense could be to decompose the SfM problem in two phases by (i) extracting and classifying, possibly from an initial measurement in the form of a point cloud, a set of primitive shapes belonging to the classes described in our work (points, spheres, cylinders or other 3D geometries), or also to other classes such as 2D planar patches made of discrete/dense sets of points; (ii) performing an (active) estimation of the whole scene structure by applying the same strategy presented here and by either sequencing the estimation of single scene components, or considering an "extended" system obtained by concatenating the observable and unobservable components corresponding to each of the basic shapes/classes.

We are also investigating how to extend our solution to problems involving the estimation of more parameters than the number of available measurements (i.e., with m < p), thus requiring to fulfil the more general observability condition [START_REF] Civera | Inverse Depth Parametrization for Monocular SLAM[END_REF]. Finally, we are also interested in the use of the proposed active strategy in the context of vision-based manipulation tasks. Some preliminary results in this context are reported in [START_REF] Spica | Coupling Visual Servoing with Active Structure from Motion[END_REF] where it is shown, and experimentally proven, that an online optimization of the estimation convergergence rate can improve the performance in executing visual servoing tasks.

APPENDIX A TIME-DERIVATIVE OF THE LIMB SURFACE PARAMETERS

FOR A SPHERICAL TARGET Differentiation of χ from (21) w.r.t. time yields

χ = Ṗ 0 K -P 0 K K 2 = Ṗ 0 K -2P 0 P T 0 Ṗ 0 K 2 (43)
which, being Ṗ 0 = -v -[ω] × P 0 and exploiting the property P T 0 [ω] × P 0 = 0, can be rewritten as

χ = - v K - [ω] × P 0 K + 2 P 0 P T 0 v K 2 = - v K -[ω] × χ + 2χχ T v.
(44) Letting s z = Z 0 /R (s z > 1), one then has

χ T χ - 1 s 2 z χ 2 z = X 2 0 + Y 2 0 + Z 2 0 K 2 - R 2 Z 2 0 Z 2 0 K 2 = 1 K . ( 45
)
This then shows how 1/K can be expressed in terms of χ and of s 2 z , with s z being directly obtainable from image measurements, see [START_REF] Corke | Spherical Image-Based Visual Servo and Structure Estimation[END_REF].

APPENDIX B ESTIMATION OF THE LIMB SURFACE PARAMETER FOR A CYLINDRICAL TARGET

In order to estimate the parameters of the limb surface associated to a cylindrical object, one could consider as measurement the 2 + 2 angle-distance parameters (θ i , ρ i ) of the straight lines resulting from the projection of the cylinder on the image plane. From [START_REF] Chaumette | Structure from controlled motion[END_REF], [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF], the interaction matrix in this case is given by:

L =     λ ρ1 c 1 λ ρ1 s 1 -λ ρ1 ρ 1 (1 + ρ 2 1 )s 1 -(1 + ρ 2 1 )c 1 0 λ θ1 c 1 λ θ1 s 1 -λ θ1 ρ 1 -ρ 1 c 1 -ρ 1 s 1 -1 λ ρ2 c 2 λ ρ2 s 2 -λ ρ2 ρ 2 (1 + ρ 2 2 )s 2 -(1 + ρ 2 2 )c 2 0 λ θ2 c 2 λ θ2 s 2 -λ θ2 ρ 2 -ρ 2 c 2 -ρ 2 s 2 -1     (46
) with s i = sin θ i , c i = cos θ i , and

λ ρi = -(χ x ρ i c i + χ y ρ i s i + χ z ) λ θi = χ y c i -χ x s i . (47)
Therefore, being (46) and (47) linear in the unknown χ, one can again apply the estimation scheme (2) with s taken as the vector of measured quantities on the image plane, i.e., s = (ρ 1 , θ 1 , ρ 2 , θ 2 ).

As for the dynamics of χ, since (21) still holds for a cylindrical object (see [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]), one can again exploit (43) with, however, in this case

Ṗ 0 = -I 3 -aa T v -[ω] × P 0 and thus χ = - 1 K I 3 -2χχ T I 3 -aa T v -[ω] × χ.
Finally, one can invoke (45) in order to express 1/K as a function of χ and s 2 z , with s z being the third element of vector s in (33).

APPENDIX C DERIVATION OF EQUATION (36)

We note that the cylinder axis a can be determined by the intersection of two planes P i : r T i X -d i = 0, i = 1, 2, with r 1 = [a] × P 0 P 0 , d 1 = 0, r 2 = -P 0 P 0 , d 2 = P 0 , (48) see Fig. 3. In particular, plane P 1 passes through the camera optical center, it is orthogonal to plane P 2 , and both planes contain the axis a passing through P 0 (by construction). Since Rs = P 0 and P 0 belongs to the cylinder axis a, we have Rr T i s -d i = 0, i = 1, 2 (the point Rs belongs to both planes P i ). Taking the time derivative of these latter constraints (with R = const), one has

r T i ṡ = 1 R ḋi -s T ṙi , i = 1, 2. (49)
Since ṙi = [r i ] × ω and ḋi = r T i v (see [START_REF] Giordano | 3D structure identification from image moments[END_REF]), eq. ( 49) can be rewritten as

r T i ṡ = 1 R r T i v -s T [r i ] × ω, i = 1, 2. (50)
Finally, from a T P 0 = 0 and P 0 = Rs we have a T s = 0 implying that

a T ṡ = -s T ȧ = -s T [a] × ω. (51)
We now note that equations ( 50) and (51) provide three linear constraints for ṡ which, by using (48), can be rearranged in matrix form as the following linear system It is easy to verify that the 3 × 3 matrix on the left hand side of (52) is orthonormal: by then solving (52) for ṡ and performing some simplifications we finally obtain the sought result ṡ = -

1 R I 3 -aa T [s] × u. (53)

Fig. 1 :

 1 Fig. 1: Optimality conditions for the camera linear velocity v as dictated by system (17). (a): when p = e 3 , vector v must be orthogonal to p and lie on the plane S spanned by p and e 3 (that is, v must belong to a specific straight line). (b): when p = e 3 , any v ⊥ e 3 is a valid solution to (17).

Fig. 2 :

 2 Fig. 2: Spherical target O s and planar limb surface L.

Fig. 3 :

 3 Fig. 3: Camera C and cylindrical target O c with the planar limb surface L and the other planes of interest P 1 and P 2

Fig. 4 :

 4 Fig. 4: Camera snapshots for the point feature (a), the sphere (b) and the cylinder (c) experiment.

Fig. 6 :

 6 Fig. 6: Experimental results for the point feature case. (a): behavior of the estimation error for case I (solid red line) and case II (solid blue line), and for an 'ideal' second order system (9) with desired poles at σ 2 max (dashed black line). The two vertical dashed lines indicate the times T 1 = 4.95 s and T 2 = 9.85 s at which the estimation error drops below the threshold of 5 mm. (b): Camera trajectories for case I (red line) and case II (blue line) with arrows indicating the direction of the camera optical axis. (c): behavior of σ 2 1 (t) for case I (red line) and case II (blue line).

Fig. 7 :

 7 Fig.7: Experimental results for the estimation of the radius of a sphere using different constant camera velocities with the same norm. (a): behavior of the estimation error z(t) for the two cases (solid blue and red lines), and for an 'ideal' second order system with poles at the desired locations (dashed black line). The vertical dashed lines indicate the times at which the estimation error z(t) drops below the threshold of 1 mm. (b): camera trajectories for case I (blue line) and case II (red line) with arrows indicating the direction of the camera optical axis.

1 c1 = c * 1 c1 = 2c * 1 Fig. 8 :

 1118 Fig.8: Experimental results for the estimation of the radius of a sphere with c 1 = c * 1 (blue line), c 1 = 2c * 1 (green line) and c 1 = 0.5c * 1 (red line). The dashed lines represent the response of an 'ideal' second order system with the corresponding poles. Note again the almost perfect match between the plots.

Fig. 9 :

 9 Fig. 9: Experimental results for the estimation of the radius of a cylinder with the following color coding: blue -case I, red -case II, green -case III. (a): behavior of σ 2 1 (t) for the three cases (coincident for cases I and III and larger than in case II). (b): behavior of z(t). The three vertical dashed lines indicate the times T 1 = 2.74 s , T 2 = 4.78 s and T 3 = 2.66 s at which the estimation error drops below the threshold of 2 mm. Note how T 1 ≈ T 3 and T 1 < T 2 as expected. (c): two views of the camera trajectories for the three cases with arrows indicating the direction of the camera optical axis.

This is due to the well-known fact that, under perspective and spherical projection, the motion in the image induced by pure rotations of the camera (i.e., when v = 0) does not depend on the structure of the scene.

The stability proof requires some additional technical assumptions on the regularity of the vanishing disturbance g (locally Lipschitz in a neighbourhood of the origin), on its growth bound w.r.t. e (which, since g(•) → 0 if (v, ω) → 0, can always be made small enough by limiting (v, ω)), and on the norm of the initial error e(t 0 ) . The interested reader can find in[START_REF] De Luca | Feature depth observation for image-based visual servoing: Theory and experiments[END_REF] a detailed derivation of the proof.

The ground truth Z 0 (t) was obtained from a previous offline estimation of the 3D position P 0 in the world frame, and by then using the information on the camera position provided by the robot forward kinematics.

The value of v resulting from this (and following) optimizations was numerically integrated so as to obtain the commanded v(t) sent to the robot low-level controller.