Jihong Zhu

David Navarro-Alarcon
email: david.navarro-alarcon@polyu.edu.hk

Robin Passama

Andrea Cherubini

Vision-based Manipulation of Deformable and Rigid Objects Using Subspace Projections of 2D Contours

Keywords: Visual servoing, sensor-based control, deformable object manipulation

This paper proposes a unified vision-based manipulation framework using image contours of deformable/rigid objects. Instead of explicitly defining the features by geometries or functions, the robot automatically learns the visual features from processed vision data. Our method simultaneously generatesfrom the same data-both visual features and the interaction matrix that relates them to the robot control inputs. Extraction of the feature vector and control commands is done online and adaptively, and requires little data for initialization. Our method allows the robot to manipulate an object without knowing whether it is rigid or deformable. To validate our approach, we conduct numerical simulations and experiments with both deformable and rigid objects.

Introduction

Humans are capable of manipulating both rigid and deformable objects. However, robotic researchers tend to consider the manipulation of these two classes of objects as separate problems. Unless otherwise mentioned, object rigidity is an implicit assumption in most manipulation tasks. On the other hand, methods designed for deformable object manipulation [START_REF] Sanchez | Robotic manipulation and sensing of deformable objects in domestic and industrial applications: a survey[END_REF], are never applied on rigid objects. This paper presents our efforts in formulating a generalized framework for vision-based manipulation of both rigid and deformable objects, which does not require prior knowledge of the object's mechanical properties.

In the visual servoing literature [START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF], vector s denotes the set of features selected to represent the object in the image.

These features represent both the object's pose and its shape. We denote the process of selecting s as parameterization. The aim of visual servoing is to minimize, through robot motion, the feedback error e = s * -s between the target s * and the current (i.e., measured) feature s.

One of the initial works on vision-based manipulation of deformable objects is presented in [START_REF] Inoue | Hand-eye coordination in rope handling[END_REF] to solve a knotting problem by a topological model. Smith et al. developed a relative elasticity model, such that vision can be utilized without a physical model for the manipulation task [START_REF] Smith | Vision based manipulation of non-rigid objects[END_REF]. A classical model-free approach in manipulating deformable objects is developed in [START_REF] Berenson | Manipulation of deformable objects without modeling and simulating deformation[END_REF]. More recent research (Lagneau et al., 2020a) and (Lagneau et al., 2020b) proposes a method for online estimation of the deformation Jacobian, based on weighted least square minimization with a sliding window. In [START_REF] Navarro-Alarcon | On the visual deformation servoing of compliant objects: Uncalibrated control methods and experiments[END_REF] and [START_REF] Navarro-Alarcon | Fourier-based shape servoing: A new feedback method to actively deform soft objects into desired 2d image contours[END_REF], the vision-based deformable objects manipulation is termed as shape servoing. An expository paper on the topic is available in [START_REF] Navarro-Alarcon | On model adaptation for sensorimotor control of robots[END_REF]. A recent work on vision-based shape servoing of plastic material was presented in [START_REF] Cherubini | Model-free vision-based shaping of deformable plastic materials[END_REF].

For a detailed survey on shape servoing we refer readers to [START_REF] Sanchez | Robotic manipulation and sensing of deformable objects in domestic and industrial applications: a survey[END_REF]. For shape servoing, commonly selected features are curvatures [START_REF] Navarro-Alarcon | On the visual deformation servoing of compliant objects: Uncalibrated control methods and experiments[END_REF], points (Wang et al.) and angles (Navarro-Alarcon and Liu, 2013). Laranjeira et al. proposed a catenary-based feature for tethered management on wheeled and underwater robots [START_REF] Laranjeira | Catenary-based visual servoing for tethered robots[END_REF][START_REF] Laranjeira | Catenary-based visual servoing for tether shape control between underwater vehicles[END_REF]. A more general feature vector is that containing the Fourier coefficients of the object contour (Navarro-Alarcon and Liu, 2018; [START_REF] Zhu | Dualarm robotic manipulation of flexible cables[END_REF]. Yet, all these approaches require the user to specify a model, e.g., the object geometry [START_REF] Wang | A unified controller for region-reaching and deforming of soft objects[END_REF]Navarro-Alarcon and Liu, 2013;[START_REF] Navarro-Alarcon | On the visual deformation servoing of compliant objects: Uncalibrated control methods and experiments[END_REF] or a function [START_REF] Laranjeira | Catenary-based visual servoing for tethered robots[END_REF][START_REF] Navarro-Alarcon | Fourier-based shape servoing: A new feedback method to actively deform soft objects into desired 2d image contours[END_REF][START_REF] Zhu | Dualarm robotic manipulation of flexible cables[END_REF] for selecting the feature. Alternative data-driven (hence, model-free) approaches rely on machine learning. Nair et al. combine learning and visual feedback to manipulate ropes in [START_REF] Nair | Combining self-supervised learning and imitation for vision-based rope manipulation[END_REF]. Li et al. approximate the deformation and camera model using a neural network [START_REF] Li | Vision-based robotic grasping and manipulation of usb wires[END_REF]. The authors of [START_REF] Hu | 3-d deformable object manipulation using deep neural networks[END_REF] em-ploy deep neural networks to manipulate deformable objects given their 3D point cloud. All these methods rely on (deep) connectionist models, which invariably require training through an extensive data set. The collected data has to be diverse enough to generalize the model learnt by this type of networks. Instead of relying on algorithmic solutions, [START_REF] She | Cable manipulation with a tactile-reactive gripper[END_REF] utilizes a vision-based tactile sensor (GelSight) for manipulating cables.

It is noteworthy that some of the above mentioned methods may apply to rigid objects. Yet, none of the previous works has investigated the possibility of this extension nor reported its experimental validation, as we do in this paper.

The trend in visual servoing, when controlling the pose of rigid objects is to find features which are independent from the object characteristics. Following this trend, [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF] proposes the use of image moments. More recently, researchers have proposed direct visual servoing (DVS) methods, which eliminate the need for user-defined features and for the related image processing procedures. The pioneer DVS works [START_REF] Collewet | Visual servoing set free from image processing[END_REF][START_REF] Collewet | Photometric visual servoing[END_REF] propose using the whole image luminance to control the robot, leading to "photometric" visual servoing. Bakthavatchalam et al. join the two ideas by introducing photometric moments [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF]. A subspace method [START_REF] Marchand | Subspace-based direct visual servoing[END_REF] can further enhance the convergence of photometric visual servoing, via Principal Component Analysis (PCA). This method was first introduced for visual servoing in [START_REF] Nayar | Subspace methods for robot vision[END_REF]. In that work, using an eye-in-hand setup, the image was compressed to obtain a low-dimensional vector for controlling the robot to a target pose. Similarly, the authors of [START_REF] Deguchi | Visual servoing using eigenspace method and dynamic calculation of interaction matrices[END_REF] transformed the image into a lower dimensional hyper surface, to control the robot position via in-hand camera feedback. However, DVS generally considers rigid and static scenes, where the robot controls the motion of the camera (eye-in-hand setup) to change only the image viewpoint, and not the environment. These constraints on the setup avoid breaking the Lambertian hypothesis that is needed, since DVS relies on the raw image luminance, which should not vary with the viewpoint. For this reason, to our knowledge, DVS was never applied to object manipulation, since changes in the pose and/or shape of the object would break the Lambertian assumption. This is not the case of feature-based methods (such as the one we present here), as long as the feature is chosen to be reliable even when the viewpoint and/or scene change.

Compared with the above-mentioned works, our paper presents the fol-lowing original contributions:

1. We propose to use a feature vector -based on PCA of sampled 2D contours -for model-free manipulation of both deformable and rigid objects.

2. We exploit the linear properties of PCA and of the local interaction matrix, to initialize our algorithm with little data -the same data for feature vector extraction and for interaction matrix estimation. 3. We report experiments using the same framework to manipulate objects with different unknown geometric and mechanical properties.

The paper is organized as follows. Sect. 2 presents the problem. Sect.

3 outlines the framework. Sect. 4 elaborates on the methods. In Sect. 5, we analyze and verify the methods by numerical simulations. Then, Sect. 6 presents the robotic experiments and we conclude in Sect. 7.

Problem statement

In this work, we aim at solving object manipulation tasks with visual feedback. We rely on the following hypotheses:

• The shape and pose of the object are represented by its 2-D contour on the image as seen from a camera fixed in the robot workspace (eyeto-hand setup). We denote this contour as

c = [p 1 • • • p K ] T ∈ R 2K , (1)
where p j = [u j v j ] ∈ I denotes the jth pixel of the contour in the image I.

• The contour is always entirely visible in the scene and there are no occlusions.

• One of the robot's end-effectors holds one point of the object (we consider the grasping problem to be already solved). At each control iteration i, its pose is r i ∈ SE (3), and it can execute motion commands δr i ∈ SE (3) that drive the robot so that r i+1 = r i + δr i .

• The target constant shape (i.e., contour) of the object, c * , is physically reachable with shaping motions of the grasping point r. To ensure this hypothesis, one can first command the robot to verify that it can move the shape to c * . Problem Statement. Given a target shape of the object, represented by a constant contour vector c * , we aim at designing a vision-based controller that generates a sequence of robot motions δr i to drive the initial contour to the target one.

The controller should work without any knowledge of the object physical characteristics, i.e., for both rigid and deformable objects. In the latter case, we assume that the deformation is homogeneous. Since rigid and deformable objects behave differently during manipulation, we set the following manipulation goals:

• Rigid objects: move them to a target pose (see Fig. 1a).

• Deformable objects: move them to a target pose with a target shape (see Fig. 1b).

The formulation of the problem is general, but due to challenges in perception (discussed in Sect. 7), we carried out the cases of study with movements in SE (2).

Preliminary

In this section, we present an overview of the proposed approach, motivated by the problem analysis. Throughout the paper, we use c to indicate the object contour and s as the feature vector obtained from the contour.

The subscript i indicates the instance of the variable at iteration i (e.g., c i is the contour at iteration i).

We can work directly on the object shape space by selecting the contour as the feature vector s ≡ c ∈ R 2K . With image and data processing, we can extract a fixed number of ordered (i.e., identified) contour points to represent the shape/pose of the object. However, this will result in an unnecessarily large dimension of the feature vector (e.g., if K = 50, s has 100 components). The high dimensional feature vector increases the computation demand and complicates the control due to the high under-actuation of the system. Therefore, instead of working on this feature vector, we work on one with smaller dimensions. To this end, we split the problem into two sub-problems: parameterization and control, see Fig. 2. Parameterization consists in representing the contour via a compact feature vector s ∈ R k , such that k << 2K. We denote this representation as s = g(c). We introduce the method for parameterization in Sect. 4.1.

Control consists in computing robot motions δr 1 , δr 2 , . . . , so that the object's representation s converges to the target s * . Control can be broken down to solving the optimization problem:

r * = arg min r (f (r) -s * ) (2)
where s = f (r) denotes the mapping between robot pose and feature vector, which is assumed to be smooth and generally nonlinear. The smoothness assumption requires that the objects' contour is at least twice differentiable with respect to the robot motion. If we know the analytic solution to f (r), we can solve (2) and obtain the target shape in a single iteration by commanding r * . A solution to this problem is to approximate the full mapping f (r) from sensor observations. Classic deep learning-based approaches typically require a long training phase to collect vast and diverse data for approximating f (r).

In some cases (for instance, robotics surgery), it is not possible to collect such data beforehand. Moreover, if the object changes, new data has to be collected to retrain the model, leading to a cumbersome process. In this paper instead, we aim at doing the data collection online, with minimum initialization. Thus, instead of estimating the full nonlinear mapping f (r), we divide it into piece-wise linear models [START_REF] Sang | Adaptive control of piecewise linear systems: the state tracking case[END_REF] at successive equilibrium points. The locality assumption refers to both the time and spatial dimensions. These models are considered time invariant in the neighbourhood of the equilibrium points. We then compute the control law for each linear model and apply it to the robot end-effector. We will dedicate Sections 4.2 and 4.3 to the local models and Sections 4.4 and 4.5 to derive the control inputs and to analyze (local) stability.

Methodology

Given a target shape c * , we define an intermediate local target c * i at each i = 1, 2, . . . (see Fig. 2). At the i th iteration, the robot autonomously generates a local mapping g i to produce the feature vector s i = g i (c i ). The robot then finds the local mapping s i = f i (r i ) online.

Consider at the current time instant i, the shape c i , the intermediate target c * i and the local parameterization g i . We can transform shape data into a feature vector by:

s i = g i (c i ), s * i = g i (c * i ). (3)
The linearized version of s = f (r) centered at (s i , r i ) is then:

δs i = L i δr i , (4)
with

L i = ∂f i ∂r | r=r i , δs i = s i+1 -s i , δr i = r i+1 -r i .
(5)

The matrix L i represents a local mapping, referred to as the interaction matrix in the visual servoing literature [START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF].

If L i can be estimated online at each iteration i, then, we can design one-step control laws to drive s i towards s * i .

After the robot has executed the motion command δr i , we update the next target to be s * i+1 , and so on, until it reaches the final target s * . Although the validity region of this local mapping is smaller than that of the original nonlinear mapping, it enables to use an online training approach that requires less data and reduced computational demand.

Figure 3 shows the building blocks of the overall framework. In this section, we focus on the red dashed part of the diagram. We will elaborate on each red block in the subsequent subsections. The blue block represents the image processing pipeline that will be discussed in Sect. 6.1.

Feature vector extraction

There are many ways to parameterize c in order to reduce its dimension. One of the prominent dimension reduction methods is Principal Component Analysis (PCA). PCA finds a new orthogonal basis for high-dimensional data. This enables projection of the data to lower dimension with the minimal sum of squared residuals. It was used in image processing [START_REF] Zhang | Two-stage image denoising by principal component analysis with local pixel grouping[END_REF] and classification [START_REF] Zeng | Color image classification via quaternion principal component analysis network[END_REF]. In visual servoing, the method was first introduced in [START_REF] Nayar | Subspace methods for robot vision[END_REF]. PCA is proven to be an effective, yet easy to implement, algorithm for dimension reduction. By projecting to the new orthogonal space, each feature component is linearly independent. Besides, by checking the explained variance of a feature, we can intuitively measure if it represents the original shape.

We apply PCA to reduce c ∈ R 2K to s ∈ R k . To find the projection, we collect M images with different shapes of the object and construct the data

matrix Γ = [c 1 c 2 • • • c M ] ∈ R 2K×M . Then, we shift the columns of Γ by the mean contour c = M i=1 c i /M : Γ = [c 1 -c c 2 -c • • • c M -c] ∈ R 2K×M . ( 6
)
We then compute the covariance matrix C = ΓΓ T , and apply Singular Value Decomposition (SVD) to it:

C = U ΣV T . ( 7
)
Once we have obtained the eigenvector matrix U ∈ R 2K×2K , we can move on to select the first k columns1 of U denoted by U (k) ∈ R 2K×k . Then, the 2K-dimensional contour c can be projected into a smaller k-dimensional feature vector s as:

s = U T (k)(c -c) ∈ R k . (8)
To assess the quality of this projection, we can compute the explained variance using the eigenvalue matrix Σ ∈ R 2K×2K in (7). By denoting the diagonal entries of Σ as σ 1 , • • • , σ 2K , the explained variance of the first k components is:

Υ(k) = k j=1 σ j 2K j=1 σ j . ( 9
)
where Υ is a scalar between 0 and 1 (since σ j > 0, ∀j), indicating to what extent the k components represent the original data (a larger Υ suggests a better representation).

Since PCA calculate features that lie on an orthonormal basis, these features are linearly independent. For controlling n DoF, at least the same number of independent visual features should be used. Therefore, we set k = n features.

Local target generation

Let us now explain how we generate a local target contour c * i given a current contour c i and final target contour c * . We also show this in Algorithm 1. The overall shape error is given by:

c e = c * -c i . ( 10
)
We define the intermediate target contour as:

c * i = c i + 1 η c e , (11)
with η = 1, 2, . . . an integer that ensures that c * i is a "good" local target for c i (i.e., the two are similar). Therefore, if we project the intermediate local data using the eigenvector matrix at the current iteration,

U i ∈ R 2K×2K
(note that we are using the full projection matrix and not just the first k columns), the projection

s p i = U i (c * i -c) ∈ R 2K should fulfil: Ψ (k) = k j=1 s p i,j 2K j=1 s p i,j ≥ , (12)
with ∈ [0; 1] a threshold and s p i,j the j-th component of the projection. Then, we select the first k components in s p i to be the local target

s * i ∈ R k .
Algorithm 1 outlines the steps for computing the local intermediate targets, so that:

• they are near the final target,

• the corresponding feature vector can be extracted with the current learned projection matrix.

Remark 1. The reachability of a local target can only be verified with a global deformation model which we want to avoid identifying in our methods. We will further discuss this issue in the Conclusion (Sect. 7).

Interaction matrix estimation

Let us consider the current contour c i and the local target c * i . In this section, we show how we can implement the PCA and model estimation

Algorithm 1 Local target generation localTargetFound = false Ψ 0 = 0 η = 1 while not localTargetFound do c * i = c i + 1 η (c * -c i ) s p i = U i c * i Ψ η = k j=1 s p i,j / 2K j=1 s p i,j if Ψ η ≥ or Ψ η < Ψ η-1 then localTargetFound = true s * i = [I 0]s p i end if η = η + 1 end while
together and online. We denote the robot motions and corresponding object contours over the last M iterations (prior to iteration i, with i ≥ M ) as:

∆R i = δr i-M +1 δr i-M +2 • • • δr i ∈ R n×M Γ i = c i-M c i-M +1 c i-M +2 • • • c i ∈ R 2K×(M +1) , (13)
with M the number of data samples collected during initialization, i.e., the size of the sliding window used for model adaptation (see Sect. 4.5).

By selecting k = n (note that n is also the number of DoFs of the robot manipulator we considered in the task execution), we compute the projection matrix U i (n) ∈ R 2K×n , from Γ i and ci via ( 6) and (7). Then, using U i (n), we project current contour c i , target contour c * i and shape matrix Γ i :

s i = U i (n) T (c i -ci ) ∈ R n , s * i = U i (n) T (c * i -ci ) ∈ R n , S i = U i (n) T Γi = s i-M s i-M +1 • • • s i ∈ R n×(M +1) . (14)
In ( 14), Γi is normalized by ci as in (6). We can then compute ∆S i from ( 5) and ( 14), by subtracting consecutive columns of S i :

∆S i = δs i-M +1 δs i-M +2 • • • δs i ∈ R n×M . ( 15
)
Using ∆S i ∈ R n×M and ∆R i ∈ R n×M we can now estimate the local interaction matrix L i ∈ R n×n at iteration i. We assume that near this iteration, the system remains linear and time invariant: L i is constant. Using the local linear model (4), we can write the following:

∆S i = L i ∆R i . ( 16
)
Our goal then is to solve for L i , given ∆S i and ∆R i . Note that this is an overdetermined linear system (with n×M equations for n 2 unknowns). Let us consider ∆R i ∈ R n×M has full row rank. Note this sufficiently implies M ≥ n. With this prerequisite, rank

(∆R i ) = n. Therefore, rank(∆R i ∆R T i ) = n,
and its inverse exists. We post multiply ( 16) by ∆R T i :

∆S i R T i = L i ∆R i ∆R T i . (17)
Then, since ∆R i ∆R T i is invertible, the L i that best fulfills ( 16) is:

Li = ∆S i ∆R T i (∆R i ∆R T i ) -1 . (18)
If, in practice, the full row rank condition of ∆R i is not satisfied, rank(∆R i ∆R T i ) < n and ∆R i ∆R T i becomes singular. Then, instead of ( 18), we can use Tikhonov regularization:

Li = ∆S i ∆R T i (∆R i ∆R T i + λI) -1 , (19)
with λ an arbitrary (generally small) scalar.

Practically, this implies that one or more inputs motions do not appear in ∆R i . Therefore, we cannot infer the relationship between these motions and the resulting feature vector changes. In this case it is better to increase M and obtain more data, so that ∆R i has full row rank.

Instead of computing the interaction matrix, it is also possible to directly compute its inverse, since this guarantees better control properties [START_REF] Lapresté | An efficient method to compute the inverse jacobian matrix in visual servoing[END_REF]. With the same data, one can re-write (16) as:

L ⊕ i ∆S i = ∆R i . ( 20
)
We can also solve (20) with Tikhonov regularization:

L⊕ i = ∆R i ∆S T i (∆S i ∆S T i + λI) -1 . (21)

Control law and stability analysis

We can now control the robot, with either of the following strategies:

δr i = -α L † i (s i -s * i ), (22)
if one estimates the interaction matrix with ( 19), where † denotes the pseudoinverse, or:

δr i = -α L⊕ i (s i -s * i ) ( 23
)
if one estimates the inverse of the interaction matrix with (21). In both equations, α > 0 is an arbitrary control gain.

Proposition 1. Consider that locally, the model (4) closely approximates the interaction matrix L i = Li . For M number of linearly independent displacement vectors δr such that the interaction matrix Li is invertible, the update rule ( 22) asymptotically minimizes the error e i = s * i -s i , where s * i is the local target.

Proof. With δs i = s i+1 -s i , we can write (4) in discretized form as

s i+1 = s i + L i δr i . (24)
From the definition of e i we have (Note here the target s * i is not updated with i since we want to prove local convergence to a constant target):

e i = s * i -s i e i+1 = s * i -s i+1 (25)
Taking ( 24) into ( 25):

e i+1 = s * i -s i+1 = s * i -s i -L i δr i = e i -L i δr i . ( 26
)
We replace δr i in (26) with the control ( 22), the error dynamic is then:

e i+1 = e i -αL i L -1 i (s * i -s i ) = e i -αe i = (1 -α)e i . ( 27
)
is asymptotically stable for α ∈ [0; 1]. This can be proved by considering the Lyapunov function

V(e) = e T e. (28)
Using the error dynamic ( 27), one can derive:

∆V = V(e i+1 ) -V(e i ) = e T i ((1 -α) 2 -1)e i < 0. ( 29
)
This proves the local asymptotic stability of the error e using our inputs.

Model adaptation

Since both the projection matrix U T (n) and the interaction matrix are local approximations of the full nonlinear mapping, they need to be updated constantly. We choose a receding window approach with window size M .

At current iteration i, we estimate the projection matrix U T i and local interaction matrix L i with M samples of the most recent data. Using the interaction matrix and the local target c * i , we can derive the one-step command δr i by ( 22). Once we execute the motion δr i , a new contour data c i+1 is obtained. We move to the next iteration i + 1. A new pair of input and shape data [δr i , c i+1 ] is obtained. We shift the window by deleting the oldest data in the window and add in the new data pair. Then, using the shifted window, we compute one step control at iteration i + 1. The receding window approach ensures that, at each iteration, we are using the latest data to estimate the interaction matrix. The overall algorithm is initialized with small random motions around the initial configuration.

First, M samples of shape data and the corresponding robot motions are collected. With this initialization, we can simultaneously solve for the projection matrix and estimate the initial interaction matrix using the methods described in Sect. 4.1 and 4.3. Using the projection matrix and the initial/target shapes, we can then find an intermediate target (see Sect. 4.2). We consider quasi-static deformation. Hence, at each iteration the system is in equilibrium and can be linearized according to (4). The data that best captures the current system are the most recent ones. The choice of M is a trade-off between locality and richness. For fast varying deformations2 , we would expect to reduce M since a larger M will hinder the locality as-

sumption. Yet, if M is too small, it affects the estimation of Li (refer to the detailed discussion in Sect. 4.3).

Simulation results

In this section, we present the numerical simulations that we ran to validate our method.

Simulating the objects

We ran simulations on MATLAB (R2018b) with two types of objects: a rigid box and a deformable cable, both constrained to move on a plane. The rigid object is represented by a uniformly sampled rectangular contour. The controllable inputs are its position and orientation. For the cable, we developed a simulator, which is publicly available at https://github.com/Jihong-Zhu/cableModelling2D. The simulator relies on the differential geometry cable model introduced in [START_REF] Wakamatsu | Static modeling of linear object deformation based on differential geometry[END_REF], with the shape defined by solving a constrained optimization problem. The underlying principle is that the object's potential energy is minimal for the object's static shape [START_REF] Wakamatsu | Modeling of linear objects considering bend, twist, and extensional deformations[END_REF]. Position and orientation constraints (imposed at the cable ends) are input to the simulator. The output is the sampled cable. Figures 456, 9, 10 show simulated shapes of cables and rigid boxes.

We choose K = 50 samples for both rigid objects and cables. The camera perspective projection is simulated, with optical axis perpendicular to the plane. In each sub-figure, the solid red lines are the initial shapes and the dashed black are the shapes resulting from 10 random motions of the right tip (translations limited to ±5% of the length, rotations limited to ±5 • ).

Selecting the feature dimension k

To check whether choosing k = n can represent the shape accurately, we simulate 6 trials with distinct initial shapes of a cable. The dimension of the robot motion vector δr is n = 3 (two translations and one rotation of the right tip), and the motions are limited: each translation to ±5% of the cable length and the rotation to ±5 • . This range of motion gives a rule of thumb, which we have used for generating random movements throughout our experiments. For each trial, we command M = 10 random motions around the initial shape using our simulator. Figure 4 shows the 6 initial cable shapes (solid red) and the resulting shapes from 10 random movements (dashed black).

For each trial, we apply PCA to map the cable contour c ∈ R 2K to feature vector s ∈ R k , as explained in Sect. 4.1. We do this for k = 1, 2 and 3 and for each of these 18 experiments, we calculate the explained variance Υ(k) with ( 9). Table 1 shows these explained variances. In all 6 trials, k = n = 3 yields explained variances very close to 1. This result confirms that choosing k = n as the dimension of the feature vector gives an excellent representation of the shape data. It is also possible to select k = 2, since the first two components can represent more than 99% of the variance. Nevertheless, the simulation is noise-free. Therefore, although Υ(k) increases little from k = 2 to k = 3, this increase is not related to noise but to an actual gain in data information.

Table 1: Explained variance Υ(k) for the 6 trials with small motion.

trial 1 trial 2 trial 3 trial 4 trial 5 trial 6 k = 1 0.727 0.795 0.871 0.847 0.847 0.705 k = 2 0.992 0.995 0.996 0.997 0.997 0.994 k = 3 0.999 0.999 0.999 0.999 0.999 0.999 At this stage, it is legitimate to ask: how does this scale to larger movements? Figure 5 illustrates 10 cable shapes generated by large movements (angle variation: [-π 2 , π 2 ], maximum translation: 106%). Again, we apply PCA (M = 10); Table 2 shows the Υ(k) resulting from various values of k. Comparing Tables 1 and2, it is noteworthy that Υ(4) with large motion is smaller than Υ(2) with small motion. There are two possible explanation here. One is that when shapes stays local, the local linear mapping L in (4) remains constant and we need less features to characterize it; the more the shape varies, the more features we need. Another possible explanation is that for larger motions, M = 10 shapes may be insufficient for PCA. Likely, the larger the changes, the larger the number of shapes M needed.

Manipulation of deformable objects

With our cable simulator, we can now test the controller to modify the shape from an initial to a target one. Again, the left tip of the cable is fixed, and we control the right tip with n = 3 degrees of freedom (two translations and one rotation). Using the methods described in Sect. 4, we choose window size M = 5, the Tikhonov factor λ = 0.01, the local target threshold = 0.8, the control gain α = 0.01. To quantify the effectiveness of our algorithms in driving the contour to c * , we define a scalar measure: the Average Sample Error (ASE). At iteration i, with current contour c i it is:

ASE = c i -c * 2 2K . (30)
A small ASE indicates that the current contour is near the target one. In Sect. 4.4, we have proved that our controller asymptotically stabilizes the feature vector, s to s * . Hence, since we have also shown that s is a "very good" representation" of c, we also expect our controller to drive c to c * , thus ASE to 0. This measure is also used in the real experiments.

Using the cable simulator, we compare the convergence of two control laws proposed in our paper ( 22) and ( 23) against a baseline algorithm in [START_REF] Zhu | Dualarm robotic manipulation of flexible cables[END_REF] which uses Fourier parameters as feature. To make methods compatible, we choose first order Fourier approximation. Note that this results in a feature vector of dimension of 6 (see [START_REF] Zhu | Dualarm robotic manipulation of flexible cables[END_REF]) which is still twice the number k used in our method. We also normalize the computed control action and then multiply by the same gain factor 0.01.

We also introduced artificial noise to the contour data, to test the robustness of our method. For a unit length cable, we add Gaussian noise of zero mean and 0.01 standard deviation to the contour sample points. Fig. 6b shows that our algorithm converges in these conditions as well. It is worth mentioning that in robotic experiments, as shape data is obtained and sampled from the images, signal noise is inevitable. Yet, our framework is robust enough to still converge to the target shape/pose (see Sect. 6 for detailed real robot experiments).

Figure 6 shows two other simulation results: on the left, a reachable target and on the right an unreachable one. In Fig. 6a, the cable shape successfully evolves towards the target thanks to our controller (23). Figure 6c shows the starting shape (blue), unreachable target (dash black) and final shape (solid black) obtained using ( 23). Note that the controller gets stuck in a local minimum.

Figure 7 compares the evolution of ASE with our methods against the Fourier-based method for the reachable target; in the same figure, we also plot the evolution of ASE using (23) for the unreachable target. We can observe that our method provides faster convergence using half the features than [START_REF] Zhu | Dualarm robotic manipulation of flexible cables[END_REF]. Also, directly computing the inverse (23) provides faster convergence than (22). It is noteworthy to point out that the Fourier-based method requires a different parameterization for closed and open contours (see [START_REF] Navarro-Alarcon | Fourier-based shape servoing: A new feedback method to actively deform soft objects into desired 2d image contours[END_REF] and [START_REF] Zhu | Dualarm robotic manipulation of flexible cables[END_REF]), whereas in our framework, the parameterization can be kept the same. Last but not least, our approach is the only one among the three, which has been validated on both rigid and deformable objects.

Comparison with the Broyden update law

The Broyden update law [START_REF] Broyden | A class of methods for solving nonlinear simultaneous equations[END_REF], has been used to update the interaction matrix in classic visual servoing [START_REF] Hosoda | Versatile visual servoing without knowledge of true jacobian[END_REF][START_REF] Jagersand | Experimental evaluation of uncalibrated visual servoing for precision manipulation[END_REF][START_REF] Chaumette | Visual servo control, part II: Advanced approaches[END_REF] and shape servoing (Navarro-Alarcon et al., 2013).

In this section, we compare it with our method for updating the interaction matrix (19), which relies on a receding horizon. We will hereby show why the Broyden update law is not applicable in our framework.

Reachable Target Unreachable Target

Figure 7: The evolution of the ASE of the simulated cable manipulation using our method against the Fourier-based method as baseline and the ASE of the unreachable target with (23).

The Broyden update is an iterative method for estimating L i at iteration i. Its standard discrete-time formulation is:

Li = Li-1 + β δs i-1 -Li δr i-1 δr T i-1 δr i-1 δr T i-1 , ∀r i-1 = 0 (31)
with β ∈ [0; 1] an adjustable gain. Using our simulator, we estimate the interaction matrix using both Broyden update (with three different values of β) and our receding horizon method ( 19). We then compare (with L estimated with either method) the one-step prediction of the resulting feature vector:

ŝi+1 = Li r i + s i , (32)
with the ground truth s i+1 from the simulator. The results (plotted in Fig. 8)

show that receding horizon outperforms all three Broyden trials. One possible reason is that the components of s fluctuate since (at each iteration) a new matrix U is used. These variations cause the Broyden method to accumulate the result from old interaction matrices, and therefore perform badly on a long term. This result contrasts with that of (Navarro-Alarcon et al., 2013), where the Broyden method performs well since there is a fixed mapping from contour data to feature vector. Another advantage of the receding horizon approach is that it does not require any gain tuning.

Manipulation of rigid objects

The same framework can also be applied to rigid object manipulation.

Consider the problem of moving a rigid object to a certain position and orientation via visual feedback. This time, the shape of the object does not change, but its pose will (it can translate and rotate). We use the same M , λ, and α as for cable manipulation. We compare the convergence of two control laws proposed in our paper ( 22) and ( 23) against a baseline using image moments [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]. The translation and orientation can be represented with image moments and the analytic interaction matrix can be computed as explained in [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]). To make methods compatible, we normalize the computed control and then multiply it by the same factor 0.01.

Figure 9 shows two simulations where our controller successfully moves a rigid object from an initial (blue) to a target (dashed black) pose using control law (23). Figure 11 compares convergence of our methods against the image moments method. We can observe that our method provides a The topmost, middle and bottom plots show the one step prediction of s 1 , s 2 and s 3 , respectively. In all plots, the dashed red curve is the ground truth from the simulator. The plots clearly show that the receding horizon approach outperforms all three Broyden trials.

slightly slower convergence. Directly computing the inverse (23) provides a convergence similar to (22). Later, we will show why our method is slower.

Yet, the fact that it can be applied on both deformable and rigid objects makes it stand out over the other techniques.

Feature analysis for rigid objects

In this section, we analyze locally what each component of the feature vector represents, in the case of rigid object manipulation. To this end, we apply M = 10 random movements (rotation range [-0.11, 0.09], maximum translation 15% of the width) to multiple rigid rectangular objects (see Fig. 10). We compute the projection matrix as explained in Sect. 4.1, and transform the contour samples to feature vectors. Then, we seek the relationshipat each iteration -between the object pose x, y, θ and the components of the feature vector s generated by PCA. To this end, we use bivariate correlation [START_REF] Feller | An introduction to probability theory and its applications[END_REF] defined by: In Table 3, we show one instance (Left first simulation in Fig. 10) of the correlation between variables, with high absolute correlations marked in red.

ρ = E[(ξ -ξ)(ζ -ζ)] σ ξ σ ζ , (33)
It is clear from the table that each component in the feature vector relates strongly to one pose parameter. We further demonstrate the correlation in Fig. 12, where we plot the evolution of object poses and feature components. Note that s 2 and θ are negatively correlated. The slower convergence could be a result of the fact that, in contrast with image moments, here the extracted features and object pose are not completely decoupled. Yet, the main contribution of our method is that it can be directly used for both rigid and deformable objects. Therefore, it can be expected to be slower than methods specifically designed for rigid objects (such as image moments).

Table 3: Correlation ρ between s 1 , s 2 , s 3 and x, y, θ.

x y θ s 1 -0.2819 -0.3343 0.9887 s 2 0.2607 -0.8547 -0.0465 s 3 0.9230 0.3629 -0.1426

Figure 12: Progression of the auto-generated feature components (row 1, 3, 5: s 1 , s 2 , s 3 ) vs. object pose (row 2, 4, 6: x, y, θ). We have purposely arranged the variables with high correlation with the same color.

Experiments

Figure 13 outlines our experimental setup. We use a KUKA LWR IV arm. We constrain it to planar (n = 3) motions δr, defined in its base frame (red in the figure): two translations δx and δy and one counterclockwise rotation δθ around z. A Microsoft Kinect V2 observes the object3 . A Linuxbased 64-bit PC processes the image at 30 fps. In the following sections, we first introduce the image processing for contour extraction, then present the experiments.

Image processing

This section explains how we extract and sample the object contours from an image. We have developed two pipelines, according to the kind of contours (See Fig. 14): open (e.g., representing a cable) and closed. We hereby describe the two. contour. The contour is then re-sampled using Algorithm 2. The starting point and the order of the samples is determined by tracking the grasping point (red dot in Fig. 16d) and the centroid of the object (blue dot). We obtain the vector connecting the grasping point to the centroid. Then, the starting sample is the one closest to this vector, and we proceed along the contour clockwise. Therefore, the extracted contour is ordered clockwise and always starting from the same point. This solves the contour ambiguity for symmetrical objects. p O (h) = p curr 21: end if

Vision-based manipulation

In this section, we present the experiments that we ran to validate our algorithms, also visible at https://youtu.be/gYfO2ZxZ5KQ. To demonstrate the generality of our framework, we tested it with:

• Rigid objects represented by closed contours,

• Deformable objects represented by open contours (cables),

• Deformable objects represented by closed contours (sponges).

We carried out different experiments with a variety of initial and target con- We obtain the target contours by commanding the robot with predefined motions. Once the target contour is acquired, the robot goes back to the initial position, and then should autonomously reproduce the target contour.

Again, we set the number of features k = n = 3, and use K = 50 samples to represent the contour c. We set the window size M = 5, both for obtaining the feature vector s and the interaction matrix L. We choose the control gain to be 0.01. A larger control gain may result in faster convergence, but could also lead to oscillation. The local target threshold is set to 0.8. A higher threshold will result in closer local target and vice versa. The Tikhonov factor used to ensure numerical stability for matrix inversion, is set to λ = 0.01. At the beginning of each experiment, the robot executes 5 steps of small4 random motions to obtain the initial features and interaction matrix. For all the experiments, we set the same termination condition at iteration i + 1 using ASE defined in (30) such that:

1. ASE i < 1 pixel and 2. ASE i+1 ≥ ASE i .

In the graphs that follow, we show the evolution of ASE in blue before the termination condition, and in red after the condition (until manual stop by the operator).

Figure 17 presents 8 experiments, one per column. Columns 1 -3, 4 -6 and 7 -8 show respectively manipulation of: cable, rigid object and sponge.

The first row presents the full RGB image obtained from Kinect V2. The second and third rows zoom in on the manipulation at the initial and final iterations. We track the end-effector in the image with a green marker for contour sampling. The target and current contours are drawn in red and blue, respectively. Figure 18 shows the decreasing trend of error ASE for each experiment. The initial increase of ASE in the experiments can be due to the random motion at the beginning of the experiments. In general, we found that ASE is more noisy for the closed than for the open contour. This discontinuity is visible in Figures 18c and18d (zigzag evolution). Such noise is likely introduced by the way we sampled the contour. Also, the noise in the contour extraction is more visible on Fig. 18g and18h. The two plots show that our framework can converge to the target in the presence of noisy image processing. When we have false contour data, the value of ASE may encounter a sudden discontinuity. Figure 19 shows examples of these false samples, output by the image processing pipeline. Despite these errors, thanks to the "forgetting nature" of the receding horizon and to the relatively small window size (M = 5), the corrupted data will soon be forgotten, and it will not hinder the overall manipulation task. Yet, the overall framework would benefit from a more robust sensing strategy, as in (Chi and Berenson). Finally, since our framework can deal with both rigid and deformable objects, we tested it in two experiments where the same object (a sponge) can be both rigid (in the free space), and deformed (when in contact with the environment). These experiments require the robot to: 1) move the object, establish contact, 2) give the object the target shape, by relying on the contact. Figure 20 presents these two original "move and shape" servoing experiments with the corresponding errors ASE plotted in Fig. 21. We use a second fixed robot arm to generate the deforming contact. As the figures and curves show, both experiments were successful.

The success of the "move and shape" task is largely dependent on the contact establishment. However, even when the initial contact has some misalignment (see Fig. 20c and20g), our framework can still reduce the ASE to give a reasonable final configuration (see Fig. 20h and Fig. 21b).

Conclusion

In this paper, we propose algorithms to automatically and concurrently generate object representations (feature vectors) and models of interaction (interaction matrices) from the same data. We use these algorithms to generate the control inputs enabling a robot to move and shape the said object, be it rigid or deformable. The scheme is validated with comprehensive experiments, including a target contour that requires both moving and shaping. We believe it is unprecedented in previous research. Our framework adopts a model-free approach. The system characteristics are computed online with visual and manipulation data. We do not require camera calibration, nor a priori knowledge of the camera pose, object size or shape.

The proposed approach has two major limitations: 1. The challenge of extending it to 6 DoF motion, 2. Global convergence cannot be guaranteed. Below, we discuss each limitation and present possible solutions.

An open question remains the management of 6 DOF motion of the robot. Indeed, while the proposed controller can be easily generalized to 6 DOF motions, it relies on a sufficiently accurate extraction of feature vectors from vision sensors. A very challenging task is to generate complete and reliable 3D feature vectors of objects from a limited sensor set, due to partial views of the object and to occlusions. To extend it, the framework should benefit from robust deformation sensing. In addition, since the approach relies on local linear models, we expect that with higher DOF, the algorithm will more likely get stuck in local minima.

The second drawback is that the representation and model of interaction are local. Thus, they cannot guarantee global convergence. In addition, our framework cannot infer whether a shape is reachable or not. This drawback is solvable by using a global deformation model for control. But as we mentioned earlier, a global model usually requires an offline identification phase which we want to avoid. In fact, for different objects, we will need to re-identify the model. There is a dilemma in using a global deformation model. Maybe one of the possible solutions to this dilemma is to have both our method and deep learning based methods run in parallel. While our scheme enables fast online computation and direct manipulation, the extracted data can be used by a deep neural network to obtain a global interaction mapping. Once a global mapping is learned, it can later be used for direct manipulation and to infer feasibility of the goal shape.

Figure 1 :

 1 Figure1: Vision-based manipulation of rigid and deformable objects. For rigid objects (left): control pose (translation and rotation). For deformable objects (right): control the pose, and also shape.

Figure 2 :

 2 Figure 2: Graphic representation of the vision-based manipulation problem, with its two sub-problems, parameterization and control.

Figure 3 :

 3 Figure 3: The block diagram that represents the overall framework.

Figure 4 :

 4 Figure 4: Six trials conducted to test various choices of feature dimension k for a cable. In each sub-figure, the solid red lines are the initial shapes and the dashed black are the shapes resulting from 10 random motions of the right tip (translations limited to ±5% of the length, rotations limited to ±5 • ).

Figure 5 :

 5 Figure 5: Ten distinctive cable shapes generated by large motion: angle variation: [-π 2 , π 2 ], maximum translation: 106% of the cable length.

  Cable manipulation with a single end-effector, moving the right tip. (a): a reachable target, the blue and black lines are the initial and intermediate shapes, respectively, and the dashed black line is the target shape. The red frame indicates the end-effector position and orientation generated by our controller. (b): Adding Gaussian noise with zero mean and 0.01 standard deviation to the shape data with a reachable target. (c): An unreachable target and the final shape obtained with our controller.

Figure 8 :

 8 Figure8: Comparison -for estimating s -of the receding horizon approach (RH, left) and of the Broyden update (right, with three values of β). The topmost, middle and bottom plots show the one step prediction of s 1 , s 2 and s 3 , respectively. In all plots, the dashed red curve is the ground truth from the simulator. The plots clearly show that the receding horizon approach outperforms all three Broyden trials.

Figure 9 :

 9 Figure 9: Manipulation of a rigid object with a single end-effector (red frame). The initial, intermediate and target contours are respectively blue, solid black and dashed black. Note that in both cases, our controller moves the object to the target pose.

Figure 10 :

 10 Figure 10: From an initial (red) pose, we generate 10 (dashed blue) random motions of a rigid object. This figure shows multiple examples of different rectangular rigid objects. where ξ and ζ are two variables with expected values ξ and ζ and standard deviations σ ξ and σ ζ . An absolute correlation |ρ| close to 1 indicates that the variables are highly correlated. All the simulations in Fig. 10 exhibit similar correlation between the computed feature vector and the object pose.

Figure 11 :

 11 Figure 11: Evolution of ASE of the simulated rigid object manipulation using our method against image moments. Left: simulation in Fig. 9, Right: simulation in Fig. 9.

  Figure 15: Image processing steps needed to obtained the sampled open contour of an object (here, a cable).

Figure 16 :p

 16 Figure 16: Image processing for getting a sampled closed contour: (a) original image, (b) image after thresholding and Gaussian blur, (c) extracted contour, (d) finding the starting sample and the order of the samples.

  tours and camera-to-object relative poses. The variety of both geometric and physical properties demonstrates the robustness of our framework. The variety of camera-to-object relative poses shows that-as usual in image-based visual servoing[START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF])-camera calibration is unnecessary. The algorithm and parameters are the same in all experiments; the only differences are in the image processing, depending on the type of contour (closed or open, see Sect. 6.1).

Figure 17 :Figure 18 :

 1718 Figure 17: Eight experiments with the robot manipulating different objects. From left to right: a cable (columns 1 -3), a rigid object (columns 4 -6) and a sponge (columns 7 and 8). The first row shows the full Kinect V2 view, and the second and the third columns zoom in to show the manipulation process at the first and last iterations. The red contour is the target one, whereas the blue contour is the current one. The green square indicates the end-effector.

  Figure 18: Evolution of e i at each iteration i, for the 8 experiments of Fig. 17. The black dashed lines indicate the threshold ASE = 1 pixel. The blue curves show e i until the termination condition, whereas the red curves show the error until manual termination by the human operator.

Figure 19 :

 19 Figure 19: False contour data from the image can cause noise in ASE.

Figure 20 :Figure 21 :

 2021 Figure20: Two "move and shape" experiments grouped into two rows. The target contour (red dotted) is far from the initial one. This requires the robot to 1) move the object, establish contact with the right -fixed -robot arm, 2) give the object the target shape, by relying on the contact. The first column shows the starting configuration, the second column presents the contact establishment, and the third column zooms in to show the alignment. The last column shows the final results.

Table 2 :

 2 Explained variance Υ(k) computed with large motion.

In the SVD algorithm, the first k columns correspond to the k largest eigenvalues of matrix Σ.

The notion of fast or slow varying depends on both the speed of manipulation, and on the objects deformation characteristics (which affect the rate of change in shapes) with regard to the image processing time.

We only use the RGB image -not the depth -from the sensor.

The notion of small is relative, and usually dependent on the size of the object the robot is manipulating. Refer to Sect.

5.2 (especially Fig.4) for a discussion on this.

Acknowledgement This work is supported in part by the EU H2020 research and innovation programme as part of the project VERSATILE under grant agreement No 731330, by the Research Grants Council (RGC) of Hong Kong under grant number 14203917, and by the PROCORE-France/Hong Kong RGC Joint Research Scheme under grant F-PolyU503/18.

Open contours

The overall pipeline for extracting an open contour is illustrated in Fig. 15 and Algorithm 2. On the initial image, the user manually selects a Region of Interest (ROI, see Fig. 15a) containing the object. In this ROI, we apply thresholding, followed by a morphological opening, to obtain a binary image as in Fig. 15b. This image is dilated to generate a mask (Fig. 15c) used to compute the new ROI for the following image. Figure 15e is the object after a small manipulation motion and 15f shows the mask (in grey color)

which contains the cable. The OpenCV findContour function is applied to binary image, then two contours are extracted based on the two known ends of the cable, both are re-sampled (with same value of K) using Algorithm 2 and finally merged (by interpolation, for each sample, between the two contours' corresponding point) into the uniformly sampled open contour c (see Fig. 15d, where the green box indicates the end-effector).

Closed contours

The procedure is shown in Fig. 16. For an object with uniform color (in the experiment blue), we apply HSV segmentation, followed by Gaussian blur of size 3, and finally the OpenCV findContour function, to get the object