Multiscale Procrustes-Based 3-D Shape Control

The shape control problem involves deforming a nonrigid object into a desired shape. In this article, we present a new 3-D shape control method that can handle large deformations of texture-less objects. In particular, we use functional maps to generate surface points' correspondences between the current and the desired shape that serve as input for our Procrustes-based multiscale control strategy. The main contribution is a multiscale analysis with a relaxed rigidity assumption that defines scales based on geodesic distances to the grippers. Our method considers that objects may deform at different scales depending on the gripper actions applied to them. For this reason, our multiscale analysis computes shape control actions defined according to gripper influence on particularly convenient scales. We present stability analysis and discuss several simulations and experiments.

Multiscale Procrustes-Based 3-D Shape Control Ignacio Cuiral-Zueco , Graduate Student Member, IEEE, and Gonzalo López-Nicolás , Senior Member, IEEE Abstract-The shape control problem involves deforming a nonrigid object into a desired shape.In this article, we present a new 3-D shape control method that can handle large deformations of texture-less objects.In particular, we use functional maps to generate surface points' correspondences between the current and the desired shape that serve as input for our Procrustes-based multiscale control strategy.The main contribution is a multiscale analysis with a relaxed rigidity assumption that defines scales based on geodesic distances to the grippers.Our method considers that objects may deform at different scales depending on the gripper actions applied to them.For this reason, our multiscale analysis computes shape control actions defined according to gripper influence on particularly convenient scales.We present stability analysis and discuss several simulations and experiments.

I. INTRODUCTION
A UTOMATING object deformation processes has multi- ple applications, such as manipulation-related tasks performed in industrial processes, surgical procedures, or home robotics.The wide variety of tasks and object types hinders the process of classifying and contextualising control strategies.This article focuses on the shape control task, which is defined in surveys of the deformable object manipulation literature (e.g., [1], [2], and [3]).Within the deformable object manipulation context, these surveys propose different criteria for the classification of objects, control strategies, or perception methods.In [2], a deformable object classification is proposed.Based on physical and shape criteria, Sanchez et al. [2] classified objects into cloth-like, linear, planar, and solid objects.Within the manipulation tasks, they include the shape control problem for planar objects and subdivide it into single point and multiple points shape control.Survey [3] includes a classification of objects according to their spatial dimensions and categorizes them into 1-D, 2-D, and 3-D objects.Regarding the deformable object perception, Herguedas et al. [1] classified the perception of deformable objects into: force-based, vision-based, and both vision and force-based perception.Focusing more on visionbased sensing, Yin et al. [3] included control-related concepts, such as state estimation, model, state-template, or parameter identification in its perception section.

A. Related Work
Some shape control methods rely on the use of Jacobian matrices that model the object deformation dynamics in a local manner.This is the case of [4], in which the stretch limits of the objects and gripper collisions are considered.In [5], the Jacobian matrix is estimated from visual features and used to control the feedback points of elastic objects.The authors in [6] and [7] exploited the as-rigid-as-possible deformation model (ARAP [8]) in order to achieve shape control.Both approaches are validated through several experiments.Another Jacobian computation approach is presented in [9], where the Jacobian matrix parameters are estimated with the use of truncated Fourier series of the 2-D object's contour (i.e., the object's 1-D closed contour embedded in 2-D).The method is validated in experiments that involve an active and a passive gripper.Using Fourier series as well, dual-arm flexible cable manipulation experiments are presented in [10].Cable manipulation is also tackled in [11], [12], and [13].In the former, the main connectivity of the cable's free configuration space is captured by a simplified precomputed descriptor that allows to numerically solve the cable's ordinary differential equation (ODE) at a reduced cost.The method in [12] presents a B-spline-based model that relies on 3-D tracking wires with the use of a particle filter.The system proposed in [13] models collisions, contacts, and frictions between cables and a plane in the workspace.
Making use of an adaptive deformation model, experiments involving materials, such as foam, meat, or plastic, are carried out satisfactorily in [14].In [15], a principal component analysis is applied to the 2-D contour of the object's silhouette.In [16], real experiments involving large isometric deformations on planar objects, are carried out.This method uses monocular perception and a shape-from-template-based algorithm.Image contour moments are used to define a sliding control strategy to control the shape of objects that range from soft to rigid (articulated) [17].Their method is validated with a stability proof and experiments with a dual-arm robot setup.
In [18], a 3-D deformable object servo control based on a Gaussian process regression online learned model is presented.The method in [19] focuses on deformable object transport Fig. 1.Control scheme of the shape control method.A 3-D object mesh is extracted from the scene's 3-D reconstruction.Using functional maps, the current shape mesh is matched to the reference (or target) shape mesh.A multiscale analysis is then performed making use of a local rigidity measurement of how rigidly the object behaves at each scale.The multiscale analysis provides an optimal scale value that serves as input for our control law.
and introduces a consensus-based deformation model for the manipulation of broad flexible objects with the use of a large number of manipulators.Two heuristic methods and a neural network for shaping nonprehensile materials with plastic deformation characteristics (not elastic) are introduced and compared in [20].They acquire desired shapes by pushing kinetic sand with adapted robotic arms.For more state of the art publications, we refer the reader to the special issue [21], which collects a variety of publications focusing on the manipulation of deformable objects.

B. Proposed Method and Contributions
In the deformable object manipulation literature, some methods focus on specific tasks (e.g., object cutting [22]) or specific object types (e.g., flexible PCB manipulation in [23]).In this article, we propose a method for shape control that tackles the manipulation of texture-less objects that may undergo large 3-D deformations (i.e., object's curvature, length, area, or volume can globally vary in 3-D space).Note that we focus on shape control, meaning that object transport (e.g., considering error with a rigid component, such as an arbitrary translation and rotation of the object) is outside our scope.The goal is to define a control system that allows to manipulate the shape of a deformable object toward a desired target shape with the use of robots and vision sensors (see Fig. 1).
Objects are assumed to be large strain (i.e., to have low Young's modulus) and we assume there is no information about their specific physical properties (density and stiffness).We focus on objects that present certain rigidity such that gravity or inertia do not dominate their behavior (e.g., clothes are not considered).Our assumptions, include proper object perception (no significant occlusions) and knowledge of the initial gripper configuration.While better results might be achieved with appropriate gripper configuration, our method does not rely on the assumption that the grippers are ideally positioned on the object.The effects of gravity and inertia are disregarded given the assumption that deformations are gradual and slow.As the objects may be texture-less, our method does not rely on visual descriptors.
Regarding the control strategy (see Fig. 1) in Section II, we propose a multiscale analysis to determine the scale on which our novel control strategy should be focused at each time instant during the deformation process (scale defined as the topological distance at which surface points lie from a gripper).We do so under what we defined as the local rigidity assumption.We validate the assumption by means of a multiscale Procrustes analysis, i.e., we measure the extent to which the object can be considered to be moving rigidly at each scale.Using this information, we define a 6-DoF control action for each one of the grippers that are manipulating the object.We provide a stability analysis (see Section III) and several simulation and experiments (see Section IV) that validate our proposed system.
In the following, we will discuss and analyse the main contributions and contextualize them in the literature.
1) We define a novel shape control framework for 3-D deformations of both planar and volumetric objects.Our method is the first 3-D shape control system that, with the use of functional maps, performs a holistic shape analysis.That is, thanks to the novel application of [24] in a real shape control setup, we define our 3-D control strategy and shape error by analysing all of the object's perceived geometry (regardless of its visual texture) and consider deformations at different scales.Most existing methods tackle 1-D (linear) objects (e.g., [12], [13], [25], [26], and [27]) or confine their analysis to 2-D contours (e.g., [9], [15], [17], and [28]).Few methods address 3-D deformations of planar and/or volumetric objects, e.g., [14], [18], and [29].Rather than addressing holistic shape control, these methods present task-oriented approaches.For example, they base their shape analysis and errors on a few discrete features, such as a segment's curvature, the object's centroid, or the position of a reduced number of feature points.2) We prove local asymptotic stability of our control system and further validate the method with experiments.Related works also study local stability [14], [18]; other works, such as [29] do not provide stability analysis nor guarantees system convergence.
3) We validate our proposed model using real data from our experiments and compare it against other existing baselines.The comparisons show that our model is suitable for shape control given its proper balance between accuracy and computational time cost.4) Our proposed control method does not require apriori information or initial exploration of the object's behavior.Although online updated, methods such as [14] or [18] require an initial exploration of the object behavior and/or a random initialisation before converging to what they refer to as a good enough initial model.

II. MULTISCALE PROCRUSTES SHAPE CONTROL
In this section, we present a shape control strategy that, for each gripper g = 1, . .., G, defines a 6-DoF gripper control action expressed as a rigid transform U g ∈ SE(3) (composed by translation action in u g ∈ R 3 and rotational action in Euler angles ω ∈ R 3 ).

A. Procrustes Operator and Shape Error Metric.
Functional maps [24], based on the spectral analysis of the mesh through the Laplace-Beltrami operator, allow us to obtain a point-to-point match between current shape and target shape mesh nodes.We denote the current shape mesh nodes positions by x m ∈ R 3 , m = 1, . .., M.These vectors are stacked in matrix X ∈ R 3×M .Each current shape point x m has an associated (matched) target point y m ∈ R 3 .These target points are columnwise stacked in Y ∈ R 3×M .
For simplicity of notation, we define the Procrustes operator (T, d P ) = P(X, Y).This operator encloses the orthogonal Procrustes problem as it takes the two column-to-column matched sets of point coordinates X, Y, and returns their Procrustes distance d P and the rigid transform T(X, Y) ∈ SE(3) that minimizes such distance Matrices X ∈ R 3×M and Ȳ ∈ R 3×M stack the column-wise mean x, ȳ (i.e., the centroid) of matrices X, Y. Matrix R ∈ SO(3) is the rotation component of T and t = ȳ − Rx is the translation component.We can apply the Procrustes operator P(X(t), Y) to obtain shape error which measures how similar shapes are (e(t) = 0 when two shapes are identical).The goal of our control strategy is to reduce the error metric e(t).Before initiating our control strategy, we apply T −1 (t 0 ), obtained from P(X(t 0 ), Y), to our target shape as to bring it closer to our current shape in the 3-D embedding.

B. Local-Rigidity Behavior (LRB) Hypothesis
Consider we were only focusing on the surface points that lie within a topological distance r from a gripper g, i.e., a set of M g object points X g (t, r) ∈ R 3×M g defined by points x m (t, r) ∈ Ω g (t, r), where Ω g (t, r) is the object's surface open domain defined by a geodesic ball of radius r centred at gripper g.Suppose the object's rigidity (unknown for us) allows for points X g (t, r) to move on a rigid manner under small gripper transforms H g , i.e., H g ≈ I 4 (being I 4 the 4 × 4 identity matrix), and the rest of points in X(t), i.e., x m (t) / ∈ Ω g (t, r), remain unaffected by H g .We denote this rigid behavior as local-rigidity behavior (LRB).In this scenario, one could benefit from a Procrustes analysis P(X g (t, r), Y g (r)), where Y g (r) ∈ Y are the points matched to those of X g (t, r).The transform T g (t, r) from the Procrustes analysis can be used to define a incremental transform where H g (t, r) belongs to the geodesic path (in SE(3)) defined from I 4 toward T g (t, r).This path is parameterized by time step Δt ∈ R, Δt ∈ [0, 1] which, when taking low values, generates H g (t, r) that meet the small action requirement for local-rigidity behaviour (LRB).Therefore, the rigid error reduction of the subset X g (t, r) can be performed with actions U g (t) = H g (t, r), and thus, the global Procrustes residual of the whole set X(t) with respect to Y, i.e., e(t), can be reduced too.Note that, as Δt → 0, (3) defines the state equation of X g (t, r) where h denotes homogeneous coordinates and T operates on the product integral generating a time-ordered product which defines the time derivative of points when they are only affected by gripper g.

C. Relaxed LRB Analysis
In the diminishing rigidity concept introduced in [4], an exponential decay of the material's rigidity with respect to gripper positions is assumed.However, an object may present diverse and time-varying behaviors depending on its shape and/or deformation state (e.g., a discontinuous rigidity function as in a mechanism).Our method does not assume any particular rigidity decay function on the object.Rather, our proposed relaxed LRB assumption allows us to evaluate and quantify on which scale (or topological distance) gripper actions are more effective in reducing the shape error.
Local-rigidity behavior (LRB) is certainly met by points grabbed by the grippers (assuming grasping stability).However, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the rest of the object points will most likely undergo deformations and thus not present LRB.For this reason, we base our control strategy on a relaxed assumption of LRB, i.e., we make use of a multiscale analysis that quantifies how close to the LRB our object behavior is for each analyzed scale.In order to perform a multiscale analysis, we establish scale r ∈ [r 0 , R(t)] being r 0 the gripper's size and R(t) the largest topological distance that can be found in the object.Our analysis quantifies the extent to which sets X g (t, r) behave rigid-like under any action U g (t).If actions U g (t) affected X(t) at scale r with ideal LRB and assuming linear action superposition (given object's isotropy and homogeneity), we could estimate the resulting shape points X(t, r) as where function δ g (m, r) allows to disregard actions for points Using the Procrustes analysis P(X(t), X(t, r)), we can obtain a measure w(t, r) of how much the object presents LRB at each scale r (i.e., at each topological distance r from the grippers) when undergoing gripper actions U g (t): Measure w(t, r) ∈ (0, 1], w(t, r) = 1 when the LRB is fully met (i.e., d P (X(t), X(t, r)) = 0).Parameter β > 0 allows to modify the relaxation of the LRB assumption.Lower values of β imply a more relaxed rigidity assumption (e.g., if β = 0.1, almost every w(t, r) ≈ 1, and thus we consider every set X g (t, r) ∀r ∈ [r 0 , R(t)] to move rigidly, even if they do not).We propose using β ≈ 1 × 10 4 , which implies a conservative assumption of LRB.

D. Procrustes-Based Locally Optimal Scale Estimation
Note that w(t, r) also quantifies the effectiveness with which the rigid error of X g (t, r) can be reduced by means of incremental transforms H g (t, r) as defined in (3) (larger w(t, r) implies more effectiveness).With this information, we seek to define gripper actions U g (t) such that they better benefit the global error reduction e(t).In order to define U g (t), we propose analyzing scenarios in r × r ∈ R × R : r, r ∈ [r 0 , R(t)].These scenarios constitute an estimation of the object evolution if it was affected by actions U g (t) = H g (t, r) (defined at scale r) but presented ideal LRB at scale r .Each estimation X(t, r, r ) is defined as In order to perform our analysis, we define an error increment estimation surface Surface ê(t, r, r ) constitutes a continuous surface that provides an insight on the effectiveness and risks in reducing error e(t) by means of actions U g (t) = H g (t, r) (generated under the ideal LRB assumption).We can incorporate our knowledge on how much the LRB is met at each scale r (i.e., w(t, r )) and define an error-reduction effectiveness q(t, r) ∈ R which, for each estimated transform H g (t, r), yields In (11), ideal error increment estimations ê(t, r, r ) are weighted by the error reduction effectiveness w(t, r ).In particular, q(t, r) estimates the global error increment that each H g (t, r) (3) can generate.The logical choice of r when defining gripper actions U g (t) = H g (t, r) would be the one that ensures the largest error reduction, i.e., U g (t) = H g (t, r * (t)), r * (t) = arg min r (q(t, r)).However, w(t, r) (present in the definition of q(t, r)) is computed based on the effects of previous actions defined according to U g (t) = H g (t, r * (t)).For q(t, r) to be reliable, updates in r * (t) should take place in the locality of r * (t).For this reason we define r * (t) as which updates r * (t) in the direction of the ∂/∂r component of the gradient for a given time instant t and thus generates a locally optimal r * (t).
Note that (12) requires q(t, r * (t)) to be continuous differentiable with respect to r.This implies e(t), w(t, r ), and ê(t, r, r ), and thus, Procrustes transforms T g and distances d P , should be continuous differentiable with respect to r.The residual d P is continuous and differentiable as it constitutes a metric in shape space [30].On the other hand, showing that the Procrustes optimisation result T g is continuous and differentiable with respect to r requires more development that can be found in Appendix A.
We now provide some intuitions on impact of β in (8) on (12).When setting a low value for β (e.g., β = 0.1), we assume that our LRB hypothesis holds even for actions and scales that lie far from our current actions U g (t) and scale r * (t).Consequently, surface q(r, t) allows for a more unrestricted evolution of r * (t), which could potentially lead to undesired system behaviors, such as slower performance or, in the case of extremely low values of β, to larger final error.On the other hand, a high β (e.g., β = 1 × 10 4 ) leads to a conservative q(r, t) surface that disregards estimates that lie far from the current actions and locally optimal scale.This conservative approach generates resistance to large changes in r * (t), confining its movement to regions where the LRB hypothesis has been validated through measurements across iterations.Very large β values will not compromise the system's effectiveness, but they will lead to slower convergence as r * (t) evolves more conservatively.

E. Control Strategy
Our control strategy makes use of the Procrustes action defined at scale r * (t), i.e., T g (t, r * (t)).We defined our control law Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
as the new term in the time-ordered product of (4) This results in the state equation of X(t) (for each individual x m (t) ∈ X(t)) with r * (t) updated as in (12).Note that the update rule in (12) needs an initial value of r * (t 0 ) and w(t 0 , r ).We propose w(t 0 , r ) = 1 ∀r , which is equivalent to assuming equal LRB at all scales r .Our initial estimation of r takes the minimum at the initial time instant, that is, r * (t 0 ) = arg min r (q(t 0 , r)).Note that, for the update of r * (t), the partial derivative of q(t, r * (t)) in ( 12) needs to be evaluated only at r * (t).This avoids the need to compute q(t, r) in ( 11) for all r except for the neighborhood of r * (t), contributing to the cost-effectiveness of our method.Furthermore, the rest of the method relies on matrix operations and singular value decomposition (SVD), further enhancing its low-cost nature, as it will be illustrated in the experiments section.

III. STABILITY ANALYSIS
In this section, we prove that the contribution of each gripper to the error derivative ė(t) presents local asymptotic stability.First, some preliminary concepts are presented.Then, we introduce several lemmas that lead to Theorem 3.4, which concludes the local asymptotic stability of the system.

A. Preface
We use Tr to refer to the matrix's trace.For clarity, the dependence on r * and t has been omitted when it could be easily inferred.Recall super-index h denotes homogeneous coordinates.The transpose of the logarithm of a rigid transform is (see [31, where the rotation term log(R g ) can be obtained from being θ ∈ (−π, π) the angle of rotation induced by R g .Transla- (17)

B. Lyapunov Stability
Note that the Procrustes residual is invariant to changes in the frame of reference, as transformations T g can be expressed in any frame of reference.For convenience, we defined G inertial frames of reference F g (t) with axes always aligned with the global reference axes.Although their orientation does not change, frames F g (t) present a changing position xg (t), being xg (t) the centroid of X g (t, r * (t)).This is especially convenient as F g (t) implies xg = 0 ∀ t (for that particular reference).For now on, subindex g implies point coordinates and transforms are defined from frame F g (t).
Lemma 3.1: Tr(X h g (X h g ) log(T g )) = 0. Proof: We can decouple this term into the translation and rotation terms enclosed in (15) Tr where we considered t = V g t g = V g ȳg and, in the last step, trace invariance with respect to matrix transposition and cyclic permutations.The second right-hand side term in the last step of ( 18) is zero since frame of reference F g (t) ensures xg = 0 ∀t and the first term from the last step involves the trace of a symmetric matrix times a skew-symmetric matrix and thus is always zero.Lemma 3.2: −Tr(Y g X g log(R g )) ≤ 0 ∀t and θ ∈ (−π, π).Proof: Considering log(R g ) = − log(R g ), we use ( 16) to obtain where A, B F = Tr(A B) is the Frobenius inner product.We define Y g = Y g − Ȳg ( Ȳg ∈ R 3×M g stacks mean vector ȳg ).Matrix R g (t, r * (t)) is obtained from the Procrustes analysis P(X g , Y g ) in ( 1) through the singular value decomposition of matrix The optimal rotation (i.e., the minimizer of the Procrustes optimisation problem) can be computed as R g = Q g W g .Premultiplying and postmultiplying all four bracket terms in (19) by Q g and W g , respectively, (both are unitary matrix, so the result is invariant to them), we obtain Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Proof: The second term on the right-hand side of the definition of V g in ( 17) is a skew-symmetric matrix and thus can be neglected (Tr(x Ax) = 0 for all A skew-symmetric matrices).The third term on the right-hand side in (17) when using (16) and applying Rodrigues' rotation formula to matrices R g yields where K(θ) is the skew-symmetric matrix from Rodrigues formula and c(θ) ≤ 1 for θ ∈ (−π, π).With these considerations where in the third line, we applied the fact that K 2 is a symmetric matrix and the eigenvalues of a squared 3 × 3 skew-symmetric matrix are always {−1, −1, 0}.Note how ( 22) is always ≤0 for θ ∈ (−π, π) as c < 1 for θ ∈ (−π, π).Theorem 3.4: The system error e(t) (2) is locally asymptotically stable with control law (13).In particular, scalar function F is Lyapunov and presents a derivative Vg ≤ 0, ∀t and θ ∈ (−π, π).
Proof: V(t) > 0, ∀ X(t) = Y and V(t) = 0 only when X(t) = Y.We consider that gripper-unaffected points present a zero derivative, i.e., dX(t)/dt = 0 for x m (t) / ∈ Ω g (t).We also consider that gripper-affected points' derivative is dX h g /dt = log(T g (t, r * (t))X h g (t, r * (t)), i.e., deriving (14) as in (5).Therefore, the derivative of V(t) yields where Applying Lemma 3.1, the first right-hand side in ( 24) is zero and, therefore, Vg yields where, for the second right-hand side term, we considered t = V g t g = V g ȳg and trace invariance with respect to matrix transposition and cyclic permutations.Considering Lemmas 3.2 and 3.3, both right-hand side terms in (25) are ≤ 0 ∀t and θ ∈ (−π, π) and thus the proof is completed.Note that computation of r * relies in the relaxed assumption of LRB.For this reason, the stability analysis remains local since there is no guarantee for a global LRB.Remark 3.5: Note that, as deformable objects constitute under-actuated systems, the local stability in Theorem 3.4 does not guarantee a zero-valued error residual.However, it does ensure convergence of the closed-loop system toward a state with a zero error derivative.

A. Simulations
We performed several simulations using the ARAP [8] deformation model.The simulation results are presented in Fig. 2. Within this figure, each column presents the results of a simulation.We normalized the color map depicting values of q(t, r) so that negative values are blue, positive values are yellow, and zero (or close to zero) values are green.One intuition for interpreting surface q(t, r) is to think of each scale r ∈ [r 0 , R(t)] (at a given instant t) as an available choice of an action to be performed under the relaxed LRB assumption.The value of q(t, r) (11) constitutes the error derivative estimation that each action will cause.Our method seeks large error reduction by choosing actions r * (t) (12) that present large negative values of q(t, r * (t)).
The first simulation in Fig. 2 constitutes a challenging case given its antisymmetry.This is especially noticeable in the last plot where surface q(t, r) presents estimations of positive error derivatives at large scales right at the initial configuration (t = 0, yellow tones around r = 7 [cm]).In other simulations, positive derivative estimations begin to appear later during the control process.Nonetheless, our initial choice of r * (t 0 ) ensures q(t, r * ) < 0 and the control task is performed properly.
The simulation in the second column in Fig. 2 constitutes a paradigmatic example of how our multiscale analysis works.Note how our strategy is able to prioritize larger elements of the object (i.e., larger r * (t) values) when their contribution to the error term is large.The gripper positioned at the small appendix will collaborate in reducing the main beam error.Whenever the main beam error is reduced, the focus shall return to the reduction of the appendix' error (note the scale transition of r * (t) between t = 5 and t = 10).The third column in Fig. 2 involves a pure twist.This can be challenging given singularities of the opposed rotation actions.However, actions are defined symmetrically and the error is smoothly reduced.The fourth column in Fig. 2 constitutes an example of pure bending process.Additional simulations are shown in the attached video.Fig. 2. Four examples, one per column, are presented.First, three rows show the initial (X(t 0 )), final (X(t)), and target (Y) shapes.Each gripper is represented with a gray cube.On the object surface, the surface point matching is represented with a color map.The fourth row contains the shape error plot (2) and a distribution of the object's relative deformation (positive values when the object stretches and negative values if the object is being compressed).The fifth and sixth rows are the action plots of the grippers: translation vector (i.e., u g ) and rotation (Euler angles w) components, respectively.The RGB color code indicates x, y, and z components for translation and roll, pitch, and yaw components for rotation.In the action-related plots (i.e., fifth and sixth rows), each line style refers to a particular gripper.The last row represents surface q(t, r) and the evolution of r * (t) through time (red line).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 3. Four experiments, from left to right, that constitute analogous deformation cases to those of Fig. 2. First block shows several snapshots of the experiments and the target shape (Y).Each gripper is represented with a grey cube.The surface point matching is represented with a color map.Next rows contain the shape error, the translation action vector u g , and rotation action components (Euler angles ω).In the action plots, each line style refers to a particular gripper, and the RGB color code indicates x, y, and z components.Last row represents surface q(k, r) and the evolution of r * (k) (red line).Fig. 4. Comparison of our proposed shape control model against two geometry-based deformation models [8] and [32], and against as a MSD physical-based model.On the top row, for each of the four experiments, the models' estimation error is presented.The estimation error is determined by evaluating the distance distribution between predicted and measured node positions across iterations.On the second row, the time cost of each process is represented in terms of processing frequency.

B. Experiments
We present several experiments (see Fig. 3) with deformation cases analogous to those shown in the simulations.Our setup consists of two ABB IRB120 industrial robots (equipped with pneumatic grippers) and an Intel Realsense D415 RGB-D camera that provides the depth information we use for the 3-D analysis.The objects to be handled are single-colored foam cutouts that favor color-based segmentation.Recall that our method is designed for relatively rigid objects that align with our local-rigidity behavior assumption.Therefore, we cannot guarantee its proper performance with softer objects, such as clothes.The analysis is carried out in MATLAB and the communication with the robots is via TCP/IP.To avoid dropping below the process frequency of 5 Hz with the previously mentioned setup, we propose discretizing r to a set of around 20 scale values.Note that, given the analysis in Appendix A, q(r, t) in ( 12) is continuous differentiable.Therefore, using discrete points to approximate partial derivatives of q(r, t) by means of classic methods (e.g., polynomial regression) holds mathematical significance.Regarding the code implementation, experiments were conducted on an Intel(R) Core(TM) i7-8565U CPU with 1.99 GHz and 16 GB of RAM. 1  All shape control problems tackled in the experiments involve 3-D surfaces.Target shapes have been predefined with robot configurations that ensure that the shape is achievable; this ensures, for example, that the robot is not required to leave its workspace.The setup, the object segmentation and several relevant time instants are shown in Fig. 3.For each of the four experiments 1 Available demo code for the method is provided in https://github.com/nachocz/Multi-scale-Procrustes-based-3D-shape-control.
several acquired RGB images are shown grayscale along with the object segmentation highlighted in blue.The object's 3-D points (corresponding to each time instant) are shown from a slightly elevated viewpoint that provides a better understanding of the 3-D configuration of the object.The grippers are represented with gray cubes, and the target shape Y with the surface mapping (color map on the object) is shown at the bottom of the sequence.The attached video presents several additional perspectives that provide a better insight of the 3-D configurations of the objects.
The first experiment in Fig. 3 involves an antisymmetric bending process along with a light twist.A trend change can be observed on the actions after k = 20, especially on the y and z components of u g (green and blue) and the y components of ω (green).This is a consequence of the evolution of r * (k) during the first 50 iterations: our system focuses first on the bending process as it is responsible for larger error reduction and then tackles the twisting process.
The second experiment consists of unfolding an object from a bent to a straightened configuration.This shape control problem is particularly challenging as, when approaching the target shape (around k = 110) the object buckles, thus, breaking the assumption of small deformations between iterations and generating unreliable estimations of q(k, r).The buckling produces a global shape change on the object, and thus, r * (k) undergoes an abrupt change toward global scale values (i.e., r * (t) = R(t)).Even so, the system shows robustness and manages to continue converging.The buckling process can be better perceived in Fig. 3 and in the top view of the experiments video: between k = 111 and k = 112, the object suddenly changes its curvature with respect to the camera's front view.
The third experiment in Fig. 3 involves a pure-twisting process.As in the analogue simulation, the alignment of the rotation axes can be problematic, as opposite rotation actions are coupled.Our control system generates symmetric rotation actions that evenly deform both sides of the object.Optimal scale r * (t) increases as the object increases its rigidity due to the accumulated torsion: the object's response to pure rotation at larger scales takes some iterations.This is not the case in the analogous simulation, as the ARAP deformation model prioritizes homogeneous deformation along the whole object from the initial time instant.The fourth experiment involves the 3-D bending of a triangle shape.The shape control task is properly managed.However, it is worth mentioning a small error increase at k = 105.This is due to a discontinuity in the surface mapping that infringes the assumption of surface continuity.This lack of map continuity can be appreciated in the last column of Fig. 3, where a color map discontinuity can be observed on the left-hand side of the triangle (k = 105).As a result, peaks (outliers) around k = 105 can be observed in all the plots of the fourth column in Fig. 3.
Fig. 4 shows the analysis of our proposed Procrustes-based deformation model regarding estimation error, computation frequency, and the volume of analyzed data.The estimation error is defined by comparing predicted and measured node positions during iterations.We compared our model against geometrybased deformation models from [8] (used in [6] and [7]) and [32], as well as a mass-spring-damper (MSD) physical-based model.
Our model competes well with these three models, achieving a desirable balance between accuracy and computation time cost.In addition, unlike the time cost of the other three models, our process' time cost also encompasses the cost of computing the control actions.The favorable balance between model accuracy and low-computation time cost makes our model an attractive choice for shape control.Note that these comparisons do not aim to constitute conventional benchmarking: some aspects, such as the type of mesh or the set maximum number of optimization iterations, may affect the performance of different deformation models in distinct manners.This makes it difficult to make fair comparisons.However, seeking fairness in the comparisons, we imposed a maximum limit of 20 iterations on the optimization process for all of the compared models.This approach allows for a reasonable evaluation of model performance in shape control, despite the inherent complexities involved in comparing diverse methods.

V. CONCLUSION
We developed a novel 3-D shape control strategy that does not rely on physical models but rather on a multiscale shape analysis that allows to prioritize those parts of the objects that not only contribute more to the error metric but also can be influenced by the grippers to an estimated extent.We provided theoretical analysis of local asymptotic stability of the control system and performed several simulations and experiments that validate the system showing robustness and proper error convergence.To the best of authors' knowledge, this is the first 3-D shape control framework involving multiscale geometric analysis of textureless objects.Regarding future work, an interesting research line would be decoupling the scale at which each gripper is acting, i.e., independently defining an r * for each gripper and thus (most likely) achieving more error reduction in certain cases.Note that this is nontrivial as the dimensionality of the problem increases by G (i.e., r 1 × . . .× r G ).

A Differentiability of T g With Respect to r
In this appendix, we study continuity and differentiability of T g (i.e., of R g , t g ) with respect to r. Variations in r imply a smooth variation on the domains Ω g being considered for the Procrustes analysis (domains are enlarged or reduced by boundary ∂Ω g ).Regarding the continuous Procrustes problem definition in [30], we can define our analogous space-continuous squared Procrustes residual for X g and Y g as where x and y represent mapped points of the shapes' surface manifolds.For now on, as x are mapped to y with area preserving maps [24], we will use dΩ g to refer to the differential elements on both the current and the target shape manifolds.The optimal translation, in this continuous formulation, is where A(r) is the area of the open domain Ω g and, in the last step, we applied the Leibniz's integral rule.Note that (28) depends on the existence of ∂r .
The optimal Procrustes rotation component R g can be obtained from where M g (r) ∈ R 3×3 is equal to In this continuous formulation, xg = Ω g xdΩ g / Ω g dΩ g and ȳg = Ω g ydΩ g / Ω g dΩ g .Equation (29), and therefore R g is continuous differentiable with respect to r given two conditions: M g being continuous differentiable and existence of M −1 g for r ∈ [r 0 , R(t)].M g (r) is differentiable with respect to r ∀ r > 0 On the other hand, M −1 g exists if there exist at least three nonaligned points x and three nonaligned points y in the domains Ω g of each manifold.This condition is certainly achieved for any 3-D object (planar or volumetric) and r > 0.
g xdΩ g Ω g dΩ g = Ω g (y − R g x)dΩ g Ω g dΩ g (27)and can be differentiated with respect to r to obtain