# Asservissement visuel d'objets déformables à partir du modèle masse-ressort

Fouad Makiyeh

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

## To cite this version:

Fouad Makiyeh. Asservissement visuel d'objets déformables à partir du modèle masse-ressort. Automaticque / Robotique. Université de rennes, 2023. Français. NNT: tel-04394544

HAL Id: tel-04394544

https://theses.hal.science/tel-04394544

Submitted on 15 Jan 2024

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

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
## THÈSE DE DOCTORAT DE

### L'UNIVERSITÉ DE RENNES

ÉCOLE DOCTORALE N° 601
Mathématiques, Télécommunications, Informatique, Signal, Systèmes, Électronique

Spécialité : « Automatique, Productique et Robotique »

Par
« Fouad MAKIYEH »

« Vision-based shape servoing of soft objects using mass-spring model »

Thèse présentée et soutenue à Rennes, le 12/12/2023
Unité de recherche : « Centre Inria de l'Université de Rennes »

### Rapporteurs avant soutenance :

Andrea CHERUBINI Professeur des Universités, Université de Montpellier
David NAVARRO-ALARCON Associate Professor, The Hong Kong Polytechnic University

### Composition du Jury :

Président : Christian DURIEZ Directeur de recherche, Centre Inria de l'Université de Lille

Examineurs : Andrea CHERUBINI Professeur des Universités, Université de Montpellier

David NAVARRO-ALARCON Associate Professor, The Hong Kong Polytechnic University

Ekrem MISIMI Senior Scientist, SINTEF Ocean (Norvège)

Dir. de thèse : Alexandre KRUPA Directeur de recherche, Centre Inria de l'Université de Rennes

Co-dir. de thèse : François CHAUMETTE Directeur de recherche, Centre Inria de l'Université de Rennes

Co-encad. de thèse : Maud MARCHAL Professeur des Universités, INSA de Rennes
<figure><img src="image_3.png" /></figure>
ACKNOWLEDGEMENT

I would like to express my heartfelt gratitude to the members of the jury for their insightful questions and valuable contributions during the defense of this thesis. Your expertise and engagement have added significant depth to the examination process. Special thanks are extended to the rapporteurs for their meticulous comments and constructive feedback on the manuscript, undoubtedly enhancing the quality of this work.

Genuine thanks go to Ekrem MISIMI for funding the GentleMan project and extending a warm invitation to Norway. Collaborating with Ekrem and the project members has been enriching, significantly contributing to the success of this research.

Deep gratitude is extended to Alexandre KRUPA, François CHAUMETTE, and Maud MARCHAL for their guidance and mentorship, instrumental in shaping the direction of this thesis. A special acknowledgment to Fabien SPINDLER for his consistent availability and support with any issues related to the robotic platform.

I would also like to express appreciation to colleagues and Inria members for their collaboration and camaraderie, which made this research endeavor both enriching and rewarding. Additional thanks are conveyed to Bangalore Ravi KIRAN for being a consistent pillar of support, inspiration, and motivation, as well as to my close friends who consistently make my days more enjoyable.

Special and profound thanks go to each member of my family, playing a unique and invaluable role in my academic journey. To my parents, your unwavering love, guidance, and sacrifices have been the bedrock of my perseverance and success, fueling my determination. To my siblings, your support and understanding have been a constant source of strength, making this journey more meaningful.

In conclusion, my gratitude extends to all who have contributed to the realization of this thesis. Your support has been invaluable.
<figure><img src="image_4.png" /></figure>
SYNTHÈSE EN FRANÇAIS

Le domaine de la robotique a connu des avancées dans la manipulation de divers objets, permettant aux robots d'accomplir des tâches complexes. Par exemple, les robots agri-coles sont capables de récolter des fruits et légumes, mais cette récolte pose un défi en raison du risque de destruction de ces produits si des forces de préhension excessives sont appliquées. Cette tâche représente un des nombreux cas où des robots entrent en contact avec des objets souples susceptibles de changer de forme. La manipulation d'objets dé-formables concerne des objets tels que des tissus, des câbles ou bien encore des organes humains qui peuvent subir des déformations non rigides. Cela nécessite des développements particuliers et des stratégies nouvelles pour une manipulation robotique efficace. En effet, la commande des robots pour saisir ou manipuler des objets souples est complexe, et les lois de commande classiques utilisées pour saisir des objets rigides ne peuvent plus être utilisées. Néanmoins, la capacité des robots à manipuler de tels objets souples peut aboutir à de nouvelles possibilités d'automatisation, à améliorer la productivité et renforcer la collaboration entre homme et robot dans divers domaines, notamment dans l'industrie pour la fabrication et l'emballage de produits, dans le domaine médical pour l'assistance aux gestes chirurgicaux et dans de nombreux autres secteurs.

Considérons par exemple le contexte de l'emballage du saumon. Dans ce scénario, le robot doit d'abord détecter le saumon, puis le saisir tout en maintenant sa forme d'origine ou en le déformant pour qu'il prenne la forme souhaitée, facilitant ainsi la tâche de découpe ultérieure. Les systèmes robotiques équipés d'outils de découpe avancés utilisent des algorithmes de découpe précis pour garantir des tailles de filets précises et uniformes, réduisant ainsi les coûts de main-d'œuvre, améliorant la productivité et préservant l'intégrité du produit. Le développement d'un système robotique pour le conditionnement du saumon constitue l'une des principales préoccupations du projet GentleMAN. Ma thèse de doctorat fait partie de ce projet, financé par le Conseil de la recherche de Norvège, dirigé par l'Institut Sintef à Trondheim, et impliquant comme partenaires le Laboratoire d'informatique et d'intelligence artificielle du MIT (États-Unis), l'Université de technologie du Queensland (Australie) et le centre Inria de l'Université de Rennes. Cette thèse se concentre donc sur la manipulation d'objets déformables à l'aide de robots, plus spé-
cifiquement de bras manipulateurs. De plus, cette recherche vise à exploiter les avancées en robotique et en vision par ordinateur pour développer des méthodologies innovantes et des stratégies de commande permettant aux robots de percevoir et de commander de manière adaptative la forme des objets déformables.

Défis

La manipulation d'objets déformables est un domaine en évolution avec divers défis que les chercheurs abordent activement, car ces objets posent d'importantes difficultés en termes de commande et de manipulation.

Tout d'abord, les interactions entre les bras robotiques et les objets déformables impliquent des incertitudes au niveau du contact et des forces de friction. Ces incertitudes peuvent entraîner une instabilité dans la saisie ou provoquer le glissement de l'objet.

Un autre défi réside dans la prédiction du comportement des objets souples. Les matériaux déformables présentent des caractéristiques non linéaires et variables dans le temps, ce qui rend difficile le développement de modèles précis. Une solution possible est d'utiliser des approches basées sur des algorithmes d'apprentissage automatique pour apprendre le comportement de ces objets. Une autre solution envisageable consiste à utiliser des modèles physiques utilisés pour représenter la dynamique de tels objets. Cependant, ces modèles peuvent être peu précis à moins que leurs paramètres soient parfaitement connus.

De plus, estimer la forme, la pose et les déformations internes des objets déformables est difficile en raison de leur nature complexe et non rigide. Les techniques de vision par ordinateur, telles que la reconstruction 3D et les algorithmes d'estimation de forme, peuvent être utilisées pour obtenir des représentations précises des objets déformables. De plus, l'incorporation de simulations basées sur des modèles physiques qui capturent les déformations et les interactions de ces objets peut améliorer les capacités de perception pour les tâches de manipulation.

En outre, la manipulation d'objets déformables nécessite des stratégies de commande robustes et adaptatives. Les techniques traditionnelles de manipulation d'objets rigides ne sont pas directement applicables en raison de la nature flexible des corps déformables. De nouvelles approches de commande qui tiennent compte des changements dynamiques de la forme et de la compliance de l'objet sont nécessaires. Ces méthodes impliquent l'incorporation de capteurs pour estimer et mettre à jour de manière continue les informations sur les propriétés de l'objet ou les conditions de l'environnement. Cela garantit que le
contrôleur peut s'adapter et répondre de manière appropriée au comportement non rigide de l'objet considéré.

En résumé, relever les défis de la manipulation d'objets déformables nécessite une approche multidisciplinaire. En utilisant des modèles basés sur les données, des modèles physiques, l'haptique, des stratégies de commande adaptatives et des techniques de perception avancées, il est possible d'améliorer la précision, la stabilité et l'efficacité de la manipulation d'objets déformables.

## Objectifs et contributions

Cette thèse contribue à l'état de l'art de la manipulation d'objets déformables en introduisant de nouvelles lois de commande basées sur un modèle physique pour la manipulation d'objets souples en temps réel. Plus précisément, nous utilisons un modèle simple qui est le modèle masse-ressort. Bien que ce modèle puisse capturer les caractéristiques essentielles des objets déformables, il peut rencontrer des difficultés à reproduire avec précision le comportement réel de l'objet. Cependant, en incorporant des techniques avancées de perception visuelle en temps réel, il est possible de corriger les décalages entre le modèle et la dynamique réelle de l'objet.

Le manuscrit de la thèse suit la méthodologie suivante :

Dans le deuxième chapitre, nous présentons une revue de l'état de l'art qui inclut différentes approches utilisées pour modéliser un objet souple, estimer les paramètres des modèles physiques utilisés pour décrire sa dynamique, différentes lois de commande utilisées pour sa manipulation, et divers algorithmes employés pour le suivre.

Dans le troisième chapitre, nous présentons une nouvelle méthode basée sur le système masse-ressort pour estimer le déplacement de chaque point du modèle en fonction des mouvements qui lui sont appliqués. Nous désignons par les points manipulés les points du modèle où les mouvements sont appliqués. Nous commençons par établir la relation analytique qui relie les déplacements des points de l'objet aux mouvements successifs des points manipulés. Cette relation forme la base de notre contribution, construite sur le modèle masse-ressort en tenant compte du retard de propagation qu'il introduit. Nous formulons les déplacements des points dans un maillage 1D en réponse à une succession de mouvements appliqués à un ou deux points manipulés, puis nous généralisons cette formulation pour un maillage 3D impliquant plusieurs points manipulés. En se basant sur ces formes analytiques, nous introduisons une loi de commande qui permet de po-
situationner indirectement plusieurs points de l'objet à leur position souhaitée en agissant sur des points manipulés éloignés. Dans ce même chapitre, nous présentons des résultats de simulations qui démontrent la robustesse du système en boucle fermée proposé. Nous prenons en compte des facteurs tels que la résolution du maillage et les incertitudes des paramètres du modèle pour évaluer l'efficacité de la loi de commande. Nous présentons également une étude comparative entre l'approche proposée et deux méthodes de l'état de l'art qui n'utilisent pas de modèle.

* Une vidéo qui présente divers résultats de simulation pour la tâche de positionnement indirecte et qui montre également la robustesse du contrôleur proposé vis-à-vis des incertitudes des paramètres du modèle ainsi que les résultats de comparaison, est disponible ici : [Vidéo-positionnement-indirect].

Le quatrième chapitre se concentre sur la validation expérimentale de la tâche de positionnement indirect avec des objets déformables. Cette tâche repose sur la méthodologie introduite dans le chapitre précédent, où nous utilisons à la fois des approches basées sur la vision (suivi) et des approches basées sur le modèle (loi de commande) pour accomplir le positionnement indirect de points caractéristiques d'objets souples, modélisés par un modèle masse-ressort. L'objectif est de démontrer l'efficacité de cette approche par des expériences réelles. Une première étape consiste à créer un modèle approximé de l'objet souple, puisque celui-ci n'est pas connu à l'avance. Le modèle est ensuite constitué par un maillage 3D décrivant la géométrie de l'objet et par un système masse-ressort pour estimer son comportement physique. Cependant, ce modèle étant une approximation, il peut y avoir un décalage entre le modèle et l'objet réel. Pour corriger cet écart, l'objet est suivi en utilisant une caméra RGB-D, puis réaligné en appliquant des contraintes externes sur l'objet qui dépendent des mesures visuelles. Une fois que le modèle de l'objet est créé et que ses paramètres sont approximés, des résultats expérimentaux sont présentés pour évaluer l'efficacité de l'approche proposée. Cette approche d'asservissement visuel combine à la fois le retour visuel et le comportement dynamique de l'objet qui est prédit par le modèle.

* Les résultats du positionnement indirect d'un point de l'objet ont été présentés dans notre première publication [MAKIYEH et al., 2022], avec une vidéo associée disponible ici : [Vidéo-ICARCV].

* Une autre vidéo présentant divers résultats expérimentaux pour le positionnement
indirect de deux points de l'objet peut également être trouvée ici : [Vidéo-positionnement-indirect].

Le cinquième chapitre présente deux méthodes de commande basées sur le système massressort pour amener des objets déformables à une forme globale désirée. Ces méthodes visent à déformer le contour pour des objets 2D et la surface complète pour des objets 3D. Pour réduire la complexité de la forme, deux approches de réduction de dimension sont proposées : les descripteurs de Fourier et les moments 3D. Ces approches sont évaluées à l'aide de simulations et d'expériences pour démontrer leur efficacité dans la commande de la forme d'objets souples. L'objectif est de permettre l'application d'une forme désirée à ces objets en utilisant un nombre limité de points manipulés.

* Les résultats de notre approche de commande de la forme d'objets en 2D et 3D basée sur les descripteurs de Fourier ont été présentés dans notre deuxième publication [MAKIYEH et al., 2023], avec une vidéo illustrant deux résultats expérimentaux disponibles ici : [Vidéo-IROS].

Le dernier chapitre offre une conclusion sur les approches proposées, ainsi qu'une discussion sur les perspectives à court et à long terme de la recherche future.
<figure><img src="image_5.png" /></figure>
# TABLE OF CONTENTS

## 1 Introduction 15

1.1 Introduction 15

1.2 Challenges 16

1.3 Thesis objectives and contributions 18

## 2 State Of The Art 21

2.1 Visual servoing 22

2.1.1 Basics of visual servoing 22

2.1.2 Classical Image-Based visual servoing 23

2.1.3 Conclusion 27

2.2 Deformation modeling 28

2.2.1 Overview on continuum mechanics 28

2.2.2 Particle-Based Models 31

2.2.3 Position-Based Dynamics 33

2.2.4 Mesh-based Models 35

2.2.5 Conclusion 38

2.3 Deformable object tracking 39

2.3.1 Model-free 39

2.3.2 Model-based 41

2.3.3 Conclusion 47

2.4 Parameter estimation methods 47

2.4.1 Estimation of non-linear heterogeneous FEM parameters 48

2.4.2 Estimation of homogeneous FEM elasticity parameters using a robotic manipulator 49

2.4.3 Estimation of MSM parameters using a FEM reference model 50

2.4.4 Estimation of isotropic co-rotational FEM parameters using a robotic manipulator 51

2.4.5 Estimation of coarse FEM parameters 51

2.4.6 Conclusion 52
TABLE OF CONTENTS

<table><tbody><tr><td>2.5</td><td>Deformable object manipulation</td><td>52</td></tr><tr><td>2.5.1</td><td>Model-free</td><td>53</td></tr><tr><td>2.5.2</td><td>Model-based</td><td>61</td></tr><tr><td>2.5.3</td><td>Conclusion</td><td>67</td></tr><tr><td>2.6</td><td>Positioning of this thesis in relation to the existing literature</td><td>69</td></tr><tr><td><b>3</b></td><td><b>Methodology</b></td><td><b>75</b></td></tr><tr><td>3.1</td><td>Deformable object servoing</td><td>76</td></tr><tr><td>3.1.1</td><td>Simple 1D MSM</td><td>77</td></tr><tr><td>3.1.2</td><td>General 3D MSM</td><td>84</td></tr><tr><td>3.1.3</td><td>Control scheme</td><td>94</td></tr><tr><td>3.2</td><td>Simulation results</td><td>97</td></tr><tr><td>3.2.1</td><td>Object parameters</td><td>97</td></tr><tr><td>3.2.2</td><td>Robustness analysis – 3D Object</td><td>98</td></tr><tr><td>3.2.3</td><td>Actuation study</td><td>102</td></tr><tr><td>3.2.4</td><td>Comparison with two state-of-the-art approaches</td><td>107</td></tr><tr><td>3.2.5</td><td>Effect of the feed-forward term</td><td>114</td></tr><tr><td>3.3</td><td>Conclusion</td><td>115</td></tr><tr><td><b>4</b></td><td><b>Indirect Positioning of Multiple Points on a Soft Object</b></td><td><b>117</b></td></tr><tr><td>4.1</td><td>Deformable object modeling</td><td>118</td></tr><tr><td>4.1.1</td><td>Model 3D mesh generation</td><td>118</td></tr><tr><td>4.1.2</td><td>Model physics parameters estimation</td><td>119</td></tr><tr><td>4.1.3</td><td>Model correction by tracking</td><td>126</td></tr><tr><td>4.2</td><td>Experiments</td><td>128</td></tr><tr><td>4.2.1</td><td>Experimental setup</td><td>128</td></tr><tr><td>4.2.2</td><td>Experimental results</td><td>131</td></tr><tr><td>4.3</td><td>Conclusion</td><td>146</td></tr><tr><td><b>5</b></td><td><b>Shape servoing of a soft object</b></td><td><b>149</b></td></tr><tr><td>5.1</td><td>Problem formulation</td><td>150</td></tr><tr><td>5.2</td><td>Shape servoing of a soft object using Fourier Descriptors</td><td>150</td></tr><tr><td>5.2.1</td><td>3D shape servoing</td><td>150</td></tr><tr><td>5.2.2</td><td>2D shape servoing</td><td>155</td></tr><tr><td>5.2.3</td><td>Simulation results</td><td>156</td></tr></tbody></table>
<table><tbody><tr><td>5.2.4</td><td>Experimental results</td><td>159</td></tr><tr><td>5.3</td><td>Shape servoing of a soft object using 3D moments</td><td>163</td></tr><tr><td>5.3.1</td><td>3D shape servoing</td><td>163</td></tr><tr><td>5.3.2</td><td>Simulation results</td><td>167</td></tr><tr><td>5.3.3</td><td>Experimental results</td><td>170</td></tr><tr><td>5.4</td><td>Conclusion</td><td>171</td></tr><tr><td><b>6</b></td><td><b>Conclusion</b></td><td><b>175</b></td></tr><tr><td>6.1</td><td>Conclusion</td><td>175</td></tr><tr><td>6.1.1</td><td>Indirect positioning of multiple points on a soft object using a simple mass-spring-model</td><td>175</td></tr><tr><td>6.1.2</td><td>Shape servoing of a soft object using a simple mass-spring-model</td><td>176</td></tr><tr><td>6.2</td><td>Short-term perspectives</td><td>177</td></tr><tr><td>6.2.1</td><td>Advancing precision in deformable object manipulation - Exploring advanced control strategies and model refinement</td><td>177</td></tr><tr><td>6.2.2</td><td>Complex scene - multiple objects</td><td>179</td></tr><tr><td>6.3</td><td>Long-term perspectives</td><td>179</td></tr><tr><td>6.3.1</td><td>Multi-Modal sensing to enhance perception and precision in soft object manipulation</td><td>179</td></tr><tr><td>6.3.2</td><td>Sensitive organ manipulation - Human-Machine collaboration</td><td>180</td></tr><tr><td><b>7</b></td><td><b>Appendix</b></td><td><b>183</b></td></tr><tr><td>7.1</td><td>Shape servoing based on 3D moments - Derivation of the interaction matrix and the feed-forward term</td><td>183</td></tr><tr><td><b>Bibliography</b></td><td></td><td><b>199</b></td></tr></tbody></table>
<figure><img src="image_6.png" /></figure>
# INTRODUCTION

## 1.1 Introduction

In recent years, the field of robotics has witnessed remarkable advancements in the manipulation of various objects, enabling robots to perform a wide range of complex tasks. For example, in packaging and assembly lines, robots within manufacturing industries are frequently utilized to grasp objects such as boxes and bottles from one location, and to accurately place them in another. However, during such tasks, these entities can be deformed, altering their shapes and potentially causing them to slip from the robotic hand. Additionally, in agriculture, agricultural robots are capable of tasks such as picking vegetables, though picking tomatoes presents a challenge due to the risk of destruction if excessive gripping forces are applied. Furthermore, in the field of medicine, surgical robots are employed to manipulate surgical tools during operations. However, these tools come into contact with delicate human organs such as the liver, necessitating careful consideration of the applied force to prevent damage.

Automating the manipulation of deformable objects with a robot presents a significant challenge due to their non-rigid and highly variable nature. Deformable object manipulation involves items like fabrics, cables, and biological tissues that can undergo non-rigid deformations, requiring distinct skills and strategies for effective robotic handling. This challenge is especially prominent within the field of robotics, as interacting with deformable objects poses unique difficulties. Nonetheless, the capacity of robots to handle such soft objects can unlock novel possibilities for automation, improve productivity, and enhance human-robot collaboration across various domains, including manufacturing, healthcare, and search-and-rescue operations.

As an illustrative example, let us consider the context of salmon packaging. In this scenario, the robot needs to first detect the salmon, then grasp it while maintaining its original shape or reshaping it into the desired forms that facilitate the subsequent cutting task. Robotic systems equipped with advanced cutting tools utilize precise cutting
algorithms to ensure accurate and uniform fillet sizes, thereby reducing labor costs, en-
hancing productivity, and maintaining product integrity. An example of this pipeline is
presented in Figure. 1.1, which is one of the main concerns of the GentleMAN (299757)
project. My PhD is part of this project, funded by The Research Council of Norway, led
by the Sintef Institute in Trondheim, Norway, and involving partners such as the MIT
Computer Science and Artificial Intelligence Laboratory (USA), Queensland University
of Technology (Australia), and the Inria center at University of Rennes.

This PhD research aims to address the fundamental challenges associated with de-
formable object manipulation using robotic manipulators. By leveraging cutting-edge
advancements in robotics and computer vision, this study seeks to develop novel method-
ologies and control strategies that enable robots to perceive and to control the shape of
deformable objects in an adaptive manner.

<figure><img src="image_7.png" /><figcaption>Figure 1.1: An example of automated robotic salmon fillet cutting, from [Einarsdóttir et al., 2022].</figcaption></figure>

1.2 Challenges

Deformable object manipulation is a dynamic and evolving field with various challenges
that researchers actively tackle, as these objects pose significant difficulties in terms of
accurate control and manipulation.
Firstly, the interactions between robotic grippers and deformable objects involve un-
certainty in terms of contact and friction. Deformable materials often exhibit complex
contact behavior, such as sticking, sliding, and rolling. This uncertainty may result in
grasping instability or object slippage.

Another challenge involves predicting the deformation behavior of soft objects. In-
deed, deformable materials exhibit non-linear and time-varying characteristics, making it
challenging to develop accurate models. One possible solution is to employ data-driven
approaches, such as machine learning algorithms, to learn the behavior of deformable ob-
jects from real-world sensor data. By training models on a diverse set of examples, these
algorithms can capture the underlying dynamics and improve manipulation performance.
Another possible solution is to consider physics-based models. However, these models
may not achieve high accuracy unless the exact physics parameters of the actual soft
objet are known. To address this limitation, incorporating real-time sensing data into the
models, the object deformations can be refined.

Additionally, estimating the shape, pose, and internal deformations of deformable ob-
jects is challenging due to their complex and non-rigid nature. Computer vision tech-
niques, such as 3D reconstruction and shape estimation algorithms, can be employed to
obtain accurate representations of deformable objects. In addition, incorporating physics-
based simulation that capture the deformations and interactions of these objects can
enhance perception capabilities for manipulation tasks.

Furthermore, deformable object manipulation requires robust and adaptive control
strategies. Traditional rigid object manipulation techniques are not directly applicable
due to the flexible nature of deformable bodies. Control approaches that account for
the dynamic changes in object shape and compliance are needed. These methods involve
incorporating sensor feedback to consistently estimate and update information about the
object properties or the environment conditions. This ensures that the controller can
adapt and respond appropriately to these changes.

In summary, addressing the challenges of deformable object manipulation requires a
comprehensive approach. By leveraging data-driven models, physics-based models, tactile
sensing, adaptive control strategies, and advanced perception techniques, it is possible to
improve the accuracy, stability, and efficiency of manipulating deformable objects.
1.3 Thesis objectives and contributions

The objective of this thesis is to develop a physics-based controller for deformable object
manipulation, focusing on real-time applications.

The research aims to investigate the effectiveness of using a simple model that is,
the mass-spring model, for shape servoing of soft objects. While this simple model can
capture the essential characteristics of deformable objects, it may struggle to accurately
replicate the real object behavior. However, by incorporating advanced real-time percep-
tion techniques, it is possible to correct any discrepancies or drift between the model and
the actual object dynamics. By doing so, this research aims to develop a computationally
efficient, interpretable, and uncertainty-aware controller.

This thesis contributes to the state-of-the-art in soft object manipulation by introduc-
ing novel control laws for controlling the shape of deformable objects. Shape servoing
involves driving specific points of the object to desired positions by deforming the object
using robotic manipulators.

To achieve this, the following chapters are included:

* Chapter 2 presents a comprehensive literature review of relevant state-of-the-art. It includes different approaches for modeling a soft object, estimating the parameters of physics-based models, various control laws used for manipulating a soft object, and diverse algorithms employed for tracking it.

* Chapter 3 outlines a novel physics-based method for estimating the displacement of any model point in function of the applied motions on the manipulated points based on the mass-spring model. This method is then employed to design a control law for the deformation task, enabling indirect positioning of object points by acting on different points. Indeed, the proposed control law addresses the challenge of considering the inertia and flexible nature of deformable objects. This chapter also includes simulation results illustrating the success and robustness of the indirect positioning task using the proposed approach. Furthermore, the chapter compares the effectiveness of our approach in achieving the indirect positioning task with two state-of-the-art model-free methods.

- A video presenting various simulation results for the indirect positioning task, including the robustness of the proposed controller against model parameters and the comparison results, can be found here: [Indirect-positioning-video]
* Chapter 4 is dedicated to validating the indirect positioning task introduced in Chapter 3 within an experimental setup, considering various soft objects. To accomplish this, we begin by outlining the different steps for reconstructing a model corresponding to an unknown soft object, where both its geometry and material properties are unknown. However, given that physics-based models are approximations of real object behavior, there may exist a gap between the object model deformation and its actual behavior. To overcome this challenge, we employ a method involving real-time tracking of the object during deformation. Whenever a drift is detected, we rectify the object model by applying external forces, thereby aligning it with the real deformation. Subsequently, we accomplish the indirect positioning task of one or multiple object points to their desired positions. This is achieved by employing our proposed method, which integrates now both vision-based (tracking process) and model-based (control law) approaches. The deformation is achieved using two robotic manipulators and a RGB-D camera.

- The results of the indirect positioning of one feature point were presented in our first publication [Makiyeh et al., 2022], with an accompanying video available here: [ICARCV-video]

- Another video showing various experimental results for the indirect positioning of two feature points can also be found here: [Indirect-positioning-video]

* Chapter 5 introduces two physics-based methods for controlling the complete shape of the soft object. The object shape is defined by the complete set of its surface points in the case of 3D objects and by its 2D contour when dealing with 2D objects. To represent the object shape effectively using a low-dimensional feature vector, we present two distinct approaches. The first approach employs Fourier descriptors, while the second utilizes 3D moments. For each of these approaches, we assess the corresponding proposed control law for the deformation task by conducting a comprehensive evaluation that combines both simulation and experimental results.

- The results of our servoing approach based on Fourier descriptors were presented in our second publication [Makiyeh et al., 2023], with an accompanying video illustrating two experimental results available here: [IROS-video]

* Chapter 6 offers a conclusion on the proposed approaches, as well as discussing short term and long-term perspectives on future research.
<figure><img src="image_8.png" /></figure>
STATE OF THE ART

Deformable object manipulation has emerged as a challenging yet promising field with applications spanning robotics, medical procedures, and virtual reality. Perceiving and manipulating objects with flexible and complex structures, such as cloth, soft tissue, or cables, require precise control strategies that account for their nonlinear and dynamic behaviors.

This research specifically concentrates on the shape servoing of soft objects using robotic manipulators. Over time, researchers have explored various control methods and techniques to address this challenge and enhance manipulation capabilities.

To manipulate soft objects, some approaches first approximate their physical behavior using physics-based models and then develop control laws to guide the robots in accomplishing the task. In these model-based approaches, object models are simulated to derive the control laws, which are then applied to the robots in either open-loop execution or within a closed-loop system using feedback measurements. Open-loop systems typically rely on highly accurate models, while closed-loop systems can adapt to models of varying accuracy, benefiting from the feedback loop to execute the task effectively. When using physics-based models, it is crucial to note that these models depend on physical parameters specific to the deformable object. If these parameters are unknown, they must be estimated. Additionally, these models are still approximations of the real behavior of the object, meaning that the simulated deformation may deviate from the actual deformation of the object. This deviation can be detected by tracking the object, and as a result, continuous correction of any deviations becomes essential if they occur.

In contrast, other approaches rely solely on visual feedback information without the need for complex modeling of soft objects to generate the corresponding control laws. These approaches, known as model-free approaches, leverage vision feedback to identify object relevant features, and formulate a control strategy for closed-loop execution. The robot adjusts its position to align its observed features with the desired ones, corresponding to the intended task. This approach is commonly referred to as a visual servoing
problem.

This chapter provides an overview of the current literature in deformable object ma-
nipulation, highlighting the strengths and limitations of existing approaches. By exploring
the recent advancements, this review sets the stage for the subsequent chapters, which
aim to propose novel model-based control strategies and algorithms to enhance the ma-
nipulation of deformable objects in various practical applications.

We start by providing an introduction to the basics of visual servoing in Section 2.1, a fundamental aspect of robotics control method that utilizes visual data from cameras to direct a robot motion. Then, in Section 2.2, we delve into the various techniques that have been investigated for simulating soft object deformations, which is essential for model-based shape servoing approaches. Following that, we explore in Section 2.3 various tracking methods aimed at minimizing the discrepancy between simulated object deformations and their real deformations. Additionally, we review in Section 2.4 different approaches for physics-based model parameter estimation, which is necessary for model-based shape servoing approaches. Furthermore, we explore in Section 2.5 a wide range of techniques that have been studied for shape servoing of soft objects, encompassing model-based and model-free approaches. Lastly, in Section 2.6, we contextualize this work in relation to the existing literature.

## 2.1 Visual servoing

This section presents fundamental principles of visual servoing, which is a technique that
allows robots to interact with both rigid and deformable objects. Visual servoing control
involves using visual information provided by a visual sensor to control robot movements.
This visual feedback information can be obtained from a camera attached to a robot,
which moves along with the robot (eye-in-hand scenario), or from a stationary camera
observing the robot within the workspace (eye-to-hand scenario).

2.1.1 Basics of visual servoing

The reader can refer to [Chaumette & Hutchinson, 2008] for a detailed presentation of the basics of visual servoing. Visual servoing is a technique that converts visual measurements into a vector of *n* visual features, denoted as *s* ∈ R*n*, enabling the robot to interact with its environment. The difference between the current value of this vector and the desired
value, $\mathbf{s}^*$, is used to calculate the task error, $\mathbf{e}$, as follows:

$$ \mathbf{e} = \mathbf{s} - \mathbf{s}^* \qquad (2.1) $$

The task error is then used to generate control laws to be applied to the robots in order to minimize its norm.

When the variations in $\mathbf{s}$ are solely influenced by the camera motion, then the relationship between the changes in the visual feature vector, $\dot{\mathbf{s}}$, and the camera velocity $\mathbf{v}_c = (\mathbf{v}_c, \omega_c) \in \mathbb{R}^6$, is expressed using an interaction matrix denoted as $\mathbf{J}_s$. Here, $\mathbf{v}_c$ represents the instantaneous linear velocity vector, and $\omega_c$ represents the instantaneous angular velocity vector. In case $\mathbf{s}^*$ remains constant, which induces $\dot{\mathbf{e}} = \dot{\mathbf{s}}$, the expression of the relation between $\dot{\mathbf{e}}$ and $\mathbf{v}_c$ is therefore as follows:

$$ \dot{\mathbf{e}} = \dot{\mathbf{s}} = \mathbf{J}_s \mathbf{v}_c = \mathbf{J}_e \mathbf{v}_c \qquad (2.2) $$

Typically, the error is chosen to exhibit an exponential decoupled decrease, represented as $\dot{\mathbf{e}} = -\lambda \mathbf{e}$, where $\lambda > 0$. This specific selection plays a crucial role in designing the camera control velocity, denoted as $\mathbf{v}_c$, which follows the formulation:

$$ \mathbf{v}_c = -\lambda \hat{\mathbf{J}}_e^+ \mathbf{e} \qquad (2.3) $$

with $\hat{\mathbf{J}}_e^+$ the Moore-Penrose pseudo-inverse of $\hat{\mathbf{J}}_e$, an approximation of $\mathbf{J}_e$.

Visual servoing can be classified into two distinct categories: Image-Based Visual Servoing (IBVS) and Pose-Based Visual Servoing (PBVS). The former leverages the image to extract various features, such as image points, image moments, and image Fourier descriptors, to name a few. On the other hand, the latter employs image data to estimate three-dimensional measurements, such as the pose of an object, which subsequently aids in the servoing process. In this thesis, our focus lies on IBVS since it is related to our chosen methodology.

## 2.1.2 Classical Image-Based visual servoing

When constructing the visual feature vector $\mathbf{s}$ from image points $p_1, p_2, \dots, p_n$, the following procedure is employed. Each image point $p_i$ corresponds to a pair of horizontal and vertical coordinates $(X_i, Y_i)$ obtained from projecting 3D points $P_i$, with corresponding 3D coordinates $(x_i, y_i, z_i)$, expressed in the camera Cartesian frame. The visual feature
vector **s** is then formed as **s** = (p₁, p₂, ..., pₙ).

For instance, when **s** is simply constructed from one image point, i.e **s** = (X₁, Y₁), the interaction matrix is obtained as follows:

$$
\mathbf{J}_{\mathbf{s}} = \begin{bmatrix} \frac{-1}{z_1} & 0 & \frac{X_1}{z_1} & X_1 Y_1 & -(1+X_1^2) & Y_1 \\ 0 & \frac{-1}{z_1} & \frac{Y_1}{z_1} & 1+Y_1^2 & -X_1 Y_1 & -X_1 \end{bmatrix} \quad (2.4)
$$

It should be noted that the matrix in (2.4) depends on the depth information of the
point. Therefore, when using a RGB camera, the depth information is not available,
and an approximation of this matrix needs to be utilized. Multiple approaches for either
estimating the point depth or directly numerically estimating the interaction matrix are
recalled in [Chaumette & Hutchinson, 2008]. We now describe two methods for numer-
ically estimating such a Jacobian matrix, the Broyden algorithm [Broyden, 2000] and
Least-Squares Minimization [de Mathelin & Lozano, 1999].

2.1.2.1 Numerical estimation of the feature Jacobian matrix using the Broy-
den algorithm

The Broyden algorithm is an iterative method used to numerically estimate the Jacobian
matrix of a function *f* [Broyden, 2000]. It is particularly useful when the function is
expensive to compute or when the analytical expression for the Jacobian is not available.
The Broyden algorithm updates an approximation of the Jacobian matrix based on the
changes in the function outputs and inputs.

Let consider a function *f*: ℝ<sup>*m*</sup> → ℝ<sup>*n*</sup>, where *m* is the dimension of the input vector
**q** ∈ ℝ<sup>*m*</sup> and *n* is the dimension of the output **s** ∈ ℝ<sup>*n*</sup>. The goal is to estimate the Jacobian
matrix **J**<sub>*f*</sub>, a *n* × *m* matrix, with its approximation denoted **\tilde{J}**<sub>*f*</sub>. The different steps to
update the Jacobian matrix in the Broyden algorithm are as follows:

1. Initialization: Start by providing $\tilde{\mathbf{J}}_f$, an initial approximation to $\mathbf{J}_f$. This can be a numerical approximation or an identity matrix, depending on the problem and available information.

2. Evaluation: Initialize an initial guess for the inputs, evaluate the function *f* at the initial guess, and compute the corresponding outputs.

3. Iterative Steps: Perform the following steps iteratively at each iteration *k* until convergence:
• Compute the difference between the function output and the previous output:

$$ \delta s = f(q_k) - s_{k-1} $$

• Compute the difference between the current input and the previous input:

$$ \delta q = q_k - q_{k-1} $$

• Update the Jacobian matrix with an update speed α:

$$ \hat{J}_{f_k} = \hat{J}_{f_{k-1}} + \alpha \frac{\delta s - \hat{J}_{f_{k-1}} \delta q}{\delta q^T \delta q} \delta q^T $$

• Update the inputs:

$$ q_{k+1} = q_k - \hat{J}_{f_k}^+ f(q_k) $$

• Update the function outputs:

$$ s_{k+1} = f(q_{k+1}) $$

• Check for convergence. If the change in the inputs and outputs is below a specified tolerance, terminate the iterations.

4. Convergence: Once the iterations have converged, the final approximation of the Jacobian matrix, $\hat{J}_f$, is obtained.

The Broyden update is computationally efficient, and supports incremental updates, making it well-suited for situations where the matrix or its inverse requires adaptation in response to changing data or conditions. However, this approach does have some limitations. For instance, it can be sensitive to numerical stability issues, particularly when dealing with ill-conditioned matrices. Furthermore, it does not guarantee global convergence, or even convergence at all. Indeed, its convergence relies on the choice of the initial guess.
### 2.1.2.2 Numerical estimation of the feature Jacobian using Least-Squares Minimization

The Least-Squares Minimization algorithm is another method used for the numerical estimation of the Jacobian matrix of a function $f$ [de Mathelin & Lozano, 1999]. This algorithm involves formulating the problem as a least-squares minimization, where the goal is to find the Jacobian matrix that minimizes the sum of squared residuals between the actual function outputs and the estimated outputs based on the Jacobian.

Let consider again a function $f: \mathbb{R}^m \to \mathbb{R}^n$. Here is the step-by-step explanation of the Least-Squares Minimization algorithm for estimating the Jacobian matrix $\hat{J}_f$, by $\hat{J}_f$:

1. Initialization: Start by providing $\hat{J}_f$ an initial approximation to $\hat{J}_f$.

2. Evaluation: Initialize an initial guess for the inputs, evaluate the function $f$ at the initial guess, and compute the corresponding outputs.

3. Iterative Steps: Perform the following steps iteratively at each iteration $k$ until convergence:

* Compute the difference between the function output and the previous output:

$$ \delta s = f(q_k) - s_{k-1} $$

* Compute the difference between the current input and the previous input:

$$ \delta q = q_k - q_{k-1} $$

* Update the Jacobian matrix approximation by minimizing the sum of squared residuals:

$$ \hat{J}_{f_k} = \arg\min \| \delta s - \hat{J}_f \delta q \|^2 $$

* Update the inputs:

$$ q_{k+1} = q_k - (\hat{J}_{f_k}^T \hat{J}_{f_k})^{-1} \hat{J}_{f_k} \delta s $$

* Update the function outputs:

$$ s_{k+1} = f(q_{k+1}) $$

* Check for convergence.
4. Convergence: Once the change in inputs and outputs is below a specified tolerance, the final approximation of the Jacobian matrix $J_f, \hat{J}_f$, is obtained.

The Least-Squares Minimization method is generally more numerically stable than the Broyden rule, especially when dealing with ill-conditioned matrices. However, it can be more computationally intensive. Moreover, this method can be sensitive to outliers in the data, which can significantly impact the estimated Jacobian.

### 2.1.3 Conclusion

This section offered a brief overview of the fundamentals of IBVS, which is used to control objects interacting with rigid scene, such as tracking a moving target [Papanikolopoulos et al., 1993]. Furthermore, we introduced two potential approaches for estimating the IBVS Jacobian matrix. These numerical estimations are useful when components of this matrix cannot be obtained directly from analytical form. The Broyden update rule has found application in the field of robotics. For example, [Hosoda & Asada, 1994] applied this approach to estimate a Jacobian matrix, allowing for the convergence of features to desired values without the need for prior knowledge of the kinematic structure of the robot system. [Jägersand & Nelson, 1996] utilized this technique within the context of visually-guided grasping. [Shahamiri & Jagersand, 2005] proposed an IBVS method, an IBVS approach without the knowledge of camera or robot calibration, to achieve positioning tasks while avoiding obstacles. Regarding the Least-Squares Minimization algorithm, [Shademan et al., 2010] employed a robust local least-squares (LLS) technique for the control of a robot in a IBVS context and for Jacobian estimation through Iteratively Reweighted Least Squares. The LLS method introduced by [Massoud Farahmand et al., 2007] enables the estimation of the Jacobian for any point from visual information within a nearby region of the point in question.

While these approaches were initially designed for robots interacting with rigid objects, which is not the primary focus of this thesis, some model-free methods have successfully adapted classical visual servoing principles for deformable object shape servoing. These methods will be further detailed in Section 2.5.

The next section presents principles of continuum mechanics and various approaches used to discretize object models, which are necessary for model-based approaches to manufacture deformable objects.
## 2.2 Deformation modeling

Continuum mechanics offers essential insights into the modeling and understanding of deformable objects. Deformable objects encompass a wide range of materials, from common substances like metals and polymers to biological tissues and geological formations, all of which can undergo deformation under various external forces.

Continuum mechanics provides a comprehensive framework for modeling the behavior of deformable materials in a continuous manner, treating materials as continuous media and encompassing concepts such as stress, strain, elasticity, and the partial differential equations governing their behavior. Therefore, this approach allows for highly accurate modeling of soft objects. To enable practical simulations of soft objects, numerical approximations based on continuum mechanics, such as the Finite Element Method (FEM) and mass-spring models (MSMs), have been developed. These methods facilitate efficient numerical solutions while introducing some level of approximation. The choice of method depends on the specific problem, the required level of accuracy, available computational resources, and whether real-time simulations are necessary.

In the subsequent sections, we will begin in Section 2.2.1 by providing a brief introduction to continuum mechanics, followed in Section 2.2.2, Section 2.2.3 and Section 2.2.4 by an exploration of various discretization methods used to approximate continuum mechanics. These methods play a crucial role in bridging the gap between the continuous mathematical descriptions of deformable objects and their practical computational implementation. Moreover, they showcase diverse levels of physical accuracy and simplicity, which are crucial considerations for the implementation of basic scenarios and adaptation to handle increasingly complex situations.

### 2.2.1 Overview on continuum mechanics

Continuum mechanics provides a framework for modeling deformable objects. Deformation refers to the alteration in shape or size of the object when subjected to external forces or loads. Deformation can be categorized as either elastic (reversible) or plastic (irreversible), depending on whether the object returns to its original shape after the forces are removed.

At time $t_0$, we consider the object in its rest state, not subject to any deformation. Let consider a point $P_i$ belonging to the object material. Its corresponding coordinates
at each time (t), $\mathbf{x}_{i(t)}$, and velocities $\dot{\mathbf{x}}_{i(t)}$, are given by:

$$ \mathbf{x}_{i(t)} = (x_{i(t)}, y_{i(t)}, z_{i(t)}) $$

$$ \dot{\mathbf{x}}_{i(t)} = (\dot{x}_{i(t)}, \dot{y}_{i(t)}, \dot{z}_{i(t)}) $$

Then, when external forces are exerted on the object, it will deform in response to these forces, causing its points to move to new positions corresponding to the deformed configuration. Let $\mathbf{u}_i$ denote the displacement of $P_i$, given by:

$$ \mathbf{u}_{i(t)} = (u_{ix}, u_{iy}, u_{iz}) = \mathbf{x}_{i(t)} - \mathbf{x}_{i(t_0)} \qquad (2.5) $$

An example of this deformation is represented in Figure 2.1.

<figure><img src="image_9.png" /><figcaption>Figure 2.1: Illustration of a deformation with the displacement vector $\mathbf{u}_i(t)$: (a) and (b) correspond to the undeformed object and deformed object, respectively.</figcaption></figure>

The object response to external displacements depends on its geometrical and physical properties. The deformation resulting from this displacement can be calculated by considering the stress tensor ($\sigma$), which represents the force applied per unit area of the object shape, and the strain ($\epsilon$), which quantifies the ratio of deformation to the original size of the object shape. Furthermore, different types of objects exhibit different responses. For instance, linear elastic objects undergo linear deformation, whereas viscoelastic objects exhibit partial resistance to deformations.

For example, to model linear elastic behaviors, the Cauchy strain tensor $\epsilon_C$ given in (2.6) is employed, but its validity is limited to small deformations. Nonetheless, it offers the advantage of fast computation, making it suitable for interactive applications. On the
other hand, the Green strain tensor $\epsilon_G$, given in (2.7), enables the modeling of nonlinear behaviors but comes with a higher computational cost. Consequently, it is less prevalent in real-time implementations.

$$ \epsilon_G = \frac{1}{2}(\nabla u + \nabla u^T) \in \mathbb{R}^{3 \times 3} \qquad (2.6) $$

$$ \epsilon_G = \frac{1}{2}(\nabla u + \nabla u^T + \nabla u^T \nabla u) \in \mathbb{R}^{3 \times 3} \qquad (2.7) $$

with

$$ \nabla u = \begin{bmatrix} \partial u_x / \partial x & \partial u_x / \partial y & \partial u_x / \partial z \\ \partial u_y / \partial x & \partial u_y / \partial y & \partial u_y / \partial z \\ \partial u_z / \partial x & \partial u_z / \partial y & \partial u_z / \partial z \end{bmatrix} \qquad (2.8) $$

Regarding the stress tensor, a common approach for linear objects is to apply Hooke's law, which states that stress is linearly proportional to strain. Mathematically, the relation can be expressed as follows:

$$ \sigma = E\epsilon \qquad (2.9) $$

with $E$ a tensor that depends on the material properties of the object, such as Young's modulus $E$ and Poisson's ratio $\nu$. Young's Modulus, also referred to as the modulus of elasticity, serves as a measure of a material resistance to deformation when subjected to an external force. It quantifies the relationship between stress and strain within the material, particularly when it experiences an axial deformation. Regarding Poisson's ratio, it expresses the ratio of transverse (lateral) strain to axial (longitudinal) strain when a material encounters uniaxial stress along one direction. In practical terms, Poisson's ratio characterizes how a material changes in dimensions, specifically in width or thickness, when subjected to stretching or compression. This dimensionless value is typically within the range of -1 to 0.5.

The dynamic behavior of the deformation using Newton's law is described as:

$$ M \ddot{u} = f(\dot{u}, u). \qquad (2.10) $$

where:

* $M$ is the mass matrix of the object;

* $(\dot{u}, \ddot{u})$ are respectively the first and second derivative of $u$.
* *f* represents both the internal and external forces acting on the object, often formulated as a set of partial differential equations.

By integrating the previous equation, the positions of the model points can be updated and used for simulating deformable objects. Some commonly integration methods include semi-implicit integration, implicit integration, Verlet integration, and Runge-Kutta integration. For a comprehensive review of numerical methods used, readers can refer to the work presented in [Hauth et al., 2003].

### 2.2.1.1 Conclusion

In conclusion, continuum mechanics serves as a fundamental framework for modeling the behavior of deformable materials, offering valuable insights into soft object deformation, stress, and strain. However, its mathematical complexity often requires the aid of numerical estimation methods to solve real-world problems efficiently. In the following sections, we explore various numerical methods and their crucial role in solving complex problems in continuum mechanics.

## 2.2.2 Particle-Based Models

Particle-based models are known for their computational efficiency and intuitive nature, however, they lack achieving high levels of physical accuracy. In particle-based methods, the material is represented as a collection of individual particles, each with properties like mass, position, velocity, and possibly other material-specific attributes.

### 2.2.2.1 Particle Systems

The particles, referred to as model points, initially exist in an equilibrium state that corresponds to their initial positions. When an external force is applied, the object undergoes deformation, resulting in the particles moving to new coordinates. The movement of the particles adheres to Newton's second law (2.10) and is governed by a time integration scheme with considerations for lifespan and collision detection.

Particles are commonly employed for modeling objects like clouds or liquids. However, there are also particle frameworks that utilize dynamically coupled particles to simulate solids [Tonnesen, 1998].

The advantage of particle systems lies in their simplicity, enabling the simulation of complex scenes with a large number of particles. Nonetheless, a limitation of particle
systems is the lack of an explicit surface definition. This drawback can pose challenges in applications such as tracking the restoration of elastic objects to their original form after a deformation was applied by a robot during a manipulation task.

### 2.2.2.2 Mass-Spring Systems

Another possible particle-based model is the mass-spring model (MSM), which approximates the behavior of a deformable object by the dynamics of a set of masses (points), linked by mass-less springs with dampers as connections to create a volumetric mesh.

Similar to particle systems, the movement of particles is simulated by employing Newton second law of motion. Moreover, a mass point $P_i$ is linked to finite neighbors $P_j$ with springs of rest length $l_{ij}^0$ and stiffness $K_{ij} = \text{diag}(k_x, k_y, k_z)$, where $(k_x, k_y, k_z)$ represent respectively the spring stiffness along the three axes $(x, y, z)$. Then (2.10) can be expressed for each point as follows:

$$m_i \ddot{x}_i = f_{si} + f_{Di} + f_{ei} \qquad (2.11)$$

where:

* $m_i$ is the mass of any point $P_i$ in a model with $N$ particles, where $i \in \mathcal{N} = [1, ..., N]$;

* $(x_i, \dot{x}_i = v_i, \ddot{x}_i = a_i)$ are respectively the 3D coordinates of $P_i$, its velocity and its acceleration;

* $f_{si}$ is the 3D force vector acting on $P_i$ due to springs with stiffness $K_{ij}$ connecting $P_i$ to its neighbors, $P_j, \forall j \in \nu_i \subset \mathcal{N}$, given by:

$$\begin{aligned} f_{si} &= \sum_{j \in \nu_i} f_{sij} = \sum_{j \in \nu_i} K_{ij} (\|x_i - x_j\| - l_{ij}^0) \frac{(x_j - x_i)}{\|x_j - x_i\|} \\ &= \sum_{j \in \nu_i} \alpha_{ij} K_{ij} (x_j - x_i) \end{aligned} \qquad (2.12)$$

* $f_{Di} = -D_v v_i$ is the 3D force vector acting on $P_i$ due to the damping $D_v$.

* $f_{ei}$ combines external forces, including gravity.

MSMs offer computational efficiency and intuitive representation, and are used for example for predicting and tracking object states during robotic manipulation, [Petit, Lippiello, Fontanelli, et al., 2017; Schulman, Lee, et al., 2013].
However, a challenge arises in tuning the spring stiffnesses and the damping values of the MSM to achieve the object desired deformation behavior in accordance with its material properties. To address this issue, some approaches proposed learning algorithms [Arriola-Rios & Wyatt, 2017; Bianchi et al., 2004; Morris & Salisbury, 2008]. On the other side, another approach derived an analytical expression for the stiffness of each tetrahedron constituting the mesh, which is proportional to the volume of the tetrahedron and the Young's modulus, [Lloyd et al., 2007]. Another limitation of MSMs is their inability to directly simulate volumetric effects, which prompted to introduce additional energy formulations to account for volume conservation [Teschner et al., 2004]. The orientation of springs also affects the behavior of a MSM, leading to introduce virtual springs to compensate for this effect [Bourguignon & Cani, 2000]. Furthermore, an innovative approach proposed incorporating extra elastic forces into the traditional MSM to accommodate complex mechanical behaviors like viscoelasticity, non-linearity, and incompressibility, [Xu et al., 2018].

### 2.2.2.3 Conclusion

In conclusion, particle systems offer a valuable tool for approximating continuum mechanics, providing a delicate balance between accuracy and computation time [Etzmuss et al., 2003]. While this section was limited on presenting particle systems and the MSM, it is important to note that various alternative particle-based approaches exist. For instance, there exists a method that uses a recurrent neural network to control MSM dynamics, effectively representing positions, velocities, and accelerations of mass points and springs within a mesh structure [Nurnberger et al., 1998]. Another alternative involves employing a convolutional neural network (CNN) to simulate mechanical loads, opting for the Poisson equation instead of a MSM [J. Zhang et al., 2019].

Next, we present Position-based dynamics, which is an approximation method in continuum mechanics that simulates physical behavior by iteratively enforcing constraints on the positions of discrete elements to achieve dynamic simulations of deformable materials.

### 2.2.3 Position-Based Dynamics

Unlike particle systems and MSMs, which rely on force-based approaches where model point velocities and positions are updated using a time integration scheme and given forces, position-based dynamics (PBD) models take a different approach. For a detailed
explanation of PBD, readers can refer to [Bender et al., 2017].

In PBD, the positions of particles are updated iteratively by satisfying a set of geometric constraints. These constraints can represent various physical properties and behaviors such as distance preservation, collision avoidance, bending, stretching, or even custom constraints specific to the simulation. Let consider a generic constraint that can be applied to any type of geometric constraint. The constraint equation can be written as:

$$C(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N) \approx 0 \qquad (2.13)$$

where C represents the constraint function. To compute the updated positions that satisfy the constraint, an iterative approach known as Gauss-Seidel projection can be used. It adjusts the positions of particles by calculating correction vectors based on the constraint function. At each iteration, the correction vector, $\delta\mathbf{x}_i$, for each particle i involved in the constraint is computed as:

$$\delta\mathbf{x}_i = -\lambda \frac{\nabla C_i}{\nabla C_i^T \nabla C_i} C_i(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N), \qquad (2.14)$$

where $\lambda$ is a stiffness parameter that controls the convergence rate, $\nabla C_i$ represents the gradient of the constraint function with respect to the position of particle i, and $C_i(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N)$ represents the value of the constraint function for particle i.

To update the positions of the particles, the correction vector is applied as follows:

$$\mathbf{x}_i = \mathbf{x}_i + \delta\mathbf{x}_i \qquad (2.15)$$

By iteratively applying these constraint equations and updating the positions of particles, PBD ensures that the geometric constraints are satisfied while simulating the desired physical behavior.

PBD methods offer several notable advantages. Firstly, they are characterized by their simplicity and computational efficiency. PBD methods do not require a mesh model, resulting in improved memory usage and computational speed. Moreover, their particle-based parallel nature allows for scalability, enabling simulations with a large number of particles. PBD has found extensive application across various interactive graphical applications [Macklin et al., 2014; Tian et al., 2013]. Its utilization has been particularly prominent in modeling the deformation of human body parts [Romeo et al., 2020; B. Zhu et al., 2008] as well as in robotic manipulation tasks [Caccamo et al., 2016; Güler et al.,
2017].

However, PBD has certain limitations. One limitation is that it does not explicitly represent the surface of the deformable object, which can be a drawback when it comes to visualizing or interacting with the object. Additionally, PBD heavily relies on geometric constraints and may require additional techniques or models to handle complex material behavior. On the other hand, MSMs offer greater control over material properties by adjusting the parameters of the springs. By using different types of springs, such as linear or nonlinear, a wide range of materials with varying stiffness, elasticity, and damping properties can be simulated, which allows for more fine-grained control over the behavior of the deformable object.

In the next section, we will explore further numerical techniques, distinct from particle systems and PBD, which are mesh-based techniques and are employed to approximate continuum mechanics in simulations.

### 2.2.4 Mesh-based Models

In contrast to the particle systems discussed in previous sections, this section delves into mesh-based models as an alternative approach to approximate continuum mechanics. Mesh-based models provide mathematical frameworks for characterizing material behavior, offering a distinct perspective on the analysis of deformable materials and structures.

#### 2.2.4.1 Finite Element Method

The finite element method (FEM) is an effective technique used to approximate the complex physical behavior exhibited by deformable objects. By dividing the deformable body into smaller and simpler finite elements, connected through nodes forming a mesh, FEM allows for accurate representation of the object in response to external forces. Unlike particle-based approaches, FEM operates on node displacements rather than particles. The deformation of a mesh containing $N$ points is determined by computing the displacement vector field $u$, which is obtained through solving (2.16):

$$Ma + Dv + Ku = f_e \qquad (2.16)$$

where $M \in \mathbb{R}^{3Nx3N}$, $D \in \mathbb{R}^{3Nx3N}$, $K \in \mathbb{R}^{3Nx3N}$ and $f_e \in \mathbb{R}^{3N}$ are respectively the mass matrix, damping matrix, stiffness matrix and external forces vector in a 3D dimension.
At every time step *t*, the computation of the matrix **K** relies on material characteristics like Young's modulus **E**, Poisson's ratio *ν*, and nodal displacement **u**.

The FEM has seen widespread use in robotics due to its capacity to generate realistic simulations and accurately represent complex deformations. It has been successfully applied in soft object tracking [Petit et al., 2018; Sengupta et al., 2019] and planning manipulation [Frank et al., 2014]. Nonetheless, a limitation of the FEM arises from the necessity to recalculate the dense matrix **K** at each time step, resulting in significant computational overhead. To address this concern, researchers have adopted an alternative approach wherein the matrix **K** remains constant. However, this method restricts the model capability to simulate solely small deformations. Alternatively, approaches such as co-rotational FEM [Müller & Gross, 2004] can be employed to reduce computational complexity while maintaining accuracy.

### 2.2.4.2 Boundary Element Method

The deformation of a surface utilizing the boundary element method (BEM) involves solving the equation (2.16) over the surface, as opposed to the volume-based approach employed by the FEM [Y.-M. Tang et al., 2006]. The surface is divided into N non-overlapping mesh elements, with node coordinates representing the centroids of these elements. In contrast to FEM, BEM offers notable speed improvements by utilizing a reduced number of nodes and elements. However, it is worth noting that the BEM is applicable only to objects with a composition of homogeneous materials inside.

Notably, BEM has found successful applications in real-time volumetric model simulation [James & Pai, 1999], and showed enhanced tracking accuracy in scenarios involving occlusions [Greminger & Nelson, 2008].

### 2.2.4.3 Modal Analysis

Modal analysis can be utilized to improve the computational efficiency of the FEM by reducing the number of degrees of freedom needed for analysis. The computational burden associated with the FEM arises from the calculation of motion using large matrices **M**, **K**, and **D** in (2.16). By incorporating modal analysis into the FEM workflow, the computational burden can be alleviated while still capturing the essential dynamic characteristics of the system [Pentland & Williams, 1989]. Here is a breakdown of the mathematical approach for enhancing FEM using modal analysis:
1. Discretization of the system: The continuous system is divided into a finite number of elements and nodes, where each element is defined by its geometry, material properties, and connectivity with neighboring elements.

2. Formulation of the equations of motion: They are derived following the same approach as in the FEM and shown in (2.16).

3. Modal analysis computation: Solve the eigenvalue problem $K\mathbf{v} = \lambda M\mathbf{v}$ to obtain the eigenvalues ($\lambda$) and corresponding eigenvectors ($\mathbf{v}$) of the system. The eigenvectors represent the mode shapes of the system, while the eigenvalues determine the natural frequencies and damping ratios.

4. Selection of reduced modes: Choose a subset of modes that significantly contribute to the system response based on criteria like mode participation factors, energy contributions, or engineering requirements. For example, by examining the eigenvalues calculated before, it becomes possible to remove the higher-frequency modes, which allows for the selective update of the most influential modes.

5. Projection of equations onto reduced mode set: Project the equations of motion onto the reduced mode set using the selected mode shapes.

6. Solution of the reduced system: Solve the reduced system of equations considering only the selected modes. This considerably reduces the computational workload as the number of degrees of freedom is reduced to the number of selected modes.

7. Computation of full response: Reconstruct the full response of the system by combining the contributions from the selected modes.

In modal analysis, assuming a constant value for $K$ can enhance computational efficiency. However, this assumption holds true only for simulating small linear deformations and introduces errors when dynamically simulating large deformations. Furthermore, modal analysis offers a robust tool for effectively simulating the characteristics of deformable objects [Basdogan, 2001; Hauser et al., 2003].

### 2.2.4.4 As-Rigid-As-Possible

The As-Rigid-As-Possible (ARAP) modeling approach is utilized to approximate the behavior of deformable objects by assuming a high level of rigidity. The main objective
of ARAP modeling is to preserve the object shape and rigidity while allowing for local
deformations [Sorkine & Alexa, 2007]. This is achieved by formulating an energy function
that enforces the rigidity constraint and encourages deformations that closely resemble
rigid motion.

The energy function typically consists of a rigidity term and a regularization term. The rigidity term penalizes non-rigid deformations by minimizing the differences in rotations and translations between neighboring vertices. The regularization term promotes smoothness and regularity in the deformation. It penalizes high curvature and encourages a smooth deformation. By minimizing the combined energy function, the ARAP modeling achieves a balance between preserving shape and rigidity and allows for localized deformations, providing a valuable tool for modeling deformable objects.

ARAP has gained significant popularity across diverse applications, such as 3D mesh
animation [Levi & Gotsman, 2014], regularization prior for tracking nonrigid scenes [Dou
et al., 2017], deformable object tracking during robotic manipulation [Han et al., 2018],
and shape servoing of soft objects [Shetab-Bushehri et al., 2022].

### 2.2.4.5 Conclusion

In summary, the selection of an appropriate modeling approach within the domain of
mesh-based models necessitates a careful consideration of the trade-off between precision
and computational efficiency. While FEM offers remarkable precision at the cost of in-
creased computational demands, Modal Analysis strikes a balance between accuracy and
efficiency, particularly for linear dynamic analysis. On the other hand, the ARAP method
excels in maintaining rigidity while prioritizing shape preservation, which is also capable
of simulating deformation with less accuracy than FEM. The ultimate choice should be
made in accordance with the specific application, the computational resources available,
and the requisite level of accuracy for the given task.

## 2.2.5 Conclusion

In this section, we presented the basics of continuum mechanics and different approaches
for approximating continuum mechanics. The choice between these modeling approaches
for deformable objects depends on the specific requirements of the application.

Particle systems such as MSM offer simplicity and real-time interactivity, making them
a valuable choice for applications where responsiveness and ease of implementation are
crucial. By representing objects as interconnected masses and springs, MSM captures the global deformations and overall structural characteristics of deformable materials efficiently.

On the other hand, PBD presents an intermediate representation between particle systems and mesh-based models, offering a simplified yet efficient approach to enforcing constraints, handling collisions, and simulating deformations in real-time. PBD is well-suited for applications that require a balance between responsiveness and accuracy.

Finally, for high-fidelity simulations with complex material behavior, mesh-based models may be more appropriate but may come at a higher computational cost. For example, by dividing objects into smaller elements, the FEM offers a detailed and accurate representation of deformable objects, making it essential in applications where realism and fidelity are paramount.

In the next section, we will present different approaches for tracking deformable objects and discuss parameter estimation for the physics-based models.

## 2.3 Deformable object tracking

The initial shape of the object is represented using a topological model before any deformation is applied to it. Then, deformable object tracking is employed to adapt the object shape online during its deformation by continuously tracking it. The process of deformable object tracking can be approached using either a physics-based model or a model-free method. Incorporating an explicit physics-based model can significantly enhance the tracking process by enabling reliable computation and prediction of internal forces experienced by the object. These forces are often not directly measurable by visual or force sensors. Conversely, model-free approaches are better suited for reconstructing soft object models rather than tracking them. Since topological model reconstruction falls outside the scope of this thesis, we will focus hereafter mainly on the most relevant methods for model-free and model-based visual tracking of deformable objects.

### 2.3.1 Model-free

Model-Free approaches do not require a predefined model of the object being tracked. Instead, they directly track the object within the camera frames without any prior knowledge about object shape or deformation properties.
2.3.1.1 Linear deformable object tracking by CPD

A framework that encompasses state estimation, task planning, and trajectory planning, all relying on the coherent point drift (CPD) [Myronenko & Song, 2010], was introduced in [T. Tang et al., 2018]. CPD is a registration method that facilitates the non-rigid mapping of one point set to another (see Figure 2.2). We particularly focus on the state estimation proposed in their work, which is applied to a deformable linear object (DLO) using stereo cameras and CPD. A real-time observer was proposed to estimate the states of deformable object by using its corresponding point cloud. Considering that the DLO is often occluded by robot arms or self-occluded during manipulation, the state estimator is designed with robust occlusion handling capabilities. The state estimation process involves acquiring the position of each node on the object by registering the previous estimation results with the new point cloud measurement. As a result, the object state can be robustly estimated in real-time, effectively handling noise, outliers, and occlusions.

The experimental results presented show that the authors approach can successfully
and efficiently track a DLO (corresponding in practice to a rope) in real-time, even when
dealing with noisy and occasionally occluded point clouds caused by obstacles.

<figure><img src="image_10.png" /><figcaption>Figure 2.2: Example of CPD registration, from [T. Tang et al., 2018].</figcaption></figure>

2.3.1.2 Occlusion-robust deformable object tracking by regularizing the CPD output

Another model-free occlusion-robust RGB-D sequence tracking framework also based on
CPD was proposed in [Chi & Berenson, 2019]. In contrast to the previous approach, the
authors stated that CPD only enforces motion coherence between consecutive RGB-D
frames, and therefore, there is a possibility of accumulating errors between the current
tracking state shape and the true state. As a result, the topology of the tracking result may deviate from the original model, even though the tracking result remains statistically correct when viewed as a point set registration problem. To address this challenge, the authors propose the use of a regularization term based on locally linear embedding (LLE), wherein LLE performs non-linear dimensional reduction while preserving the local neighborhood structure. This regularization term ensures topological consistency with respect to the original model.

Several experiments were presented, involving the tracking of rope and 2D cloth, to demonstrate the improved robustness of the approach against occlusion. These experiments included both simulation and real-world data.

### 2.3.1.3 Conclusion

The presented model-free approaches enable the state estimation of 1D and 2D deformable objects by aligning the last step state estimation with the current point cloud measurement, without requiring a physics-based model. Furthermore, these approaches have demonstrated their effectiveness and robustness against object occlusion.

## 2.3.2 Model-based

Model-based approaches rely on predefined models to describe the object topology, its physics-behavior and its initial undeformed state. Then, they use registration techniques to estimate the update state of the object and optimization techniques to estimate parameters that define an object deformation.

### 2.3.2.1 Deformable object tracking for model registration and parameter tuning

An approach focusing on tracking and modeling deformation in soft objects using a data-driven approach was presented by [Wang et al., 2015]. Their method involves an iterative framework with two main components: physics-based deformation tracking and an optimization of deformation parameters. The soft objects are modeled using the co-rotational linear Finite Element Method. The optimized deformation model enhances the accuracy of tracking results, leading to improved deformation parameter estimation in subsequent iterations. Notably, their proposed optimization considers not only material elasticity parameters and damping coefficients but also incorporates a reference shape. It refers to the
object geometry subject to zero external forces and can be different from the static shape,
which corresponds to the equilibrium geometry under gravity. For a visual representation
of the complete pipeline, refer to Figure 2.3.

The process begins with scanning and reconstructing the static shape (Figure 2.3(a)). To capture the dynamic motion, three depth cameras from different viewpoints are utilized, resulting in a time sequence of point clouds (Figure 2.3(b)). Employing a physics-based probabilistic approach and using the static shape as a template, the dynamic point cloud sequence is tracked (Figure 2.3(c)). The result is represented as 3D coordinates of a tetrahedral mesh, closely enclosing the surface template of the static shape. Subsequently, the parameter estimation component optimizes the elastic material parameters, damping coefficients, and reference shape (Figure 2.3(d-e)).

Experimental results validated the ability to track soft objects while simultaneously
recovering the reference shape and estimating the deformation parameters. However, this
approach necessitates a high accurate initial mesh of the object.

<figure><img src="image_11.png" /><figcaption>Figure 2.3: The framework of a soft object tracking and estimation of its model parame-
ters, from [Wang et al., 2015].</figcaption></figure>

2.3.2.2 Real-time 3D deformable object tracking using RGB-D camera

Another model-based tracking approach using RGB-D camera, which involves a real-time
tracking method designed for 3D soft objects that experiencing significant deformations,
including elastic deformations and rigid motions was proposed by [Petit et al., 2015].
The approach considers point cloud data obtained from an RGB-D sensor and utilizes
the Finite Element Method for modeling the soft object. It incorporates a consistent
mesh representation and takes into account known material properties such as the Young
modulus and Poisson ratio. The process involves initially registering the segmented point
cloud in a rigid manner, followed by non-rigid fitting of the mesh. Geometrical point-to-
point correspondences are used to calculate external forces to be applied on the mesh.
The entire workflow is illustrated in Figure 2.4.
The object displacement can be decomposed into a rigid transformation and a pure
deformation, with rigid translation and rotation transformations being associated with
the rigid transformation. To achieve this, a rigid Iterative Closest Point (ICP) process
is applied to estimate a rigid transformation from the resulting segmented point cloud
to the mesh. Subsequently, the latter is rigidly aligned to the point cloud data. Then,
a deformable registration problem is addressed to fit this point cloud data to the rigidly
aligned tetrahedral mesh. This deformable registration involves determining external
forces exerted by the point cloud on the mesh and integrating these forces, along with the
internal forces computed using the physical model. The computation of external forces
relies on geometrical point-to-point nearest neighbor correspondences by employing K-d
tree searches from the segmented point cloud to the visible surface of the rigidly aligned
mesh and vice versa. Processing both sets of correspondences is crucial because relying
solely on a single match from one set to another might lead to inconsistent matches.

Simulation and experimental results confirm that the approach proposed by the au-
thors shows promise as a real-time tracking method capable of effectively managing diverse
deformations and motions by combining the FEM, an efficient segmentation method, and
conventional point cloud registration techniques.

An extension of the previous work was developed in [Petit, Lippiello, Fontanelli, et al., 2017]. This study involves testing various physical models, such as MSM, FEM, and co-rotational FEM, and comparing the tracking accuracy between these models. Additionally, the researchers applied their method to track the deformation of a pizza in a challenging scenario, where the pizza was manipulated by a humanoid pizza chef robot through stretching and tossing. Another extension was proposed in [Petit et al., 2018], which enables the registration of multiple objects within a single scene. This was not possible before, and it is achieved by effectively handling collision detection model between the different objects of the scene, which are consistently simulated using the co-rotational FEM.

2.3.2.3 Real-time 3D soft tissue tracking in 3D ultrasound images

Another model-based tracking approach to achieve real-time tracking of deformable struc-
tures (targets) in 3D ultrasound sequences was introduced by [Royer et al., 2017]. This
method involves obtaining the target rigid and non-rigid displacements by minimizing a
visual error criteria relying on dense ultrasound visual observation and physics-based sim-
ulation. The approach integrates a physics-based model based on a MSM and a 3D object
<figure><img src="image_12.png" /><figcaption>Figure 2.4: The framework of a soft object tracking based on FEM, from [Petit et al., 2015].</figcaption></figure>
mesh to estimate the target motions across consecutive 3D images. For this purpose, the
mesh vertices displacements are online updated by iteratively summing their displace-
ments resulted from the MSM internal forces with external displacements that minimize
a visual error based on dense ultrasound feedback. This method also includes an ultra-
sound shadow detection process that increases the tracking robustness by a confidence
map used as a weighting mechanism in the visual error criterion. The computational flow
of the method is summarized in Figure 2.5.

The primary goal of this approach is to iteratively estimate both the external and inter-
nal displacements of the mesh. For computing the external displacements, they adopted
an intensity-based method, where the cost function is assessed using various dissimilarity
functions such as Sum of Squared Differences, Weighted Sum of Squared Differences, Sum
of Conditional Variance, and the newly proposed Sum of Confident Conditional Variance.
Concerning the estimation of internal displacements, it is accomplished using the MSM.

The results illustrate that this approach, utilizing the dissimilarity measure Sum of Confident Conditional Variance, provides accurate motion estimation and performs well on ultrasound sequences affected by various challenges, such as speckle noise, large shadows, and ultrasound gain variation. However, it should be noted that the proposed method may encounter inaccuracies when the target experiences significant deformation beyond the limits defined by the provided elastic parameters.

<figure><img src="image_13.png" /><figcaption>Figure 2.5: The computational flow of the tracking method based on a MSM, from [Royer et al., 2017].</figcaption></figure>
### 2.3.2.4 Deformable object tracking using coarse physics-based model

Later, [Sengupta et al., 2019] employed a model-based method for tracking deformable objects using a RGB-D camera. They utilized a coarse co-rotational Finite Element Method to represent the physical behavior of the object, and minimize a point-to-plane distance-based geometric error between the point cloud and the mesh by considering virtual forces. The proposed approach does not require an exact physical model of the object or precise properties like Young modulus and Poisson ratio for accurate tracking.

The registration process starts with the same initial step as previous approaches [Petit et al., 2018; Petit, Lippiello, Fontanelli, et al., 2017; Petit et al., 2015], which is rigid registration. For that, two features, namely depth-based geometric error and keypoint-based feature tracking [Trinh et al., 2018], take place. The second step is non-rigid registration. In this part, only a subset of points of interest on the object surface, that are referred to as clusters, are considered to determine the location where a limited number of external forces needed to be applied for deforming the entire mesh. The positions of these clusters correspond to the parts of the mesh that deviate most from the current point cloud. The forces applied on these clusters are then computed to minimize the geometric point-to-plane distance between the RGB-D point cloud and the entire model surface mesh. This is accomplished by establishing a numerical Jacobian that correlates the variation in geometric error with the variation of force applied. Indeed, an Iteratively Reweighted Least Square formulation was employed to minimize the geometric error using this Jacobian.

Simulation and experimental results demonstrate the effectiveness of the authors proposed real-time tracking approach. In contrast to [Petit et al., 2018], the approach achieves accurate tracking of the deforming object even when coarse parameters are considered.

### 2.3.2.5 Conclusion

The first step of the model-based approach involves generating a 3D tetrahedral mesh model associated with the soft object. Once this model is defined, the approach integrates a physics-based model with the previously described mesh to estimate the target motions across consecutive 3D images. For this purpose, the vertices displacements are calculated by iteratively summing the internal displacements estimated from the mechanical component and the external displacements computed using different methods, such as point-to-point correspondences and an intensity-based approach combined with a shadow
detection process, among others.

### 2.3.3 Conclusion

In this section, we discussed various approaches to deformable object visual tracking, ranging from model-free methods to model-based techniques. Moreover, deformable object tracking offers diverse applications, one of which is reducing the disparity between real object deformation and its simulated deformation. Reducing this gap enhances the accuracy of the simulated object and allows for fine-tuning of model parameters, as we discuss in more details in the following section.

## 2.4 Parameter estimation methods

Careful selections of material models and their corresponding parameters is crucial to achieve lifelike soft deformations. The main concept behind the parameter estimation methods is to iteratively adjust the parameters of a simulation system until the deformations obtained in simulation closely match those measured on the real object. In this section, we focus on estimating the elasticity parameters of both MSM and FEM models.

A common preliminary step for all subsequent methods is to use dense 3D reconstruction for object surface mesh generation, comprising N points to represent the undeformed object. Afterward, the model parameters are estimated through an optimization process aimed at identifying the values of the model elasticity parameters $p$ by applying external forces to the object ($f_e$) and reducing the distance between the model points and a reference model:

$$ \hat{\mathbf{p}} = \arg \min_{\mathbf{p}} \sum_{i=1}^{N} \| \mathbf{x}_i(\mathbf{p}, \mathbf{f}_e) - \bar{\mathbf{x}}_i \|^2 \quad (2.17) $$

with $N = [1, 2, \dots, N]$, and $\mathbf{x}_i(\mathbf{p}, \mathbf{f}_e)$ represents the coordinates of model point $P_i$, which depend on both the model parameters, $p$, and the applied force, $f_e$. In contrast, $\bar{\mathbf{x}}_i$ corresponds to the measured point coordinates obtained from the reference object. This reference model that serves as the ground truth can be either an object captured using an RGB-D sensor or a pre-existing, accurately modeled object. In the latter case, $\mathbf{x}_i$ belongs to a coarse model.
## 2.4.1 Estimation of non-linear heterogeneous FEM parameters

A data-driven approach for approximating parameters of a non-linear heterogeneous soft tissue was presented by [Bickel et al., 2009]. The object is modeled using a FEM and in order to capture material non-linearity, the authors defined the parameter values of the mesh-based model through data interpolation in the strain-space. A gradual load $f_e$ is applied to the object using a probe. At each loading step, with knowledge of the applied force, a sample of the stress-strain relationship is estimated, and the material parameters are computed. By considering measured displacements and applied forces, spatially varying material parameters are computed using an updated version of (2.17) which is given by:

$$ \hat{\mathbf{p}} = \arg \min_{\mathbf{p}} \left\{ \sum_{i=1}^{N} \| \mathbf{x}_i(\mathbf{p}, \mathbf{f}) - \bar{\mathbf{x}}_i \|^2 + \gamma \| \mathbf{L} \mathbf{p} \|^2 \right\} \quad (2.18) $$

The spatial smoothness of parameters is ensured by employing the sparse Laplacian matrix $\mathbf{L}$, while the non-linear residual equation presented in (2.18) is iteratively solved using the Levenberg-Marquardt algorithm [Moré, 2006].

The experimental results validated the approach of estimating the object material elasticity using a data acquisition system that comprised force probes and a marker-based trinocular stereo system. Nevertheless, it is important to note that the approach has certain limitations. For instance, the minimization process may encounter local minima, potentially resulting in inaccurate deformation synthesis. Additionally, the method does not account for viscosity and plasticity parameters.

Another method for estimating the Young's modulus of heterogeneous soft tissues modeled by FEM was proposed by [Coevoet et al., 2015]. The primary innovation involved employing an interactive inverse real-time simulation combined with manual registration of a restricted set of points. This approach relies on utilizing a non-linear FEM within a constraint-based framework. By receiving a limited number of registered points from the user, the method performs real-time optimization to adjust FEM parameters and achieve appropriate geometric deformations. The goal is to estimate the Young modulus from an observed deformation, but for accurate estimation, the applied forces on the model need to be known. To address this, the FEM model equations are projected into the constraint space using Schur complement, effectively reducing the optimization space.

The validation procedure involves several steps. Firstly, an arbitrary deformation of a deformable object is created. Next, few relevant points from the model are selected,
and their positions are recorded when equilibrium is achieved. Then, the simulation is reset without the deformation, and the inverse simulation approach is used to register the chosen points. Finally, a comparison is made between the actual values of parameters used in the process and the values estimated through the inverse method. This comparison is utilized to update the parameter estimation.

Simulation results are utilized to validate the proposed approach for estimating the heterogeneous material, consisting of three distinct Young moduli, one per object part. Various forces were applied to these parts of the object, and the approach effectively estimated the three Young moduli with an error of less than 3%.

### 2.4.2 Estimation of homogeneous FEM elasticity parameters using a robotic manipulator

Another approach, which focuses on determining the elasticity parameters of a FEM, involves the use of a robotic manipulator equipped with a force sensor and was introduced by [Frank et al., 2010]. The authors estimated the object elasticity parameters, the Young's modulus and the Poisson's ratio, by establishing a relationship between applied forces and resulting surface deformations. The robot exerted a sequence of forces on the soft object while simultaneously observing both the surface of the model object and the real object using a depth camera. Next, an appropriate error function was defined (similar to (2.17)), reflecting the differences between the measured and simulated surfaces observed by the robot. To establish correspondences between the simulated and segmented object surfaces, a 3D registration technique based on point-clouds was utilized. Subsequently, they employed a numerically approximated gradient-descent-based error minimization technique, effectively reducing the disparities between the actual and simulated deformations with respect to the elasticity parameters.

Simulation and experiments confirmed the accuracy of estimating appropriate parameters for various soft objects. To reduce occlusions resulting from the manipulator deforming the object, a thin wooden stick was attached to the end-effector. It is important to note that this approach was restricted to a linear, isotropic, and homogeneous deformation model.
### 2.4.3 Estimation of MSM parameters using a FEM reference model

A different approach estimating the parameters of an isotropic MSM by leveraging known parameters of the same object modeled with FEM was presented by [Natsupakpong & Çavuşoğlu, 2010]. The authors relied on the premise that FEMs derive accurate material properties of the object, which are assumed to be known. On the other hand, MSMs do not directly approximate the physical behavior of the object. Consequently, the proposed approach aimed to approximate a FEM by optimizing the MSM parameters, while minimizing the matrix norm of the error between the MSM and linear FEM stiffness matrices of the same object.

The approach starts by creating a MSM with the same size and geometry as the FEM, where the elements of the MSM are fully connected. The nodal mass values of the elements in the MSM are assigned to match the nodal elements mass in the FEM, approximating its dense mass matrix by a diagonal one. Concerning the stiffness matrix of each element in the MSM, it is defined based on the unknown spring constants. On the other hand, the stiffness matrix for the FEM element is calculated using a linear elastic model and expressed in terms of the known FEM parameters. Due to the assumption of isotropic material, the MSM element only necessitates two parameters to achieve symmetry, whereas the FEM element possesses two independent elastic parameters. To identify the parameters of the MSM element, both models are subjected to identical displacements. Subsequently, the discrepancy in nodal forces between the two models is assessed using the previously mentioned stiffness matrices. To minimize this error (2.17), an optimization process is undertaken, employing an appropriate matrix norm. It is crucial to note that achieving a zero error is not feasible due to the inherent structural differences between the stiffness matrices of the FEM and MSM elements.

Simulation results demonstrate the efficiency of estimating MSM parameters by using a FEM as a reference model, and minimizing the error between their stiffness matrices through an optimization process. The validity of this approach was further confirmed by validating it across four different configurations. Planar object examples employed triangular and quadrilateral meshes, while three-dimensional volumetric object examples utilized tetrahedral and hexahedral meshes. Despite its effectiveness, the approach does have several limitations. The method is not suitable for arbitrarily connected MSMs, and it calculates the stiffness matrix per element. Furthermore, while there is an alternative
of conducting the optimization at the whole object level instead of per element, it is not feasible due to the high computational cost and the potential risk of encountering local minima problems.

### 2.4.4 Estimation of isotropic co-rotational FEM parameters using a robotic manipulator

Another approach to estimate the Young Modulus and Poisson ratio parameters within a co-rotational FEM was proposed by [Petit, Ficuciello, et al., 2017]. The authors assumed a continuous isotropic material. The object elastic parameters are determined by tracking its deformations using a RGB-D camera, and incorporating the force measurements obtained from a force sensor. This estimation process follows a similar approach to the one proposed in [Frank et al., 2010], involving the minimization of a fitting error between the real deformations obtained from the RGB-D camera and the simulated ones.

The parameter estimation problem is tackled by minimizing the deviation between the simulated deformations and the observed ones (same form as (2.17)). The innovative aspect of this approach lies in its consideration of distance as the summation of both the distance from the simulated model to the real segmented object and vice versa. Addressing the non-linear optimization problem with respect to the model parameters, the evaluation of the objective function becomes expensive, and computing its gradients becomes non-trivial, thereby making gradient-based optimization methods impractical. To overcome this, the authors employed the gradient-free Nelder-Mead method [Nelder & Mead, 1965].

Experimental tests were conducted to assess the efficacy of the parameters estimation technique. During the experimentation with different initial configurations, convergence of the model parameter estimation was achieved, albeit the presence of local minima was observed.

### 2.4.5 Estimation of coarse FEM parameters

Another approach to track the deformation of soft objects and estimate their elasticity parameters was presented by [Sengupta et al., 2020]. This method involves fusing visual data obtained from an RGB-D camera with interactive FEM simulations of the object. The study primarily focused on estimating the Young's modulus, assuming that the Poisson ratio was known beforehand. To estimate the elasticity parameters at a contact point of the deformable object, the proposed algorithm relies on the tracking results and external
measurements of the deformation forces obtained at the same point.

The elasticity parameter estimation process involves minimizing the error between a tracked object and a simulated object that is deformed by a robot equipped with a force sensor. External forces are applied at a specific point, and their measurements are stored. Once a steady-state in the deformation tracking is achieved, these measurements are transmitted to the estimation algorithm. A cost-function is computed based on the Euclidean distance error between the object points, where the forces were applied and obtained from the camera, and their corresponding model points at the end of the simulation, where the same forces are replicated. The cost-function is minimized using the Levenberg-Marquardt algorithm. The authors proposed calculating a Jacobian matrix of the cost-function to update the elasticity parameters. This matrix is numerically estimated by conducting multiple simulations with varying parameter values.

Experimental results provide evidence of the effectiveness of their method in estimating
the elasticity parameters using the applied forces and a coarse 3D mesh of the object.

**2.4.6 Conclusion**

In this section, we presented various approaches to estimate the elasticity parameters of deformable objects. These objects can be modeled using discrete models such as MSM or mesh-based models like FEM. The primary goal of parameter estimation is to minimize discrepancies between real object deformations and their simulated deformations when working with real data. Additionally, it serves to reduce disparities between an accurate model of the object and a simplified, potentially less computationally one. Furthermore, parameter estimation proves particularly advantageous in model-based deformable object manipulation, especially when dealing with objects of unknown parameters, as will be discussed in the following section.

**2.5 Deformable object manipulation**

Robotic manipulation of rigid objects has been studied and successfully implemented in
several traditional application fields. However, the manipulation of deformable objects is
still an open problem since the object response to an external mechanical action cannot
be predicted in the same simple manner as for rigid objects. Manipulating a cable,
tearing a 2D cloth, controlling the 3D shape of biological tissues are some examples of
challenging soft object manipulations. Achieving precise and controlled shape changes in these objects is important for a wide range of applications, such as surgical and dressing assistance [Shin et al., 2019; Yu et al., 2017], food handling [Lehnert et al., 2017], and the automotive industry [Shah & Shah, 2016], to name a few. The automation of these tasks combines many areas such as control, perception and modeling/learning, where several issues complicate it: i) the material properties of the object (stiffness, mass, and viscosity), ii) the object geometry (dimension and shape), iii) the forces to be applied (localization, direction and magnitude).

Ensuring stable and dexterous manipulation of soft objects without causing damage is an ongoing challenge. To address it, researchers have developed a variety of shape servoing techniques, which can be broadly categorized into two categories: model-based and model-free techniques. Additionally, the object shape can be represented by some points on its surface, which can be defined by their 3D positions or their pixel coordinates corresponding to the projection of their 3D coordinates in the camera frame. Alternatively, the shape can be defined by a combination of selected object surface points or by considering all surface points.

In the following sections, we begin by reviewing model-free approaches that focus on deformable object manipulation. We then proceed to discuss model-based approaches, which rely on a physics-based model corresponding to the deformable object. Model-based approaches incorporate this physics-based model, while model-free approaches rely solely on visual information.

### 2.5.1 Model-free

This section concerns the model-free approaches. These methods are mainly based on visual servoing and machine learning, and do not require prior information on deformation behavior deduced from physics-based models. In this thesis, we will specifically focus on the approaches based on visual servoing.

#### 2.5.1.1 Vision-based object 2D deformation features positioning

One of the pioneering works in the field of model-free approaches was introduced by [Navarro-Alarcon et al., 2013]. They presented a novel method for controlling robot manipulators interacting with unknown elastic objects based on visual feedback. Their approach involved estimating the deformation Jacobian matrix in real-time using Broy-
den update rule [Broyden, 2000]. This matrix, as the features Jacobian matrix in classical
visual servoing, captures the relationship between the visual feedback features and the
robot manipulators motions, relying solely on the visual information. Their closed-loop
dynamic-state feedback velocity control law, implemented within a passivity-based frame-
work, ensured input-to-state stability in the presence of imperfect estimations and external
disturbances.

The authors introduced four distinct types of features presented in Figure 2.6. These
include point-based deformation, which explicitly focuses on displacing singular/multiple
2D points. Distance-based deformation regulates the distances between point features.
Angle-based deformation focuses on the relative angle of a line of interest constructed from
two 2D points. Finally, curvature-based deformation involves measuring and controlling
the contour of interest defined from three 2D points. These feature vectors were defined
as nonlinear functions of the visual feedback points.

Afterwards, the change in these feature vectors was represented in relation to the
robot manipulators using a Jacobian matrix. Utilizing this Jacobian matrix, a control law
based on energy-based dynamic-state feedback velocity was derived. The objective of this
control law was to convert the deformation behavior into a non-conservative Hamiltonian
dynamical system, aiming to achieve the desired deformation by minimizing an energy-like
functional. The authors stated that even with imperfect estimations of the deformation
Jacobian matrix, the closed-loop system exhibited input-to-state stability and complete
dissipativity in the presence of external disturbances.

The experimental results demonstrate the effectiveness of their control method in
handling real-time deformation tasks without the need for camera calibration parameters
or modeling of the deformation. However, it should be noted that the proposed method
has certain limitations. It requires slow motion of the robot manipulators and low-pass
filtering of the observation signals. Additionally, the method is applicable only to quasi-
static and elastic objects.

Subsequently, two vision-based control methods for shaping soft objects, both of which
did not rely on the identification of the deformation model or camera parameters, was
introduced in [Navarro-Alarcon et al., 2014]. The first controller enables real-time esti-
mation of the deformation Jacobian matrix. The second approach introduces an adaptive
deformation controller and combines a minimal amount of offline information with online
adaptation of variable parameters. Prior to conducting the experiments, the method re-
quires gathering some test points, which involves non-collinear end-effector displacement
<figure><img src="image_14.png" /><figcaption>Figure 2.6: Different types of feature-based deformation, from [Navarro-Alarcon et al., 2013]: (a), (b), (c), and (d) correspond to point-based deformation, distance-based deformation, angle-based deformation, and curvature-based deformation, respectively.</figcaption></figure>

vectors and visual feedback vectors. These vectors are used to construct a regression matrix, which remains constant throughout the online deformation test. This regression matrix is subsequently used to create an estimated projection model for a feature point onto the image plane (see Figure 2.7). The parameter adaptation rule is determined by a combination of two terms. The first term represents the difference between the current position and desired position of the test points. The second term accounts for the error in the estimated projection model (see Figure 2.7).

The effectiveness of both controllers was evaluated by implementing them on a robot manipulator. However, there are still several limitations associated with these proposed controllers. The first approach is sensitive to image noise, but the authors addressed this issue by employing a low-pass filter on the signals. On the other hand, the second control method is less affected by noise, but it may not be suitable for tasks that require a broad range of nonlinear deformations. This method is therefore better suited to the applications of small or localized deformations to the object.
<figure><img src="image_15.png" /><figcaption>Figure 2.7: Geometric representation of the projection estimation error, from [Navarro-Alarcon et al., 2014].</figcaption></figure>

### 2.5.1.2 Simultaneous object positioning and shaping

In a subsequent work, a model-free method for automatic 3D shape servoing of soft objects using robotic manipulators was proposed by [Navarro-Alarcon et al., 2016]. The object shape is represented by a deformation feature vector consisting of two distinct components: position and shape components. This approach enables the manipulation of an object to reach a desired position while controlling its relative deformation. The pipeline of this method is depicted in Figure 2.8.

The authors stated that the deformation model can be parameterized linearly using an adaptive vector of unknown deformation parameters. Subsequently, they developed an algorithm that dynamically estimates this vector and iteratively approximates the deformation model, which does not necessarily correspond to the true deformation model. Afterwards, the estimated vector is used to approximate the deformation Jacobian matrix, which is necessary for computing the velocity control input with conventional controllers. Furthermore, the visual feedback consists of two components: a positioning feedback term and a shaping feedback term. The first term involves the location of explicit feature points or the geometric center of multiple feature points, and is used to control the object position. The second term evaluates the internal displacements of the object, allowing for control of relative deformations. It encompasses the compression distance, which measures the relative distance between two feature points. Additionally, a folding angle of the object is defined as the measured angle between two intersecting lines computed using feature points and a normalized curvature that quantifies the roundness of the object. All these shape geometric descriptors are then servoed by the controller to reach desired values.

The method was tested using two manipulators and the proposed controller to ma-
nipulate the shape of a soft object and drive it towards a desired shape. However, this
method does have certain limitations. For instance, the method is specifically designed
for slow movements of the robots and may not be effective in controlling rapidly changing
shapes. Another limitation is the assumption of the object model being purely elastic.

<figure><img src="image_16.png" /><figcaption>Figure 2.8: The pipeline of a 3D shape servoing proposed in [Navarro-Alarcon et al., 2016].</figcaption></figure>

### 2.5.1.3 2D contour manipulation

In a later study, a model-free approach for automatically deforming soft objects into desired two-dimensional shapes using robot manipulators was proposed by [Navarro-Alarcon & Liu, 2017]. The shape of the object is represented by its contour, and a compact feedback characterization of deformable object shapes is developed based on a truncated Fourier series. This method differs from previous approaches that typically rely on displacements of feature points to define the servoing task. Instead, it introduced closed-loop deformation tasks by utilizing the truncated Fourier series of 2D image contours as a form of visual feedback. The closed-loop controller is depicted in Figure 2.9.

The object contour can be reconstructed using the shape parameters derived from the Fourier approximation. Therefore, the shape error corresponds to the difference between the Fourier coefficients of the current shape representation and the one corresponding to the desired shape. Regarding the control law, it was designed using a Jacobian-based approach. To estimate the Jacobian matrix and to overcome the challenges associated with the high dimensionality and nonlinear behavior of the object model, the authors proposed an online algorithm to approximate the deformation properties. This approach avoids the requirement of full parametric identification of the object model, as seen in previous approaches. The proposed method employs a local deformation model that is recalibrated at different local regions. As the robot moves into a new configuration, the model is recalibrated using data points collected along the local trajectory. The objective
of the iterative model estimator is to compute a local approximation of the Jacobian
matrix. However, it is important to note that since the estimator computes a locally valid
linearized model, it cannot capture the entire properties of the object in general.

The approach was validated through both numerical simulations and experimental
results. It is important to note that the proposed method is specifically designed to control
the 2D image projections of object contours and solely involves translational movements.
Additionally, it is worth mentioning that the system run into difficulties if the desired
contour is unreachable or if the number of Fourier coefficients is insufficient to accurately
describe the object shape. In such cases, the system may converge to a local minimum
configuration.

<figure><img src="image_17.png" /><figcaption>Figure 2.9: The pipeline of a model-free shape servoing using Fourier series represented in [Navarro-Alarcon & Liu, 2017].</figcaption></figure>

Another approach, which focuses on controlling the shape of soft objects, was intro-
duced by [Qi et al., 2021]. In this approach, the object shape is also represented through
its contour extracted from a 2D image of the observed object. The object contour was
parametrized using a set of 10 parameters. These parameters include the object cen-
troid coordinates (2 parameters), the logarithmic function of the 7 Hu moments [M.-K.
Hu, 1962], and the orientation of the contour principal axis in the plane [Zheng et al.,
2019]. These parameters are subsequently normalized. Similar to previous methods, this
approach establishes a relationship between changes in the object shape parameters and
the robot motions through the use of a deformation Jacobian matrix. Subsequently, a
sliding mode controller was proposed to automatically reshape the object to a desired
configuration while simultaneously estimating the Jacobian matrix. To assess the effec-
tiveness of this approach, experiments were conducted to validate the proposed control
and reshape various elastic and plastic objects to achieve desired contours. It is worth
noting, however, that this method comes with certain limitations. It is not suitable for
purely plastic objects, and the object deformation process must occur slowly.
2.5.1.4 Fourier-based 1D shape control

An alternative method employing Fourier series for model-free shape servoing was pro-
posed by [J. Zhu et al., 2018]. The study focused on manipulating a 1D linear cable, where
the shape of the cable was approximated using Fourier series. The manipulation was car-
ried out using two robotic manipulators. The objective of this research was to establish a
framework for multi-arm robots to effectively deform a flexible cable into a desired shape
through vision-based control. The local deformation model was estimated online using
shape parameters, allowing for real-time adaptation. The local deformation model is sim-
ilar to the one proposed in [Navarro-Alarcon & Liu, 2017], with the distinction being that
this method incorporates multiple objects instead of just one, along with the inclusion of
rotational action. A small displacement of the robot results in a minor alteration in the
cable shape and subsequently affects the feature parameters. Exploiting this observation,
the authors were able to linearize the deformation model. With the deformation model, a
velocity control law was applied to the robot, enabling the deformation of the cable into
the desired shape. The effectiveness of the method was successfully confirmed through
validation using two grippers that manipulated different cables to achieve the desired
shape.

2.5.1.5 3D positioning of feature points

Another approach for shape servoing of soft objects involving an online estimation of
the deformation Jacobian was introduced by [Lagneau et al., 2020a]. The shape of the
object is represented through 3D points of interest and the Jacobian estimation employs
weighted least-squares minimization within a sliding window. Subsequently, the obtained
deformation Jacobian is employed in the control law to manipulate the object shape.

In comparison to alternative model-free approaches, the study conducted by [Lagneau
et al., 2020a] employed a limited number of parameters for online estimation of the defor-
mation Jacobian. By contrasting their method with the approach utilizing the Broyden
update rule, the authors concluded that their technique exhibits greater resistance to
noise. This is achieved by updating each row of the deformation Jacobian matrix based
on a confidence criterion derived from the observed deformations. Specifically, in their
approach, the Jacobian is updated only for degrees of freedom (DOFs) that meet the
confidence criterion. This confidence criterion serves as a filter for noisy measurements
when the system approaches the desired deformation. The update of the estimated defor-
mation Jacobian is determined by a user-defined confidence threshold, and the confidence
criterion is determined using eigenvalue decomposition to ascertain how the Jacobian will
be updated.

Multiple experiments were conducted to verify the capability of controlling the shape of soft objects. In [Lagneau et al., 2020a], they presented various shape-serving tasks involving 3D objects, where they derived specific points of interest on the object surface to achieve desired positions. The authors demonstrated that their technique remains independent of the quantity of points of interest and external disruptions. The experiments revealed that the proposed method achieved comparable accuracy to state-of-the-art methods, all while relying on a smaller number of parameters to tune. Furthermore, they validated in [Lagneau et al., 2020b] their approach by actively controlling the 3D shape of a wire, using two robotic manipulators to control its two extremities and shape it to reach a desired 3D form.

2.5.1.6 Conclusion

In this section, we presented various model-free approaches that focus on deformable
object manipulation. The common idea among these methods is their attempt to esti-
mate the Jacobian matrix, which relates the variations of the visual feedback features
to the manipulator motions, purely from visual observation and robot odometry. These
approaches do not rely on a physics-based model of the deformable object, and in some
cases, they do not require knowledge of the camera calibration used to capture the feed-
back information. However, these methods can be sensitive to noise since their Jacobian
matrix is directly derived from inherently noisy visual data. To address this issue, some
approaches employed data filtering techniques before constructing the Jacobian matrix,
while others design robust controllers to mitigate measurement noise.

In addition to the mentioned Jacobian-based approaches, there exist numerous learning-
based approaches that have been suggested for deriving the deformation model or manip-
ulation strategies directly from data. Here, we will mention a few of them. For instance,
[Schulman, Gupta, et al., 2013] used depth images to encode the status of a manipulation
task of soft objects and learning from demonstration to allow a robot to execute multi-
step manipulation strategies. This technique has been expanded upon using reinforcement
learning [Hadfield-Menell et al., 2015] and tangent space mapping [T. Tang et al., 2016].
A deep learning-based end-to-end framework has also been proposed in [Yang et al., 2016].
Furthermore, [Clegg et al., 2017] used reinforcement learning for haptic manipulation. In
this case, the use of a high-quality simulator was crucial for the success of the learning
process. Finally, [Z. Hu et al., 2019] employed deep learning to estimate the deformation
Jacobian matrix. This matrix was subsequently used to generate a control law aimed at
controlling both the position and shape of deformable objects.

2.5.2 Model-based

Model-based deformable object manipulation refers to the process of manipulating objects using their physics-based models. A physics-based model is used to simulate the actual behavior of the object when exposed to external forces. In this section, we review different approaches on soft object shape servoing based on object models.

2.5.2.1 Indirect positioning of object feature points

Some pioneering studies have addressed the indirect positioning of soft objects, where
multiple selected object points, denoted feature points, are driven to their desired positions
indirectly by manipulating other points of the object, denoted as manipulated points.

One of the first works that proposed a control approach for simultaneously position-
ing multiple points on deformable textile fabrics was presented by [Wada et al., 1998].
The operation points in Figure 2.10(a) correspond to the manipulated points, while the
positioned points refer to the feature points. Their method involved creating a simplified
physics-based model using a two-dimensional MSM system with different spring stiff-
nesses and initial lengths along horizontal, vertical, and diagonal axes, as presented in
Figure 2.10(b).

By employing visual sensors and a linearized model in a quasi-static equilibrium state,
their control method effectively guided the positioned points to desired locations, even
without precise knowledge of the fabric physical properties.

Subsequently, in the work conducted by Kinio and Patriciu [Kinio & Patriciu, 2012],
the object was modeled using the FEM and a H∞ controller was employed for the purpose
of accomplishing the positioning task. This study involved two manipulated points and
two feature points. The authors demonstrated via simulation and experimental results
that in this case the H∞ controller presents a more advantageous design option compared
to PID (Proportional-Integral-Derivative) controllers.
<figure><img src="image_18.png" /><figcaption>Figure 2.10: Indirect positioning of multiple feature points, from [Wada et al., 1998]: (a) and (b) depict respectively the distribution of mesh points and the spring model of the 2D object.</figcaption></figure>

### 2.5.2.2 Two-phase strategy for object modeling and shaping in one dimension

A two-phase strategy for shape servoing, which is depicted in Figure 2.11(a), was introduced by [Higashimori et al., 2010]. In the first phase, the object is modeled using a four-element viscoelastic model that included two elements for elasticity and two for viscosity, then its parameters are estimated. The estimation process involved applying a contact force to the object in a step-wise manner, as shown in Figure 2.11(b). The force was increased up to a threshold level to avoid damaging the object, maintained for a specific duration, gradually reduced to zero, and finally, the contact was maintained without applying any force indefinitely. The resulting deformations during non-zero force application represented the elastic part (red plot), while the maintained deformations after unloading represented the plastic part (green plot). In the second phase, the authors actively controlled the plastic deformation response to achieve the desired final deformation. By using the estimated viscoelastic and plastic parameters, they determined the necessary force required to achieve the desired shape of the object.

Experimental validation was conducted to deform a 3D parallelepipedic object in one direction.

### 2.5.2.3 2D object contour servoing using MSM

[Das & Sarkar, 2011] proposed an approach for manipulating the shape of deformable objects, using a system of manipulators. The initial and final shapes of the object are defined by boundary curves. The object is modeled using a MSM where the connections between nodes are modeled as a Voigt model, as shown in Figure 2.12(a). The Voigt
<figure><img src="image_19.png" /><figcaption>Figure 2.11: 1D shape deformation of an unknown soft object, from [Higashimori et al., 2010]. In (a), the two-phase strategy for shape servoing of an unknown rheological object is presented. In (b), the decomposition of the object deformation response is illustrated, with $x_1$ representing the elastic response and $x_2$ representing the plastic response.</figcaption></figure>

model consists of a parallel arrangement of a spring and a damper. To achieve the desired shape change, an optimization-based planner is designed, which minimizes an energy-like criterion to determine the optimal locations of the contact points along the final shape curve. Importantly, each manipulator operates independently without the need for communication between them. Additionally, a robust controller was derived to handle modeling uncertainties.

Two simulation results provided evidence of the effectiveness of the proposed method. The first simulation involved applying the shape controller to deform a circular shape into either an ellipse or a square, using varying numbers of actuation points. The outcomes indicated that a higher number of actuation points resulted in improved system accuracy. In the second simulation, a deformable object was manipulated by a three-finger robotic hand without causing deformation. The robust controller effectively achieved this objective. Figure 2.12(b) illustrates an example of the first simulation task.

### 2.5.2.4 FEM-based soft robot control

An alternative approach for controlling the shape of soft objects was proposed in the context of soft robots. For example, the robot needs to deform itself to maintain specific configurations and accomplish given tasks. In this context, an approach that includes two controllers (see Figure 2.13) for the manipulation of soft robots, utilizing visual tracking and a simulation-based predictor, was introduced in [Z. Zhang et al., 2017].
<figure><img src="image_20.png" /><figcaption>Figure 2.12: 2D shape deformation of a soft object presented in [Das & Sarkar, 2011]. In (a), a model of a 2D object is presented. In (b), the shape transformation from a circle to an ellipse is represented using 12 manipulators. The initial shape, desired shape, and final shape are indicated by blue dashed, black dotted, and red solid lines, respectively.</figcaption></figure>

The first controller, depicted in Figure 2.13(a), is designed to guide the soft robot end effector to a desired position by applying a control law. This controller utilizes a simulation-based predictor that calculates a Jacobian matrix using the kinematic model of the soft robot, where a real-time FEM is employed. The open-loop controller consists of two parts: closed-loop control of the simulation model and open-loop control of the soft robot. The closed-loop control ensures that the Jacobian matrix of the simulation model remains close to that of the real robot, while the open-loop control ensures that the soft robot end effector tracks the desired trajectory. It is worth noting that the first controller is still considered open-loop, even with the presence of closed-loop control, as it does not require feedback from the soft robot. In this approach, both the simulation model and the soft robot utilize the same control input, which is computed based on the contribution of the actuators derived from the simulation model.

The second controller (see Figure 2.13(b)) utilizing infrared cameras and feedback is developed to correct the position of the end effector. This controller is implemented to actuate the simulation model, aligning the end effector of the simulation model with that of the real robot. Consequently, both the simulation model and the soft robot exhibit similar configurations.

To evaluate these methods, the authors conducted experiments on a soft robot actuated by cables.
<figure><img src="image_21.png" /><figcaption>Figure 2.13: The implementation schema of the two controllers proposed in [Z. Zhang et al., 2017]: (a) and (b) correspond to the implementation of the open-loop controller and closed-loop controller, respectively.</figcaption></figure>

### 2.5.2.5 FEM-based 3D deformable object elasticity parameters estimation and open-loop shape control

Another approach based on FEM to generate open-loop forces on anthropomorphic fingertips and achieve desired displacement on a soft object was proposed by [Ficuciello et al., 2018]. The methodology is depicted in Figure 2.14. Instead of using objects with known parameters, they employed a vision system and a force sensor to estimate the Young's modulus *E* and Poisson's ratio *ν*. This estimation relied on the tracking technique introduced by [Petit, Lippiello, Fontanelli, et al., 2017], which minimizes the fitting error between simulated and RGB-D sensor-captured deformations. The authors tackled the nonlinear optimization problem of parameter estimation using the gradient-free Nelder-Mead method [Nelder & Mead, 1965].

Once the model parameters were estimated, the objective is to minimize the deviation between end-effector positions and desired positions. This task was accomplished using an online inverse simulation technique, employing convex optimization, for adjusting the applied forces at the fingertips during contact and achieving the desired deformation.

Additionally, since the deformation was performed in an open-loop fashion, there was an inherent error between the expected and actual deformations on the real object. The authors identified three sources contributing to this error: convex optimization, model parameter approximations, and hardware limitations such as sensor uncertainties. The experimental results presented the deformation achieved by the underactuated anthropomorphic SCHUNK 5-Finger Hand (S5FH) on a cylindrical object.
<figure><img src="image_22.png" /><figcaption>Figure 2.14: The pipeline for an open-loop deformation control using FEM, from [Ficu-
ciello et al., 2018].</figcaption></figure>

2.5.2.6 ARAP-based planar object shape servoing

An alternative method for shape servoing, wherein a geometric mesh was utilized to represent the object surface instead of a physical mesh was introduced in [Shetab-Bushehri et al., 2022]. The authors stated that in numerous practical scenarios of deformable object manipulation, the object tends to naturally maintain local rigidity. Based on this observation, the authors introduced a shape-servoing approach based on the As-Rigid-As-Possible (ARAP) modeling [Sorkine & Alexa, 2007]. The authors further justified their motivation by highlighting that this type of modeling eliminates the need to estimate mechanical parameters of the object. The authors proposed a control scheme that aimed to guide points of a planar object to desired positions through a closed-loop system using the complete object surface points as a feedback measurement. To obtain the object points, the authors therefore tracked the object during the deformation.

The proposed control scheme integrates an estimation of the deformation Jacobian
matrix based on the ARAP model. The Jacobian was approximated by establishing a
connection between the robot end-effector frame and the ARAP mesh. In this approach,
the changes in the features correspond to the variations in the model surface points, while
the velocities correspond to the applied motions on the object surface by the robot end-
effectors. Then, to obtain an initial approximation of the Jacobian matrix, the undeformed
shape of the object was taken as a starting point, and a perturbation was simulated in
one degree of freedom (DOF). This process results in a new stable shape for the object
after the perturbation. Subsequently, a standard forward finite differences method was utilized to estimate the column of the Jacobian associated with the perturbed DOF. By repeating this procedure for all DOFs, an initial approximation of the Jacobian matrix was obtained.

Regarding the shape tracking pipeline, it was developed using the approach introduced in [Aranda et al., 2020]. This pipeline utilizes monocular vision and is built upon the principles of shape-from-template. It relies on the presence of visual texture on the object surface for accurate matching and assumes that the object surface undergoes isometric deformation.

The authors made the assumption that the object shape remains statically stable, achieving quasi-static equilibrium, and responds smoothly to robot motions. They conducted several shape servoing tasks using two Franka Emika robot arms and four different objects made of various materials to validate the effectiveness of their proposed scheme. An example of a shape servoing task accomplished in this work is presented in Figure 2.15. To ensure that the desired shape could be achieved, each experiment started by manually manipulating the two arms holding the object to record the desired shape. Afterward, the deforming object was manually adjusted to its initial shape using the two robot arms before starting the task.

### 2.5.2.7 Conclusion

In this section, we presented a review of various model-based approaches for manipulating soft objects, along with a related work concerning the deformation of soft robots. The control laws developed to accomplish these tasks depend on the object models employed to approximate the physical behavior of the object. These approaches either assume that the model parameters are known in advance or involve estimating them before the manipulation process.

## 2.5.3 Conclusion

Tables 2.1 and 2.2 present a summary of different methods presented in this section and used to deform soft objects.

Model-free approaches can be categorized into two distinct groups. In the first category, known as numerically estimated Jacobian-based approaches, some authors argue that there is no need to estimate the object model or its physical parameters. Neverthe-
<figure><img src="image_23.png" /><figcaption>Figure 2.15: Illustration of a shape servoing task using the ARAP model, from [Shetab-Bushehri et al., 2022]: (a) and (b) represent the camera view of the object before and after the deformation task, respectively. (c) and (d) illustrate the object surface points before and after the deformation task, respectively. The blue points correspond to the actual object surface points, and the red ones correspond to their desired positions.</figcaption></figure>
less, other methods within this category attempted to estimate these parameters during the deformation process. These numerically estimated Jacobian-based models primarily consider the quasi-static equilibrium of the object, limiting their applicability to dynamic objects. Furthermore, since these approaches are mainly based on the numerical estimation of a deformation Jacobian from data obtained from sensors, they are particularly sensitive to measurement noise. In the second category, machine learning techniques have garnered attention, but they require large datasets for effective model training.

In contrast to model-free approaches, there are model-based methods that rely on physics-based models. These approaches not only take into account the object geometric shape, often described by a 3D mesh, but also consider its mechanical properties. These models aim to replicate the actual behavior of the object when subjected to external forces. However, it is worth noting that these approaches require precise knowledge of the object model parameters. These parameters have to be known in advance, or accurately identified before the manipulation task.

## 2.6 Positioning of this thesis in relation to the existing literature

In the preceding sections, we provided an overview of the state-of-the-art in deformable object modeling, tracking, and deformation. Our purpose in presenting these topics was to enquire the reader with recent developments in the literature and to provide a clearer introduction to the objectives of this thesis. In this thesis, we aim to introduce novel physics-based methods for deformable manipulation, utilizing a simple mass-spring model.

Unlike the previously mentioned approaches in Section 2.5, we introduce a novel vision-based and physics-based control law for shaping deformable objects, aiming to employ a simple model that has the advantage of a low computational cost to enable fast shape servoing tasks. The motivation for combining vision-based and physics-based approaches arises from advancements in computer graphics, where various physics-based models were proposed to represent deformable objects, and from computer vision, which aids in tracking such objects, ultimately reducing the gap between the object simulated deformation and its real one. In our approach, the vision component plays a crucial role, where the object is tracked during deformation. If the object model deviates from the actual object deformation, we adjust the object model by applying constraints to it. Regarding the object model, it imitates the physical behavior of the real object, with its associated
<table><caption>Table 2.1: Summary of deformable object manipulation approaches (Part I)</caption><thead><tr><th>Approach</th><th>Category</th><th>Dimensionality</th><th>Feedback information</th><th>Hypotheses/Notes</th></tr></thead><tbody><tr><td>[Wada et al., 1998]</td><td>Model-based (MSM)</td><td>Planar objects</td><td>Object 2D points</td><td>Quasi-static assumption</td></tr><tr><td>[Higashimori et al., 2010]</td><td>Model-based (four-element viscoelastic model)</td><td>3D objects</td><td>Object length</td><td>Deformation in one direction</td></tr><tr><td>[Das & Sarkar, 2011]</td><td>Model-based (MSM)</td><td>Planar objects</td><td>Object contour points</td><td></td></tr><tr><td>[Kinio & Patriciu, 2012]</td><td>Model-based (FEM)</td><td>Planar objects</td><td>Object 2D points</td><td></td></tr><tr><td>[Navarro-Alarcon et al., 2013]</td><td>Model-free</td><td>3D objects</td><td>2D deformation features</td><td>Quasi-static assumption</td></tr><tr><td>[Navarro-Alarcon et al., 2014]</td><td>Model-free</td><td>Linear and planar objects</td><td>2D deformation features</td><td>Quasi-static assumption and small deformations</td></tr><tr><td>[Navarro-Alarcon et al., 2016]</td><td>Model-free</td><td>3D objects</td><td>2D position and deformation features</td><td>Quasi-static assumption and slow deformations</td></tr><tr><td>[Navarro-Alarcon & Liu, 2017]</td><td>Model-free</td><td>3D objects</td><td>Object 2D contour parameterized using Fourier series</td><td>Quasi-static assumption</td></tr><tr><td>[Z. Zhang et al., 2017]</td><td>Model-based (FEM)</td><td>Soft robots</td><td>Robot configuration</td><td>Quasi-static assumption</td></tr><tr><td>[Ficuciello et al., 2018]</td><td>Model-based (FEM)</td><td>3D objects</td><td>End-effector positions</td><td>Quasi-static assumption and open-loop deformation</td></tr></tbody></table>
<table><caption>Table 2.2: Summary of deformable object manipulation approaches (Part II)</caption><thead><tr><td>Approach</td><td>Category</td><td>Dimensionality</td><td>Feedback information</td><td>Hypotheses/Notes</td></tr></thead><tbody><tr><td>[J. Zhu et al., 2018]</td><td>Model-free</td><td>Linear objects</td><td>Cable pixel coordinates parameterized using Fourier series</td><td></td></tr><tr><td>[Lagneau et al., 2020a]</td><td>Model-free</td><td>3D objects</td><td>Object 3D points</td><td></td></tr><tr><td>[Lagneau et al., 2020b]</td><td>Model-free</td><td>Linear objects</td><td>Cable 3D geometric features based on B-splines</td><td>Valid for small deformations and adapted for larger ones</td></tr><tr><td>[Qi et al., 2021]</td><td>Model-free</td><td>Linear and 3D objects</td><td>Object 2D contour parameterized using image moments</td><td>Quasi-static assumption and slow deformations</td></tr><tr><td>[Shetab-Bushehri et al., 2022]</td><td>Model-based (ARAP)</td><td>Planar objects</td><td>Object 3D points</td><td>Quasi-static assumption and shape local rigidity.</td></tr></tbody></table>
parameters obtained prior to the deformation process through a few simple steps.

In contrast to approaches that solely rely on models, particularly model-based meth-
ods, our approach utilizes a simple mass-spring model and does not require precise knowl-
edge of the exact parameters of the soft object. Compared to model-free approaches, we
analytically develop the formulation that provides the variation of the object deformation
in function of the motion of one or multiple robotic manipulators using this simple mass-
spring model. Moreover, our approach does not rely on the assumption of quasi-static
equilibrium that is generally considered in the existing methods of the state-of-the-art. As
a result, our developed control law will be analytically derived, taking into consideration
the dynamic nature of the object.

Additionally, our approach differs from recent model-free methods in how object fea-
ture points and shape are represented thanks to the use of a RGB-D camera. Regarding
feature points, instead of relying solely on their 2D projections in the sensor image (pixels),
we use their 3D positions. Concerning object shape representation, we avoid describing it
solely through the projection of its contour, which heavily depends on the camera pose rel-
ative to the object model, or by using specific 3D points from the object surface. Instead,
we represent 3D objects using their complete set of surface points, while for 2D objects,
we define the shape through their 2D contour points. Then, we represent the object shape
using a low-dimensional feature vector. Unlike the use of 2D Fourier descriptors or 2D
moments for 3D objects, as proposed in the state-of-the-art approaches, we employ 3D
Fourier descriptors and 3D moments for 3D objects and 2D Fourier descriptors for 2D
objects. To the best of our knowledge, defining the shape of an object by its complete
set of 3D surface points and considering a dynamic object without the assumption of
quasi-static equilibrium has not been explored by other authors. Furthermore, to the
best of our knowledge, the uses of 3D Fourier descriptors and 3D moments have not been
exploited before in the context of deformable manipulation.

Finally, Figure 2.16 illustrates a general overview of our proposed closed-loop visual
control system. It incorporates a simulator emulating the dynamics of the soft object
(green box), a tracking process to correct the modeled deformation (blue box), and a
closed-loop controller (pink box).
<figure><img src="image_24.png" /><figcaption>Figure 2.16: Block diagram of the approach: visual tracking part in blue, physics-based simulator in green, and closed-loop control scheme in red.</figcaption></figure>
<figure><img src="image_25.png" /></figure>
# METHODOLOGY

This chapter presents the derivation of our first approach for deformable object manipulation, which consists of indirectly positioning multiple object points to their desired positions by acting on distant manipulated points. Within this approach, the model points of interest are referred to as feature points, and the object is approximated using a mass-spring model (MSM). An illustration of the desired task is depicted in Figure 3.1. Here, we develop a control law to be applied to the manipulated points to guide the feature points towards their desired positions.

We begin in Section 3.1 by establishing the analytical relationship that relates the displacements of the feature points to the successive motions of the manipulated points. This relationship forms the foundation of our contribution, built upon the MSM while accounting for the propagation delay introduced by it. We start in Section 3.1.1 by formulating the displacements of points within a 1D mesh in response to a series of motions applied to one or two manipulated points. Subsequently, in Section 3.1.2, we expand our formulation scope to encompass a 3D mesh involving multiple manipulated points. Following this, Section 3.1.3 introduces the designed control law.

Then, Section 3.2 presents simulation results that demonstrate the robustness of the proposed closed-loop system. Since our approach is model-based, several factors can influence the efficiency of the control law. One critical consideration is the mesh resolution, which refers to the number of points used to represent the model, which is a significant consideration in any numerical methods. Another important aspect is the model parameters, as it may not always be possible to predict the exact parameters of the soft objects that enable a perfect imitation of the real deformation. To account for these factors, we will include these considerations in our analysis of the proposed approach robustness. Furthermore, this section includes additional results obtained by applying the proposed control law to the system, whether it is fully actuated or not. It also includes comparisons between the proposed approach (PA) and other model-free methods [Navarro-Alarcon et al., 2013; J. Zhu et al., 2018].
Numerical simulations are used to validate the proposed approach for the indirect positioning task, where multiple feature points are considered. In addition, the proposed approach is applicable for the indirect positioning of either single or multiple feature points in real experiments involving real soft objects and robotic manipulators, as outlined in Chapter 4. Furthermore, the pipeline can be adjusted for comprehensive shape servoing, involving the deformation of the entire object surface, as explored in Chapter 5.

<figure><img src="image_26.png" /><figcaption>Figure 3.1: An example of indirect positioning of 3 feature points using 3 manipulated points.</figcaption></figure>

## 3.1 Deformable object servoing

In this modeling section, we assume that the object model is known with its corresponding parameters.

Moreover, the object is modeled using the mass-spring-model (MSM), wherein its dynamic behavior is described by (2.11). In all following analytical derivations related to the modeling and design of the control law, we neglect the external forces $f_e$ since they are a priori unknown. This exclusion includes the gravity force, which, although known, is not considered. This is because the object model is created when the object is in its initial rest state, thereby compensating for both the ground force exerted by the supporting surface and the effects of gravity. However, $f_e$ is involved in the real experiments, presented in Section 4.2.2, and we will see that these external forces do not perturb the stability and robustness of our system. Consequently, (2.11) is simplified to:

$$m_i \ddot{x}_i = f_{si} + f_{Di} \qquad (3.1)$$
where $f_D$ represents the damping force with damping value $D_v$. Regarding $f_s$, it corresponds to the spring force acting on $P_i$, which depends on the stiffness $K_{ij}$ connecting $P_i$ to its neighbors, $P_j$, $\forall j \in \nu_i \subset \mathcal{N}$, and expressed as:

$$f_{si} = \sum_{j \in \nu_i} K_{ij} (\|\mathbf{x}_i - \mathbf{x}_j\| - l_{ij}^0) \frac{(\mathbf{x}_j - \mathbf{x}_i)}{\|\mathbf{x}_j - \mathbf{x}_i\|} = \sum_{j \in \nu_i} \alpha_{ij} K_{ij} (\mathbf{x}_j - \mathbf{x}_i) \quad (3.2)$$

Subsequently, we can determine the model point positions at any time $(t + dt)$ by utilizing the semi-implicit Euler integration method, given that we have knowledge of their positions at time $(t)$, with $dt$ representing the simulation time step. We opt for the semi-implicit Euler integration method due to its simplicity and computational efficiency, as detailed in [Bhasin & Liu, 2006]. Finally, the model update is presented as follows:

$$\mathbf{x}_{i(t+dt)} = \mathbf{x}_{i(t)} + (dt - \frac{dt^2}{m_i} D_v) \dot{\mathbf{x}}_{i(t)} + \frac{dt^2}{m_i} f_{si}(t) \quad (3.3)$$

with

$$\dot{\mathbf{x}}_{i(t)} = (\mathbf{x}_{i(t)} - \mathbf{x}_{i(t-dt)})/dt \quad (3.4)$$

### 3.1.1 Simple 1D MSM

An example of a 1D mesh with $N$ nodes is represented in Figure 3.2. In a first part, we consider that the first point, $P_1$, is manipulated while the other ones obey the dynamics provided in (3.1). In a second part, we will consider the case of two manipulated points.

<figure><img src="image_27.png" /><figcaption>Figure 3.2: An example of a 1D mesh with $N$ nodes connected by springs, including one manipulated point depicted in red.</figcaption></figure>

#### 3.1.1.1 One manipulated point

For illustrative purpose, we first derive the displacements of $P_2$ due to successive displacements applied on $P_1$. This will necessitate to also determine the position of $P_3$. We then generalize the obtained results to determine the position and velocity of any point in the mesh.
Let us consider that $P_1$ is manipulated such that its successive positions are given by:

$$x_{1(t_1)} = x_{1(t_0)} + \Delta_{1(t_1)} \text{ with } t_1 = t_0 + dt,$$

$$x_{1(t_2)} = x_{1(t_1)} + \Delta_{1(t_2)} \text{ with } t_2 = t_1 + dt,$$

$$x_{1(t_3)} = x_{1(t_2)} + \Delta_{1(t_3)} \text{ with } t_3 = t_2 + dt,$$

where $x_{1(t_0)}$ is the position of $P_1$ in its initial rest state and $\Delta_{1(t)}$ are the successive displacements applied on $P_1$ at each time ($t$) (see the red circles in Figure 3.3).

<figure><img src="image_28.png" /><figcaption>Figure 3.3: An illustration of the propagation delay introduced by the MSM in a 1D mesh when one manipulated point is considered.</figcaption></figure>

From (3.3) and (3.2), the positions of $P_2$ and $P_3$ are expressed as follows:

$$x_{2(t+dt)} = x_{2(t)} + \left(dt - \frac{dt^2}{m_2}D_v\right)\dot{x}_{2(t)} + \frac{dt^2}{m_2}f_{s21(t)} + \frac{dt^2}{m_2}f_{s23(t)} \quad (3.5)$$

$$x_{3(t+dt)} = x_{3(t)} + \left(dt - \frac{dt^2}{m_3}D_v\right)\dot{x}_{3(t)} + \frac{dt^2}{m_3}f_{s32(t)} + \frac{dt^2}{m_3}f_{s34(t)} \quad (3.6)$$

We thus deduce, as illustrated in Figure 3.3, that:

**For $t = t_1$:** $P_2$ stays at rest with $x_{2(t_1)} = x_{2(t_0)}$ and $\dot{x}_{2(t_1)} = 0$ since $\dot{x}_{2(t_0)} = f_{s21(t_0)} = f_{s23(t_0)} = 0$. Similarly, $P_3$ also stays at rest. As shown in Figure 3.3, the MSM induces a propagation delay characterized by the fact that $P_2$ starts moving at $t_2 = t_1 + dt$ due to the motion applied to $P_1$ at $t_1$, whereas $P_3$ starts moving at time $t_3 = t_1 + 2dt$, etc.
For $t_2 = t_1 + dt$: $P_3$ is still at rest while $P_2$ starts moving since $f_{s21(t_1)} \neq 0$. More precisely, we have from (3.2)

$$f_{s21(t_1)} = -K_{12}(l_0 - \Delta_{1(t_1)} - l_0) = K_{12}\Delta_{1(t_1)}$$

from which we obtain:

$$x_{2(t_2)} = x_{2(t_1)} + \gamma_{21(t_2)}^{(0)} \Delta_{1(t_1)}, \quad (3.7)$$

$$\text{with } \gamma_{21(t_2)}^{(0)} = \frac{dt^2}{m_2} K_{12}$$

The term $\gamma_{ij(t)}^{(q)}$ will always appear in the following to represent the propagation coefficient denoting the ratio between the displacement of any node $P_i$ at time $t$ and the displacement of any manipulated point $P_j$ at time $t_1 + qdt$, whatever the mesh and number of manipulated points.

For $t_3 = t_2 + dt$: we obtain from (3.2) and (3.7):

$$\begin{aligned} f_{s21(t_2)} &= -K_{12}[l_0 + \gamma_{21(t_2)}^{(0)}\Delta_{1(t_1)} - (\Delta_{1(t_1)} + \Delta_{1(t_2)}) - l_0] \\ &= K_{12}[(1 - \gamma_{21(t_2)}^{(0)})\Delta_{1(t_1)} + \Delta_{1(t_2)}] \end{aligned} \quad (3.8)$$

$$\begin{aligned} f_{s23(t_2)} &= K_{23}(2l_0 - (l_0 + \gamma_{21(t_2)}^{(0)}\Delta_{1(t_1)}) - l_0) \\ &= -K_{23}\gamma_{21(t_2)}^{(0)}\Delta_{1(t_1)} \end{aligned} \quad (3.9)$$

Then, by using (3.5) and (3.7) in (3.4), we get:

$$\begin{aligned} x_{2(t_3)} = x_{2(t_2)} &+ \left(1 - \frac{dt}{m_2}D_v\right) \gamma_{21(t_2)}^{(0)} \Delta_{1(t_1)} + \frac{dt^2}{m_2} K_{12} \left[(1 - \gamma_{21(t_2)}^{(0)})\Delta_{1(t_1)} + \Delta_{1(t_2)}\right] \\ &- \frac{dt^2}{m_2} K_{23} \gamma_{21(t_2)}^{(0)} \Delta_{1(t_1)} \end{aligned} \quad (3.10)$$

which can be rewritten as:

$$x_{2(t_3)} = x_{2(t_1)} + \gamma_{21(t_3)}^{(0)} \Delta_{1(t_1)} + \gamma_{21(t_2)}^{(1)} \Delta_{1(t_2)} \quad (3.11)$$
with

$$ \gamma_{21}^{(0)}(t_3) = \gamma_{21}^{(0)}(t_2) + \left(1 - \frac{dt}{m_2}D_v\right)\gamma_{21}^{(0)}(t_2) + \frac{dt^2}{m_2}K_{12}\left(1 - \gamma_{21}^{(0)}(t_2)\right) - \frac{dt^2}{m_2}K_{23}\gamma_{21}^{(0)}(t_2) \quad (3.12) $$

$$ \gamma_{21}^{(1)}(t_2) = \frac{dt^2}{m_2} K_{12} \quad (3.13) $$

Similarly, since $f_{s32}(t_2) = -f_{s23}(t_2)$, we obtain for $P_3$:

$$ \begin{aligned} x_{3(t_3)} &= x_{3(t_1)} + \frac{dt^2}{m_3} K_{23} \gamma_{21}^{(0)}(t_2) \Delta_1(t_1) \\ &= x_{3(t_1)} + \gamma_{31}^{(0)}(t_3) \Delta_1(t_1) \end{aligned} \quad (3.14) $$

with

$$ \gamma_{31}^{(0)}(t_3) = \frac{dt^2}{m_3} K_{23} \gamma_{21}^{(0)}(t_2) \quad (3.15) $$

For $t_4 = t_3 + dt$: we obtain by following the same reasoning as before:

$$ x_{2(t_4)} = x_{2(t_1)} + \gamma_{21(t_4)}^{(0)} \Delta_1(t_1) + \gamma_{21(t_3)}^{(1)} \Delta_1(t_2) + \gamma_{21(t_2)}^{(2)} \Delta_1(t_3) \quad (3.16) $$

with

$$ \begin{aligned} \gamma_{21}^{(0)}(t_4) &= \gamma_{21}^{(0)}(t_3) + \left(1 - \frac{dt}{m_2}D_v\right)(\gamma_{21}^{(0)}(t_3) - \gamma_{21}^{(0)}(t_2)) + \frac{dt^2}{m_2}K_{12}\left(1 - \gamma_{21}^{(0)}(t_3)\right) \\ &\quad - \frac{dt^2}{m_2}K_{23}(\gamma_{21}^{(0)}(t_3) - \gamma_{31}^{(0)}(t_3)) \end{aligned} \quad (3.17) $$

$$ \gamma_{21}^{(1)}(t_3) = \gamma_{21}^{(1)}(t_2) + \left(1 - \frac{dt}{m_2}D_v\right)\gamma_{21}^{(1)}(t_2) + \frac{dt^2}{m_2}K_{12}\left(1 - \gamma_{21}^{(1)}(t_2)\right) - \frac{dt^2}{m_2}K_{23}\gamma_{21}^{(1)}(t_2) \quad (3.18) $$

$$ \gamma_{21}^{(2)}(t_2) = \frac{dt^2}{m_2} K_{12} \quad (3.19) $$

Indeed, we now have using (3.11) and (3.14):

$$ \begin{aligned} f_{s21}(t_3) &= -K_{12}[l_0 + \gamma_{21}^{(0)}(t_3)\Delta_1(t_1) + \gamma_{21}^{(1)}(t_2)\Delta_1(t_2) - (\Delta_1(t_1) + \Delta_1(t_2) + \Delta_1(t_3)) - l_0] \\ &= K_{12}[\Delta_1(t_3) + (1 - \gamma_{21}^{(0)}(t_3))\Delta_1(t_1) + (1 - \gamma_{21}^{(1)}(t_2))\Delta_1(t_2)] \end{aligned} $$

$$ \begin{aligned} f_{s23}(t_3) &= K_{23}[2l_0 + \gamma_{31}^{(0)}(t_3)\Delta_1(t_1) - (l_0 + \gamma_{21}^{(0)}(t_3)\Delta_1(t_1) + \gamma_{21}^{(1)}(t_2)\Delta_1(t_2)) - l_0] \\ &= -K_{23}[(\gamma_{21}^{(0)}(t_3) - \gamma_{31}^{(0)}(t_3))\Delta_1(t_1) + \gamma_{21}^{(1)}(t_2)\Delta_1(t_2)] \end{aligned} $$
which allows getting (3.16) using (3.3) and (3.11) in (3.4) as before.

For any time t written as $t = t_1 + kdt$: we can determine the position of $P_2$ in function of the successive displacements of $P_1$ by generalizing the above equations. We obtain:

$$x_{2(t)} = x_{2(t_1)} + \gamma_{21(t)}^{(0)} \Delta_{1(t_1)} + \gamma_{21(t-dt)}^{(1)} \Delta_{1(t_1+dt)} + \dots + \gamma_{21(t_2)}^{(k-1)} \Delta_{1(t-dt)} \quad (3.20)$$

where, as general form for $\gamma_{21(t)}^{(q)}$:

$$\gamma_{21(t)}^{(q)} = \gamma_{21(t-dt)}^{(q)} + \left(1 - \frac{dt}{m_2} D_v\right) \delta_{21(t-dt)}^{(q)} + \frac{dt^2}{m_2} K_{12} \left(1 - \gamma_{21(t-dt)}^{(q)}\right) \\ - \frac{dt^2}{m_2} K_{23} \left(\gamma_{21(t-dt)}^{(q)} - \gamma_{31(t-dt)}^{(q)}\right) \quad (3.21)$$

by denoting

$$\delta_{ij(t)}^{(q)} = \gamma_{ij(t)}^{(q)} - \gamma_{ij(t-dt)}^{(q)} \quad (3.22)$$

with $\gamma_{21(t_2)}^{(q)} = \frac{dt^2}{m_2} K_{12}, \forall q$ (which is coherent with (3.7), (3.12) and (3.19)) and $\gamma_{21(t_1)}^{(q)} = 0$ for initializing the iterative form.

Then, the velocity of $P_2$ can be deduced by evaluating (3.20) at time $t$ and $t - dt$:

$$\dot{x}_{2(t)} = (x_{2(t)} - x_{2(t-dt)})/dt \\ = \delta_{21(t)}^{(0)} \dot{x}_{1(t_1)} + \delta_{21(t-dt)}^{(1)} \dot{x}_{1(t_2)} + \dots + \delta_{21(t_2)}^{(k-1)} \dot{x}_{1(t-dt)}$$

which can be rewritten with the simple form:

$$\dot{x}_{2(t)} = \sum_{q=0}^{k-1} \delta_{21(t-qdt)}^{(q)} \dot{x}_{1(t_1+qdt)} \quad (3.23)$$

Using the same methodology outlined above, the velocity of any arbitrary point $P_n$ belonging to the 1D mesh presented in Figure 3.2 can be calculated in function of its corresponding propagation coefficient and the motion applied on $P_1$. We obtain by analogy with (3.23):

$$\dot{x}_{n(t)} = \sum_{q=0}^{k-1} \delta_{n1(t-qdt)}^{(q)} \dot{x}_{1(t_1+qdt)} \quad (3.24)$$

with $(jr_i - 1)$ the number of intermediate points between the manipulated point $P_j$ and the model point $P_i$ (for instance, $^1r_2 = 1$ since $P_2$ is directly attached to $P_1$ while $^1r_3 = 2$),
and, as a general form of the propagation coefficient $\gamma_{n1}^{(q)}$ that is involved in $\delta_{n1}^{(q)}$:

$$ \gamma_{n1(t)}^{(q)} = \gamma_{n1(t-dt)}^{(q)} + \left(1 - \frac{dt}{m_n} D_v\right) \delta_{n1(t-dt)}^{(q)} + \frac{dt^2}{m_n} K_{n,n} \left(\gamma_{n_1(t-dt)}^{(q)} - \gamma_{n1(t-dt)}^{(q)}\right) \\ - \frac{dt^2}{m_n} K_{n,n} \left(\gamma_{n1(t-dt)}^{(q)} - \gamma_{n_1(t-dt)}^{(q)}\right) \qquad (3.25) $$

where $P_{n_a}$ and $P_{n_b}$ are the previous and next neighboring points of $P_n$ (see Figure 3.2) and $\gamma_{n1(t_n)}^{(q)} = \frac{dt^2}{m_n} K_{n,n} \gamma_{n_1(t_{n-1})}^{(q)}$, with $t_n = t_1 + {}^1r_n dt$. For example, we have $\gamma_{41(t_4)}^{(q)} = \frac{dt^2}{m_2} \frac{dt^2}{m_3} \frac{dt^2}{m_4} K_{12} K_{23} K_{34}, \forall q$.

Note that in case it is no more the point $P_1$ that is manipulated but the other end-point $P_N$, $\dot{x}_{n(t)}$, $\delta_{nN(t)}^{(q)}$ and $\gamma_{nN(t)}^{(q)}$ are obtained by just replacing in (3.24) and (3.25) the terms $\dot{x}_1$, $\delta_{n1}^{(q)}$, $\gamma_{n1}^{(q)}$, $\gamma_{n_1}^{(q)}$, $\gamma_{n_2}^{(q)}$ and ${}^1r_n$ by $\dot{x}_N$, $\delta_{nN}^{(q)}$, $\gamma_{nN}^{(q)}$, $\gamma_{n_2N}^{(q)}$, $\gamma_{n_1N}^{(q)}$ and ${}^N r_n$.

The previous results are valid when only one manipulated point is considered in the 1D mesh. When multiple points are used to deform the object, the velocity of the model points is affected by all the external motions applied. This case is discussed further in the following part.

### 3.1.1.2 Two manipulated points

<figure><img src="image_29.png" /><figcaption>Figure 3.4: An example of a 1D mesh with 4 nodes connected by springs, including two manipulated points depicted in red.</figcaption></figure>

In the example depicted in Figure 3.4, $P_1$ and $P_4$ are the manipulated points whose successive positions are given at time $t = t_1 + kdt$ by:

$$ x_{1(t)} = x_{1(t_0)} + \sum_{q=0}^{k} \Delta_{1(t_1+qdt)} $$

$$ x_{4(t)} = x_{4(t_0)} + \sum_{q=0}^{k} \Delta_{4(t_1+qdt)} $$

At $t_2 = t_1 + dt$, $P_2$ moves only due to the displacement $\Delta_{1(t_1)}$ of $P_1$ while $P_3$ moves only due to the displacement $\Delta_{4(t_1)}$ of $P_4$ because of the propagation delay introduced by the
MSM. From the results obtained in the previous section, we directly obtain:

$$x_{2(t_2)} = x_{2(t_1)} + \gamma_{21}^{(0)}(t_2)\Delta_{1(t_1)}$$

$$x_{3(t_2)} = x_{3(t_1)} + \gamma_{34}^{(0)}(t_2)\Delta_{4(t_1)}$$

where $\gamma_{21}^{(0)}(t_2)$ and $\gamma_{34}^{(0)}(t_2)$ are obtained from the general form given in (3.25) (using $N = 4$ instead of $n$ and the other permutations mentioned above for computing $\gamma_{34}$).

At $t_3 = t_2 + dt$, the position of $P_2$ changes due to the displacements $\Delta_{1(t_1)}$ and $\Delta_{1(t_2)}$ of $P_1$ but also due to the displacement $\Delta_{4(t_1)}$ of $P_4$ that has modified the position of $P_3$ at $t_2$ (and vice versa for the position of $P_3$). This can be seen from the force $f_{s23(t_2)}$ obtained using (3.2) that now contains a supplementary term compared to (3.9):

$$\begin{aligned} f_{s23(t_2)} &= K_{23}[2l_0 + \gamma_{34}^{(0)}(t_2)\Delta_{4(t_1)} - (l_0 + \gamma_{21}^{(0)}(t_2)\Delta_{1(t_1)}) - l_0] \\ &= -K_{23}(\gamma_{21}^{(0)}(t_2)\Delta_{1(t_1)} - \gamma_{34}^{(0)}(t_2)\Delta_{4(t_1)}) \end{aligned}$$

Since the supplementary term is nothing but a sum proportional to $\Delta_{4(t_1)}$, by injecting this force and the other ones in (3.5) and (3.6), we obtain a very similar form as before for the position of $P_2$ and $P_3$. More precisely, we obtain:

$$x_{2(t_3)} = x_{2(t_1)} + \gamma_{21}^{(0)}(t_3)\Delta_{1(t_1)} + \gamma_{21}^{(1)}(t_2)\Delta_{1(t_2)} + \gamma_{24}^{(0)}(t_3)\Delta_{4(t_1)}$$

$$x_{3(t_3)} = x_{3(t_1)} + \gamma_{34}^{(0)}(t_3)\Delta_{4(t_1)} + \gamma_{34}^{(1)}(t_2)\Delta_{4(t_2)} + \gamma_{31}^{(0)}(t_3)\Delta_{1(t_1)}$$

By following exactly the same reasoning and computations as in the previous section, it is quite easy to generalize to a 1D mesh whose two endpoints are manipulated. For instance, when $P_1$ and $P_N$ in Figure 3.2 are deforming the object simultaneously, we obtain the velocity of any point $P_n$ under the form (note the analogy with (3.24)):

$$\dot{x}_{n(t)} = \sum_{q=0}^{k-1} \delta_{n1(t-qdt)}^{(q)} \dot{x}_{1(t_1+qdt)} + \sum_{q=0}^{k-N} \delta_{nN(t-qdt)}^{(q)} \dot{x}_{N(t_1+qdt)} \quad (3.26)$$

What we have developed so far reveals that the model point velocities are connected to the applied manipulated point motions through the proposed propagation coefficients. This holds true whether a single or two manipulated points deform the object simultaneously. In the subsequent paragraph, we explore the general case of a three-dimensional mesh with several manipulated points that deform the object simultaneously.
### 3.1.2 General 3D MSM

We now consider a general 3D mesh with M manipulators deforming the soft object simultaneously. Let us denote $P_{m_l}$, $l \in [1, ..., M]$ one of these manipulators. Its motion applied at time $t = kdt$ is given by $\Delta_{m_l(t)} = \dot{x}_{m_l(t)} dt$, and we consider $t_1 = 0$ for simplicity.

Following the methodology outlined in the previous section, we first consider a single manipulated point and determine the motion of mesh points in its neighborhood. We then generalize this result to all points in the mesh. Finally, we extend this formulation when all M manipulated points are active. As will be shown below, the main difficulty with respect to the previous 1D case lies in the fact that the forces $f_{s_{ij}}$ are no more linear with respect to $x_i$ and $x_j$, which will be solved using first-order Taylor decomposition.

#### 3.1.2.1 One manipulated point

For $t = dt$: For any point $P_i$ in the neighborhood of $P_{m_l}$, we obtain its position using (3.3) as follows:

$$ \mathbf{x}_{i(dt)} = \mathbf{x}_{i(0)} + \left( dt - \frac{dt^2}{m_i} D_v \right) \dot{\mathbf{x}}_{i(0)} + \frac{dt^2}{m_i} \sum_{j \in \nu_i} \mathbf{f}_{s_{ij}(0)} \quad (3.27) $$

where $\nu_i$ contains two parts: the manipulated point $P_{m_l}$ and the other mesh points $P_j$ in $\nu_i$ ($\nu_i - \{m_l\}$). By exploiting $f_{s_{ij}}$ in (3.27) we deduce:

$$ \sum_{j \in \nu_i - \{m_l\}} f_{s_{ij}(0)} = \sum_{j \in \nu_i - \{m_l\}} f_{s_{ij}}(\mathbf{x}_{i(0)}, \mathbf{x}_{j(0)}) \quad (3.28) $$

For all points that do not belong to the neighborhood of $P_{m_l}$, they remain at their initial positions, and consequently, the spring forces corresponding to these points are zero. This leads to:

$$ \sum_{j \in \nu_i - \{m_l\}} f_{s_{ij}(0)} = 0 \quad (3.29) $$

$$ \sum_{j \in \nu_i} f_{s_{ij}(0)} = f_{s_{im_l}}(\mathbf{x}_{i(0)}, \mathbf{x}_{m_l(0)} + \Delta_{m_l(0)}) \quad (3.30) $$

Since $f_{s_{im_l}}(\mathbf{x}_{i(0)}, \mathbf{x}_{m_l(0)}) = 0$, using a first-order approximation of $f_{s_{im_l}}$ in (3.30), we obtain:

$$ \sum_{j \in \nu_i} f_{s_{ij}(0)} = \frac{\partial f_{s_{im_l}}}{\partial \mathbf{x}_{m_l}} \Delta_{m_l(0)} \quad (3.31) $$
where, by denoting $\mathbf{d}_{ij} = (\mathbf{x}_j - \mathbf{x}_i)$, we know from (3.2) that:

$$ \frac{\partial \mathbf{f}_{s_{ij}}}{\partial \mathbf{x}_j} = \frac{\partial \mathbf{f}_{s_{ij}}}{\partial \mathbf{d}_{ij}} \frac{\partial \mathbf{d}_{ij}}{\partial \mathbf{x}_j} \left( = - \frac{\partial \mathbf{f}_{s_{ij}}}{\partial \mathbf{x}_i} \right) = \mathbf{K}_{ij} \mathbf{I}_{3 \times 3} - \mathbf{K}_{ij} l_0 \frac{\| \mathbf{d}_{ij} \|^2 \mathbf{I}_{3 \times 3} - \mathbf{d}_{ij} \mathbf{d}_{ij}^T}{\| \mathbf{d}_{ij} \|^3} \quad (3.32) $$

Therefore, using (3.31) in (3.27), we obtain:

$$ \mathbf{x}_{i(dt)} = \mathbf{x}_{i(0)} + \frac{dt^2 \partial \mathbf{f}_{s im_l}(\mathbf{x}_{i(0)}, \mathbf{x}_{m_l(0)})}{m_i \partial \mathbf{x}_{m_l}} \Delta_{m_l(0)} = \mathbf{x}_{i(0)} + \gamma_{im_l(dt)}^{(0)} \Delta_{m_l(0)} \quad (3.33) $$

where

$$ \gamma_{im_l(dt)}^{(0)} = \frac{dt^2 \partial \mathbf{f}_{s im_l}(\mathbf{x}_{i(0)}, \mathbf{x}_{m_l(0)})}{m_i \partial \mathbf{x}_{m_l}} \quad (3.34) $$

and $\gamma_{im_l(dt)}^{(j)}$ denotes the propagation coefficient that relates the motion of $P_i$ relative to $\Delta_{m_l(dt)}$ of $P_{m_l}$ at any time $t$.

Furthermore, for the other points $P_j \in \nu_i$ so that $m_l r_j = 2$, we have $\mathbf{x}_{j(dt)} = \mathbf{x}_{j(0)}$ and $\gamma_{j m_l(dt)}^{(0)} = \mathbf{0}$.

For $t = 2dt$: we have for $P_i$ in the neighborhood of $P_{m_l}$:

$$ \mathbf{x}_{i(2dt)} = \mathbf{x}_{i(dt)} + \left( dt - \frac{dt^2}{m_i} D_v \right) \dot{\mathbf{x}}_{i(dt)} + \frac{dt^2}{m_i} \sum_{j \in \nu_i} \mathbf{f}_{s_{ij}}(dt) \quad (3.35) $$

where, using (3.4) and (3.33):

$$ \dot{\mathbf{x}}_{i(dt)} = (\mathbf{x}_{i(dt)} - \mathbf{x}_{i(0)})/dt = \gamma_{im_l(dt)}^{(0)} \Delta_{m_l(0)}/dt \quad (3.36) $$

Concerning $\nu_i$, it is now decomposed in three parts (see Figure 3.5):

Case 1: for $j \in \nu_i - \{m_l\}$ and $m_l r_j = 2$, we have:

$$ \begin{aligned} \mathbf{f}_{s_{ij}(dt)} &= \mathbf{f}_{s_{ij}}(\mathbf{x}_{i(0)} + \gamma_{im_l(dt)}^{(0)} \Delta_{m_l(0)}, \mathbf{x}_{j(0)}) \quad (3.37) \\ &= \mathbf{f}_{s_{ij}}(\mathbf{x}_{i(0)}, \mathbf{x}_{j(0)}) + \frac{\partial \mathbf{f}_{s_{ij}}(\mathbf{x}_{i(0)}, \mathbf{x}_{j(0)})}{\partial \mathbf{x}_i} \gamma_{im_l(dt)}^{(0)} \Delta_{m_l(0)} \\ &= - \frac{\partial \mathbf{f}_{s_{ij}}(\mathbf{x}_{i(0)}, \mathbf{x}_{j(0)})}{\partial \mathbf{x}_j} \gamma_{im_l(dt)}^{(0)} \Delta_{m_l(0)} \end{aligned} $$
Case 2: for $j \in \nu_i - \{m_l\}$ and $^{m_l}r_j = 1$, we have:

$$
\begin{align}
\boldsymbol{f}_{s_{ij}(dt)} &= \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(0)} + \boldsymbol{\gamma}_{im_l(dt)}^{(0)} \boldsymbol{\Delta}_{m_l(0)}, \boldsymbol{x}_{j(0)} + \boldsymbol{\gamma}_{jm_l(dt)}^{(0)} \boldsymbol{\Delta}_{m_l(0)}) \tag{3.38} \\
&= \frac{\partial \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(0)}, \boldsymbol{x}_{j(0)})}{\partial \boldsymbol{x}_j} (\boldsymbol{\gamma}_{jm_l(dt)}^{(0)} - \boldsymbol{\gamma}_{im_l(dt)}^{(0)}) \boldsymbol{\Delta}_{m_l(0)} \nonumber
\end{align}
$$

Case 3: for $j = m_l$, we have:

$$
\begin{align}
\boldsymbol{f}_{s_{ij}(dt)} &= \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(0)} + \boldsymbol{\gamma}_{im_l(dt)}^{(0)} \boldsymbol{\Delta}_{m_l(0)}, \boldsymbol{x}_{j(0)} + \boldsymbol{\Delta}_{m_l(0)} + \boldsymbol{\Delta}_{m_l(dt)}) \tag{3.39} \\
&= \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(0)} + \boldsymbol{\gamma}_{im_l(dt)}^{(0)} \boldsymbol{\Delta}_{m_l(0)}, \boldsymbol{x}_{j(0)} + \boldsymbol{\Delta}_{m_l(0)}) + \frac{\partial \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(dt)}, \boldsymbol{x}_{j(dt)})}{\partial \boldsymbol{x}_j} \boldsymbol{\Delta}_{m_l(dt)} \nonumber \\
&= \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(0)}, \boldsymbol{x}_{j(0)}) + \frac{\partial \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(dt)}, \boldsymbol{x}_{j(0)} + \boldsymbol{\Delta}_{m_l(0)})}{\partial \boldsymbol{x}_j} \boldsymbol{\Delta}_{m_l(dt)} \nonumber \\
&\quad + \frac{\partial \boldsymbol{f}_{s_{ij}}(\boldsymbol{x}_{i(0)}, \boldsymbol{x}_{j(0)})}{\partial \boldsymbol{x}_j} (\boldsymbol{I}_{3 \times 3} - \boldsymbol{\gamma}_{im_l(dt)}^{(0)}) \boldsymbol{\Delta}_{m_l(0)} \nonumber
\end{align}
$$

<figure><img src="image_30.png" /><figcaption>Figure 3.5: An example of a 2D mesh consisting of 7 nodes. By selecting $Pm_l$ as $P_1$ and $P_i$ as $P_3$, we have $\nu_i = \{1, 2, 4, 5, 6, 7\}$, $^1r_2 = ^1r_3 = ^1r_4 = 1$ and $^1r_5 = ^1r_6 = ^1r_7 = 2$.</figcaption></figure>

By using (3.36), (3.37), (3.38) and (3.39) in (3.35), we obtain:

$$
\boldsymbol{x}_{i(2dt)} = \boldsymbol{x}_{i(0)} + \boldsymbol{\gamma}_{im_l(2dt)}^{(0)} \boldsymbol{\Delta}_{m_l(0)} + \boldsymbol{\gamma}_{im_l(dt)}^{(1)} \boldsymbol{\Delta}_{m_l(dt)} \quad (3.40)
$$
where:

$$ \gamma_{iml(dt)}^{(1)} = \frac{dt^2}{m_i} \frac{\partial \mathbf{f}_{sij}(\mathbf{x}_{i(dt)}, \mathbf{x}_{m_l(0)} + \mathbf{\Delta}_{m_l(0)})}{\partial \mathbf{x}_{m_l}} \quad (3.41) $$

$$ \gamma_{iml(2dt)}^{(0)} = \gamma_{iml(dt)}^{(0)} + \left(1 - \frac{dt}{m_i} D_v\right) \gamma_{iml(dt)}^{(0)} + \frac{dt^2}{m_i} \frac{\partial \mathbf{f}_{siml}(\mathbf{x}_{i(0)}, \mathbf{x}_{m_l(0)})}{\partial \mathbf{x}_{m_l}} (\mathbf{I}_{3 \times 3} - \gamma_{iml(dt)}^{(0)}) \\ + \frac{dt^2}{m_i} \sum_{j \in \nu_i - m_l} \frac{\partial \mathbf{f}_{sij}(\mathbf{x}_{i(0)}, \mathbf{x}_{j(0)})}{\partial \mathbf{x}_j} (\gamma_{jml(dt)}^{(0)} - \gamma_{iml(dt)}^{(0)}) \quad (3.42) $$

since $\gamma_{jml(dt)}^{(0)} = 0$ if $m_l r_j = 2$.

Furthermore, for the other points $P_j \in \nu_i$ so that $m_l r_j = 2$:

$$ \mathbf{x}_{j(2dt)} = \mathbf{x}_{j(dt)} + \left(dt - \frac{dt^2}{m_j} D_v\right) \dot{\mathbf{x}}_{j(dt)} + \frac{dt^2}{m_j} \sum_{o \in \nu_j} \mathbf{f}_{sjo(dt)} \quad (3.43) $$

We now have to consider the different cases for $\mathbf{f}_{sjo(dt)}$:

1. if $m_l r_o = 3$, we have:

$$ \mathbf{f}_{sjo(dt)} = \mathbf{f}_{sjo}(\mathbf{x}_{j(0)}, \mathbf{x}_{o(0)}) = \mathbf{0} $$

2. if $m_l r_o = 2$, we have:

$$ \mathbf{f}_{sjo(dt)} = \mathbf{f}_{sjo}(\mathbf{x}_{j(0)}, \mathbf{x}_{o(0)}) = \mathbf{0} $$

3. if $m_l r_o = 1$, we have:

$$ \mathbf{f}_{sjo(dt)} = \mathbf{f}_{sjo}(\mathbf{x}_{j(0)}, \mathbf{x}_{o(0)} + \gamma_{oml(dt)}^{(0)} \mathbf{\Delta}_{m_l(0)}) = \frac{\partial \mathbf{f}_{sjo}(\mathbf{x}_{j(0)}, \mathbf{x}_{o(0)})}{\partial \mathbf{x}_o} \gamma_{oml(dt)}^{(0)} \mathbf{\Delta}_{m_l(0)} $$

Therefore, by denoting $\mathbf{o}$ the set of points belonging to $\nu_j$ so that $m_l r_o = 1$, we obtain:

$$ \begin{aligned} \mathbf{x}_{j(2dt)} &= \mathbf{x}_{j(0)} + \frac{dt^2}{m_j} \sum_{o \in \mathbf{o}} \frac{\partial \mathbf{f}_{sjo}(\mathbf{x}_{j(dt)}, \mathbf{x}_{o(0)})}{\partial \mathbf{x}_o} \gamma_{oml(dt)}^{(0)} \mathbf{\Delta}_{m_l(0)} \\ &= \mathbf{x}_{j(0)} + \gamma_{jml(2dt)}^{(0)} \mathbf{\Delta}_{m_l(0)} \end{aligned} \quad (3.44) $$

Finally, for all other points $P_k$ in the mesh so that $m_l r_k > 2$, we have of course $\mathbf{x}_{k(2dt)} = \mathbf{x}_{k(0)}$.

For any $t = kdt$: From (3.33), (3.35) and (3.44), and by analogy with the 1D case,
we can show that for any model point $P_n$ in the 3D mesh:

$$ \mathbf{x}_{n(t)} = \hat{\mathbf{x}}_{n(t)} + \gamma_{nm_l}^{(k-m_l r_n)} \mathbf{\Delta}_{m_l(t-m_l r_n dt)} \quad (3.45) $$

The intermediate equations that leads to the derivation of (3.45) are provided at the end of this paragraph, where $\hat{\mathbf{x}}_{n(t)}$ denotes the position of $P_n$ without considering the last effective displacement $\mathbf{\Delta}_{m_l(t-m_l r_n dt)}$, which is given by

$$ \hat{\mathbf{x}}_{n(t)} = \mathbf{x}_{n(0)} + \sum_{q=0}^{k-(m_l r_n + 1)} \gamma_{nm_l((k-q)dt)}^{(q)} \mathbf{\Delta}_{m_l(qdt)} \quad (3.46) $$

and the propagation coefficients are given by

$$ \gamma_{nm_l(t)}^{(q)} = \gamma_{nm_l(t-dt)}^{(q)} + \frac{dt^2}{m_n} \sum_{j \in \nu_n} \tilde{f}_{s_{nj}(k,q)} + (1 - \frac{dt}{m_n} D_v) \delta_{nm_l(t-dt)}^{(q)} \quad (3.47) $$

with

$$ \tilde{f}_{s_{nj}(k,q)} = \frac{\partial \mathbf{f}_{s_{nj}}(\tilde{\mathbf{x}}_n((m_l r_n + q - 1)dt), \tilde{\mathbf{x}}_j((m_l r_n + q - 1)dt))}{\partial \mathbf{x}_j} (\gamma_{j m_l ((k-1-q)dt)}^{(q)} - \gamma_{n m_l ((k-1-q)dt)}^{(q)}) \quad (3.48) $$

$$ \delta_{nm_l(t-dt)}^{(q)} = \gamma_{nm_l(t-dt)}^{(q)} - \gamma_{nm_l(t-2dt)}^{(q)} \quad (3.49) $$

and

1. $\tilde{\mathbf{x}}_j = \hat{\mathbf{x}}_j$ and $\tilde{\mathbf{x}}_n = \mathbf{x}_n$ if $m_l r_j < m_l r_n$

2. $\tilde{\mathbf{x}}_j = \hat{\mathbf{x}}_j$ and $\tilde{\mathbf{x}}_n = \hat{\mathbf{x}}_n$ if $m_l r_j = m_l r_n$

3. $\tilde{\mathbf{x}}_j = \mathbf{x}_j$ and $\tilde{\mathbf{x}}_n = \hat{\mathbf{x}}_n$ if $m_l r_j > m_l r_n$

4. $\gamma_{m_l m_l}^{(q)} = I_{3 \times 3}, \forall q$ and $\forall t$

We recall that the index $q$ is related to the displacement applied on $P_{m_l}$ at time $qdt$. It has to be noted that if $t \le m_l r_n dt$ and $t \le m_l r_j dt$, $\tilde{f}_{s_{nj}(k,q)} = 0$ because $\mathbf{f}_{s_{nj}(t)} = 0$. Indeed, we remind that the effect of $\mathbf{\Delta}_{m_l(t)}$ applied on $P_{m_l}$ at time $t$ will not take place instantaneously on $P_n$, but at time $m_l t_n$ that depends on the position of $P_n$ with respect to $P_{m_l}$. For all $t < m_l t_n$, we have $\gamma_{m_l m_l}^{(q)} = 0$. Finding $m_l t_n$ for a 3D mesh is easy since
the mesh can be considered as a graph, with $P_{ml}$ as a root and its neighbors as the leaves. Thus, $^m t_n$ depends on the number of layers between the root and $P_n$. If this number is $^m r_n - 1$, then $^m t_n = ^m r_n dt$, similarly as for the 1D mesh considered in the previous section.

By comparing (3.47) to (3.25) and exploiting (3.48) and (3.32), we can notice that the form of the propagation coefficient $\gamma_{ij}$ given in (3.25) for a 1D mesh can be obtained from the form (3.47) obtained for a 3D mesh by neglecting the second term of $\frac{\partial \mathbf{f}_{s_{ij}}}{\partial \mathbf{x}_j}$. This implies that the forces $\mathbf{f}_{s_{ij}}$ are linear with respect to $\mathbf{x}_j$ and $\mathbf{x}_n$, i.e., $\frac{\partial \mathbf{f}_{s_{nj}}}{\partial \mathbf{x}_j} = \mathbf{K}_{nj} \mathbf{I}_{3 \times 3}$.

### 3.1.2.2 Demonstration of (3.45), (3.46) and (3.47)

In what follows, we demonstrate by induction the recurrence relations (3.45), (3.46) and (3.47).

First, by instantiating (3.46) and (3.45) for $t = dt$ and $t = 2dt$ and for any point $P_i$ within the vicinity of $P_{ml}$, which means $^m r_i = 1$, $k = 1$ for $t = dt$, and $k = 2$ for $t = 2dt$, we have:

$$ \begin{aligned} \hat{\mathbf{x}}_{i(dt)} &= \mathbf{x}_{i(0)} \\ \mathbf{x}_{i(dt)} &= \mathbf{x}_{i(0)} + \gamma_{iml(dt)}^{(0)} \Delta_{ml(0)} \end{aligned} \qquad (3.50) $$

$$ \begin{aligned} \hat{\mathbf{x}}_{i(2dt)} &= \mathbf{x}_{i(0)} + \gamma_{iml(2dt)}^{(0)} \Delta_{ml(0)} \\ \mathbf{x}_{i(2dt)} &= \hat{\mathbf{x}}_{i(2dt)} + \gamma_{iml(dt)}^{(1)} \Delta_{ml(dt)} \\ &= \mathbf{x}_{i(0)} + \gamma_{iml(2dt)}^{(0)} \Delta_{ml(0)} + \gamma_{iml(dt)}^{(1)} \Delta_{ml(dt)} \end{aligned} \qquad (3.51) $$

where (3.50) and (3.51) respectively correspond to (3.33) and (3.40). Then, by instantiating (3.47), we obtain

$$ \gamma_{iml(dt)}^{(0)} = \frac{dt^2}{m_i} \sum_{j \in \nu_i} \frac{\partial \mathbf{f}_{s_{ij}}(\tilde{\mathbf{x}}_{i(0)}, \tilde{\mathbf{x}}_{j(0)})}{\partial \mathbf{x}_j} \gamma_{jml(0)}^{(0)} \qquad (3.52) $$

$$ \gamma_{iml(dt)}^{(1)} = \frac{dt^2}{m_i} \sum_{j \in \nu_i} \frac{\partial \mathbf{f}_{s_{ij}}(\tilde{\mathbf{x}}_{i(dt)}, \tilde{\mathbf{x}}_{j(dt)})}{\partial \mathbf{x}_j} \gamma_{jml(0)}^{(0)} \qquad (3.53) $$

$$ \begin{aligned} \gamma_{iml(2dt)}^{(0)} = \gamma_{iml(dt)}^{(0)} &+ \left(1 - \frac{dt}{m_i} D_v \right) \gamma_{iml(dt)}^{(0)} \\ &+ \frac{dt^2}{m_i} \sum_{j \in \nu_i} \frac{\partial \mathbf{f}_{s_{ij}}(\tilde{\mathbf{x}}_{i(0)}, \tilde{\mathbf{x}}_{j(0)})}{\partial \mathbf{x}_j} (\gamma_{jml(dt)}^{(0)} - \gamma_{iml(dt)}^{(0)}) \end{aligned} \qquad (3.54) $$
At time $t = dt$, the displacement exerted on $P_{m_l}$ impacts only the points $P_i$ in its neighborhood, since all other points remain at rest. Consequently, the only point in $\nu_i$ that has to be considered in $f_{s_{ij}}$ for (3.52) and (3.53) is $P_j = P_{m_l}$. We thus have $\gamma_{m_l m_l(0)}^{(0)} = I_{3 \times 3}$, $\tilde{x}_{i(0)} = x_{i(0)}$, $\tilde{x}_{m_l(0)} = x_{m_l(0)}$, $\tilde{x}_{i(dt)} = x_{i(dt)}$ and $\tilde{x}_{m_l(dt)} = x_{m_l(0)} + \Delta_{m_l(dt)}$, from which we deduce:

$$ \gamma_{im_l(dt)}^{(0)} = \frac{dt^2 \partial f_{s_{im_l}}(x_{i(0)}, x_{m_l(0)})}{m_i \partial x_{m_l}} \quad (3.55) $$

$$ \gamma_{im_l(dt)}^{(1)} = \frac{dt^2 \partial f_{s_{im_l}}(x_{i(dt)}, x_{m_l(0)} + \Delta_{m_l(dt)})}{m_i \partial x_{m_l}} \quad (3.56) $$

which respectively correspond to (3.34) and (3.41).

At time $t = 2dt$, by decomposing $\nu_i$ as $m_l$ and $\nu_i - \{m_l\}$, we obtain from (3.54):

$$ \gamma_{im_l(2dt)}^{(0)} = \gamma_{im_l(dt)}^{(0)} + \left(1 - \frac{dt}{m_i} D_v\right) \gamma_{im_l(dt)}^{(0)} + \frac{dt^2 \partial f_{s_{im_l}}(x_{i(0)}, x_{m_l(0)})}{m_i \partial x_{m_l}} (I_{3 \times 3} - \gamma_{im_l(dt)}^{(0)}) \\ + \frac{dt^2}{m_{ij \in \nu_i - m_l}} \sum_{j \in \nu_i - m_l} \frac{\partial f_{s_{ij}}(x_{i(0)}, x_{j(0)})}{\partial x_j} (\gamma_{j m_l(dt)}^{(0)} - \gamma_{i m_l(dt)}^{(0)}) $$

which corresponds to (3.42).

We have validated (3.45), (3.46) and (3.47) for the points in the vicinity of the manipulated point at $t = dt$ and $t = 2dt$. Now, we consider the points $P_{v_i}$ that are not in the neighborhood of $P_{m_l}$. The points such that $m_l r_{v_i} > 2$ are still at rest for $t = dt$ and $t = 2dt$, so we have just to concentrate on points such that $m_l r_{v_i} = 2$. Let us denote $v$ as the set of these points. By instantiating (3.46), (3.45) and (3.47) at $t = dt$ and at $t = 2dt$ for any point $P_{v_i}$ belonging to $v$, we have:

$$ \begin{aligned} \hat{\mathbf{x}}_{v_i(dt)} &= \mathbf{x}_{v_i(dt)} = \mathbf{x}_{v_i(0)} \\ \hat{\mathbf{x}}_{v_i(2dt)} &= \mathbf{x}_{v_i(0)} \\ \mathbf{x}_{v_i(2dt)} &= \hat{\mathbf{x}}_{v_i(2dt)} + \gamma_{v_i m_l(2dt)}^{(0)} \Delta_{m_l(0)} = \mathbf{x}_{v_i(0)} + \gamma_{v_i m_l(2dt)}^{(0)} \Delta_{m_l(0)} \end{aligned} \quad (3.57) $$

with

$$ \gamma_{v_i m_l(2dt)}^{(0)} = \frac{dt^2}{m_{v_i}} \sum_{j \in \nu_{v_i}} \tilde{f}_{s_{v_i j}(2,0)} $$

The only points in $\nu_{v_i}$ that need to be considered are the points $P_j$ such that $m_l r_j = 1$. We denote the set of these points as $\mathbf{o}$. For points that belong to $\nu_v$ but do not belong
to $\mathbf{o}$, the force $\mathbf{f}_s$ is indeed zero. Therefore, we obtain $\tilde{\mathbf{x}}_{o(dt)} = \hat{\mathbf{x}}_{o(dt)} = \mathbf{x}_{o(0)}$ for $\mathbf{o} \in \mathbf{o}$ and $\tilde{\mathbf{x}}_{v_i(dt)} = \mathbf{x}_{v_i(dt)}$ for $v_i \in \mathbf{v}$, which implies that:

$$ \gamma_{v_i m_l (2dt)}^{(0)} = \frac{dt^2}{m_{v_i}} \sum_{o \in \mathbf{o}} \frac{\partial \mathbf{f}_{s_{v_i o}}(\mathbf{x}_{v_i(dt)}, \mathbf{x}_{o(0)})}{\partial \mathbf{x}_o} \gamma_{o m_l (dt)}^{(0)} \quad (3.58) $$

By injecting (3.58) in (3.57), we directly obtain (3.44), as expected.

So far, we have validated (3.45), (3.46) and (3.47) for all the model points at $t = dt$ and $t = 2dt$.

Let us now assume that the position of any point $P_n$ is given by (3.45) at time $t = kdt$, the goal is to obtain the same form at time $t + dt$. As usual, from (3.3), we have:

$$ \mathbf{x}_{n(t+dt)} = \mathbf{x}_{n(t)} + \left(1 - \frac{dt}{m_n}\right) D_v(\mathbf{x}_{n(t)} - \mathbf{x}_{n(t-dt)}) + \frac{dt^2}{m_n} \sum_{j \in \nu_n} \mathbf{f}_{s_{nj}}(\mathbf{x}_{n(t)}, \mathbf{x}_{j(t)}) \quad (3.59) $$

Using successive first-order approximations of $\mathbf{f}_{s_{nj}}$, similarly as done in (3.37)-(3.39), we have to distinguish three cases:

Case 1: $m_l r_j = m_l r_n + 1$ ($P_j$ is further from $P_{m_l}$ than $P_i$):

$$ \begin{aligned} \mathbf{f}_{s_{nj}}(\mathbf{x}_{n(t)}, \mathbf{x}_{j(t)}) &= \mathbf{f}_{s_{nj}}(\hat{\mathbf{x}}_{n(t)}, \mathbf{x}_{j(t)}) + \frac{\partial \mathbf{f}_{s_{nj}}(\hat{\mathbf{x}}_{n(t)}, \mathbf{x}_{j(t)})}{\partial \mathbf{x}_j} \gamma_{n m_l (m_l r_n dt)}^{(k-m_l r_n)} \Delta_{m_l (t-m_l r_n dt)} \\ &= \mathbf{f}_{s_{nj}}(\hat{\mathbf{x}}_{n(t)}, \mathbf{x}_{j(t)}) + \tilde{\mathbf{f}}_{s_{nj}(k+1, k-m_l r_n)} \Delta_{m_l ((k-m_l r_n)dt)} \\ &= \sum_{q=0}^{k-m_l r_n} \tilde{\mathbf{f}}_{s_{nj}(k+1, q)} \Delta_{m_l (qdt)} \end{aligned} $$

which is obtained from (3.48).

Case 2: $m_l r_j = m_l r_n$ ($P_j$ and $P_i$ are at the same distance of $P_{m_l}$):

$$ \begin{aligned} \mathbf{f}_{s_{nj}}(\mathbf{x}_{n(t)}, \mathbf{x}_{j(t)}) &= \mathbf{f}_{s_{nj}}(\hat{\mathbf{x}}_{n(t)}, \hat{\mathbf{x}}_{j(t)}) \\ &\quad + \frac{\partial \mathbf{f}_{s_{nj}}(\hat{\mathbf{x}}_{n(t)}, \hat{\mathbf{x}}_{j(t)})}{\partial \mathbf{x}_j} (\gamma_{j m_l (m_l r_n dt)}^{(k-m_l r_n)} - \gamma_{n m_l (m_l r_n dt)}^{(k-m_l r_n)}) \Delta_{m_l (t-m_l r_n dt)} \\ &= \sum_{q=0}^{k-m_l r_n} \tilde{\mathbf{f}}_{s_{nj}(k+1, q)} \Delta_{m_l (qdt)} \end{aligned} $$
Case 3: $^{m_l}r_j = ^{m_l}r_n - 1$ ($P_j$ is closer from $P_{m_l}$ than $P_i$):

$$
\begin{align*}
\boldsymbol{f}_{\boldsymbol{s}nj}(\boldsymbol{x}_{n(t)}, \boldsymbol{x}_{j(t)}) &= \boldsymbol{f}_{\boldsymbol{s}nj}(\boldsymbol{x}_{n(t)}, \hat{\boldsymbol{x}}_{j(t)}) + \frac{\partial \boldsymbol{f}_{\boldsymbol{s}nj}(\boldsymbol{x}_{n(t)}, \hat{\boldsymbol{x}}_{j(t)})}{\partial \boldsymbol{x}_j} \boldsymbol{\gamma}_{j m_l}^{(k-m_l r_n+1)} \binom{m_l r_j dt}{(m_l r_j dt)} \Delta_{m_l} (t-m_l r_n dt+dt) \\
&= \sum_{q=0}^{k-m_l r_n} \tilde{\boldsymbol{f}}_{\boldsymbol{s}nj(k+1,q)} \Delta_{m_l}(qdt) + \frac{\partial \boldsymbol{f}_{\boldsymbol{s}nj}(\boldsymbol{x}_{n(t)}, \hat{\boldsymbol{x}}_{j(t)})}{\partial \boldsymbol{x}_j} \boldsymbol{\gamma}_{j m_l}^{(k-m_l r_n+1)} \binom{m_l r_j dt}{(m_l r_j dt)} \Delta_{m_l} (t-m_l r_n dt+dt)
\end{align*}
$$

By using (3.45), (3.46), (3.47) and the three previous equations in (3.59), we obtain

$$
\begin{align}
\mathbf{x}_{n(t+dt)} &= \mathbf{x}_{n(0)} + \sum_{q=0}^{k-(m_l r_n+1)} \gamma_{nm_l((k-q)dt)}^{(q)} \Delta_{m_l(qdt)} + \gamma_{nm_l(m_l r_n dt)}^{(k-m_l r_n)} \Delta_{m_l((k-m_l r_n)dt)} \tag{3.60} \\
&\quad + \frac{dt^2}{m_{nj \in \nu_n}} \delta^{[m_l r_j=m_l r_n-1]} \frac{\partial \mathbf{f}_{snj}(\mathbf{x}_{n(t)}, \hat{\mathbf{x}}_{j(t)})}{\partial \mathbf{x}_j} \gamma_{j m_l}^{(k+1-m_l r_n)} \binom{m_l r_j dt}{(m_l r_j dt)} \Delta_{m_l((k+1-m_l r_n)dt)} \nonumber \\
&\quad + (1-\frac{dt}{m_n} D_v) \sum_{q=0}^{k-m_l r_n} \delta_{nm_l(t-qdt)}^{(q)} \Delta_{m_l(qdt)} + \frac{dt^2}{m_{nj \in \nu_n}} \sum_{q=0}^{k-m_l r_n} \tilde{\mathbf{f}}_{snj(k+1,q)} \Delta_{m_l((q)dt)} \nonumber \\
&= \mathbf{x}_{n(0)} + \left[ \gamma_{nm_l(kdt)}^{(0)} + (1-\frac{dt}{m_n} D_v) \delta_{nm_l(kdt)}^{(0)} + \frac{dt^2}{m_{nj \in \nu_n}} \sum_{j \in \nu_i} \tilde{\mathbf{f}}_{snj(k+1,0)} \right] \Delta_{m_l(0)} + \dots \nonumber \\
&\quad + \left[ \gamma_{nm_l(m_l r_n dt)}^{(k-m_l r_n)} + (1-\frac{dt}{m_n} D_v) \delta_{nm_l(m_l r_n dt)}^{(k-m_l r_n)} + \frac{dt^2}{m_n} \sum_{j \in \nu_n} \tilde{\mathbf{f}}_{snj(k+1,k-m_l r_n)} \right] \Delta_{m_l((k-m_l r_n)dt)} \nonumber \\
&\quad + \frac{dt^2}{m_n} \sum_{j \in \nu_n} \left[ \delta^{[m_l r_j=m_l r_n-1]} \frac{\partial \mathbf{f}_{snj}(\mathbf{x}_{n(t)}, \hat{\mathbf{x}}_{j(t)})}{\partial \mathbf{x}_j} \gamma_{j m_l}^{(k+1-m_l r_n)} \binom{m_l r_j dt}{(m_l r_j dt)} \Delta_{m_l((k+1-m_l r_n)dt)} \right] \nonumber
\end{align}
$$

where $\delta_{[c]} = 1$ if the condition c is true and $\delta_{[c]} = 0$ otherwise.

Recalling that $t = kdt$ and rearranging the terms in (3.60), we obtain as expected:

$$
\mathbf{x}_{n(t+dt)} = \hat{\mathbf{x}}_{n(t+dt)} + \gamma_{nm_l}^{(k+1-m_l r_n)} \binom{m_l r_n}{(m_l r_n dt)} \Delta_{m_l(t+dt-m_l r_n dt)}
$$

with, by identifying the terms in (3.60), we get:

$$
\hat{\mathbf{x}}_{n(t+dt)} = \mathbf{x}_{n(0)} + \sum_{q=0}^{k+1-(m_l r_n+1)} \gamma_{nm_l((k-q)dt)}^{(q)} \Delta_{m_l(qdt)}
$$

and

$$
\gamma_{nm_l(t+dt)}^{(q)} = \gamma_{nm_l(t)}^{(q)} + \frac{dt^2}{m_n} \sum_{j \in \nu_n} \tilde{f}_{snj(k+1,q)} + \left(1 - \frac{dt}{m_n} D_v\right) \delta_{nm_l(t)}^{(q)} \quad (3.61)
$$
Finally, by identifying $\gamma_{nm_l}^{(k+1-m_l r_n)}$ as the term associated to $\Delta_{m_l((k+1-m_l r_n)dt)}$ in (3.60), we see that

$$ \gamma_{nm_l}^{(k+1-m_l r_n)} \Big|_{(m_l r_n dt)} = \frac{dt^2}{m_n} \sum_{j \in \nu_n} \left[ \delta_{[m_l r_j = m_l r_n - 1]} \frac{\partial \boldsymbol{f}_{sj}(\boldsymbol{x}_n(t), \hat{\boldsymbol{x}}_{j(t)})}{\partial \boldsymbol{x}_j} \gamma_{jml}^{(k+1-m_l r_n)} \Big|_{(m_l r_j dt)} \right], \quad (3.62) $$

which can be written under the form given by (3.61).

Indeed, setting $q = k + 1 - m_l r_n$ and $t = (m_l r_n - 1)dt$ in (3.61), and using (3.48), we obtain from (3.61):

$$ \gamma_{nm_l}^{(k+1-m_l r_n)} \Big|_{(m_l r_n dt)} = \frac{dt^2}{m_n} \sum_{j \in \nu_n} \left[ \frac{\partial \boldsymbol{f}_{sj}(\tilde{\boldsymbol{x}}_n(kdt), \tilde{\boldsymbol{x}}_{j(kdt)})}{\partial \boldsymbol{x}_j} \left( \gamma_{jml}^{((k+1-m_l r_n))} \Big|_{((m_l r_n - 1)dt)} - \gamma_{jml}^{((k+1-m_l r_n))} \Big|_{((m_l r_n - 1)dt)} \right) \right] \quad (3.63) $$

Note that $\gamma_{nm_l}^{(k+1-m_l r_n)} \Big|_{((m_l r_n - 1)dt)} = \delta_{nm_l}^{(k+1-m_l r_n)} \Big|_{((m_l r_n - 1)dt)} = 0$ since $\gamma_{nm_l}^{(k+1-m_l r_n)} \Big|_{(t)} = 0$, $\forall t < m_l r_n dt$. Moreover, we can notice that only $\gamma_{jml}^{((k+1-m_l r_n))} \Big|_{((m_l r_n - 1)dt)}$ corresponding to $P_j$ with $j \in \nu_n$ such that $m_l r_j = (m_l r_n - 1)$ are not equal to zero, which corresponds to Case 3 above. As a result, we have $\tilde{\boldsymbol{x}}_n(kdt) = \boldsymbol{x}_n(kdt)$, $\tilde{\boldsymbol{x}}_{j(kdt)} = \hat{\boldsymbol{x}}_{j(kdt)}$. Since $\gamma_{jml}^{((k+1-m_l r_n))} \Big|_{((m_l r_n - 1)dt)} = 0$, this allows (3.63) to be presented in the form given in (3.62), which ends our demonstration.

### 3.1.2.3 Multiple manipulated points

What we have developed so far focused on a single manipulated point. As we showed previously in the 1D scenario and especially in (3.26), if M manipulated points are used to deform the object simultaneously, the motions of the model points can be determined by summing up the effects of the manipulated points using the propagation coefficients. Therefore, we can show by similarity with the 1D case that:

$$ \boldsymbol{x}_n(t+dt) = \boldsymbol{x}_n(0) + \sum_{l=1}^{M} \gamma_{nm_l}^{(k-m_l r_n)} \Big|_{(m_l r_n dt)} \boldsymbol{\Delta}_{m_l(t-m_l r_n dt)} + \sum_{l=1}^{M} \sum_{q=0}^{k-(m_l r_n+1)} \gamma_{nm_l}^{(q)} \Big|_{((k-q)dt)} \boldsymbol{\Delta}_{m_l(qdt)} \quad (3.64) $$

and

$$ \dot{\boldsymbol{x}}_n(t+dt) = \sum_{l=1}^{M} \sum_{q=0}^{k-m_l r_n} \boldsymbol{\delta}_{nm_l(t-qdt)}^{(q)} \dot{\boldsymbol{x}}_{ml}(qdt) \quad (3.65) $$

where $\gamma_{nm_l}^{(q)}$ and $\delta_{nm_l}^{(q)}$ are respectively given by (3.47) and (3.49).
### 3.1.3 Control scheme

In this section, $Q$ feature points are considered, and the goal is to bring them to some desired positions by acting on $M$ manipulated points through a kinematic controller whose outputs are the velocities of these manipulated points. We denote $[P_{f_1}, P_{f_2}, \dots, P_{f_Q}]$ the feature points and $[P_{m_1}, P_{m_2}, \dots, P_{m_M}]$ the manipulated points.

The measurements from a RGB-D camera are provided with a quite low frequency (typically 33 Hz). On the other hand, the model is updated with a higher frequency than the control law to better represent its dynamic behavior. Denoting $\rho$ the ratio between these two frequencies (in practice, we used $dt = 1.5$ ms so that $\rho = 20$), the control law is updated every $\rho dt$, which means that the velocity of the manipulated points remains constant between $h\rho dt$ and $(h+1)\rho dt$, $\forall h \in \mathbb{N}$. Thanks to (3.65), it is possible to determine the link between the velocity of the features points at time $(h+1)\rho dt$ and the velocity of the manipulated points at time $h\rho dt$, from which the control law will be easily designed.

We start again by considering a single manipulated point, $P_{m_l}$. From (3.65), we have for any feature point $P_{f_i}$:

$$ \dot{\mathbf{x}}_{f_i((h+1)\rho dt)} = \sum_{q=0}^{((h+1)\rho - m_l r_{f_i} - 1)} \delta_{f_i m_l (t-qdt)}^{(q)} \dot{\mathbf{x}}_{m_l (qdt)} \quad (3.66) $$

Furthermore, we know that $\dot{\mathbf{x}}_{m_l(t)} = \dot{\mathbf{x}}_{m_l(h\rho dt)}$, $\forall t \in [h\rho dt; (h+1)\rho dt]$, then (3.66) can be rewritten as follows:

$$ \dot{\mathbf{x}}_{f_i((h+1)\rho dt)} = \sum_{q=h\rho}^{((h+1)\rho - m_l r_{f_i} - 1)} \delta_{f_i m_l (t-qdt)}^{(q)} \dot{\mathbf{x}}_{m_l (h\rho dt)} + \sum_{q=0}^{(h\rho - 1)} \delta_{f_i m_l (t-qdt)}^{(q)} \dot{\mathbf{x}}_{m_l (qdt)} $$

which leads to the very simple form:

$$ \dot{\mathbf{x}}_{f_i((h+1)\rho dt)} = \mathbf{A}_{f_i m_l} \dot{\mathbf{x}}_{m_l(h\rho dt)} + \mathbf{b}_{f_i m_l} $$

with:

$$ \mathbf{A}_{f_i m_l} = \sum_{q=h\rho}^{((h+1)\rho - m_l r_{f_i} - 1)} \delta_{f_i m_l (t-qdt)}^{(q)} \quad (3.67) $$

which combines the effects of the most recent velocity $\dot{\mathbf{x}}_{m_l(h\rho dt)}$ applied on the manipulated
point $P_{ml}$, and

$$ \mathbf{b}_{f_i m_l} = \sum_{q=0}^{(h\rho-1)} \delta_{f_i m_l (t-qdt)}^{(q)} \dot{\mathbf{x}}_{m_l (qdt)} \quad (3.68) $$

which represents the effects of the velocities previously applied on $P_{ml}$.

By now considering $M$ manipulated points and $Q$ feature points, we finally obtain:

$$ \dot{\mathbf{x}}_f = \mathbf{A} \dot{\mathbf{x}}_m + \mathbf{b}, \quad (3.69) $$

with $\dot{\mathbf{x}}_f = (\dot{\mathbf{x}}_{f_1}, \dot{\mathbf{x}}_{f_2}, \dots, \dot{\mathbf{x}}_{f_Q})$, $\dot{\mathbf{x}}_m = (\dot{\mathbf{x}}_{m_1}, \dot{\mathbf{x}}_{m_2}, \dots, \dot{\mathbf{x}}_{m_M})$, and where the coefficients of matrix $\mathbf{A}$ and vector $\mathbf{b}$ are given by:

$$ \mathbf{A} = \begin{bmatrix} \mathbf{A}_{f_1 m_1} & \mathbf{A}_{f_1 m_2} & \dots & \mathbf{A}_{f_1 m_M} \\ \mathbf{A}_{f_2 m_1} & \mathbf{A}_{f_2 m_2} & \dots & \mathbf{A}_{f_2 m_M} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{A}_{f_Q m_1} & \mathbf{A}_{f_Q m_2} & \dots & \mathbf{A}_{f_Q m_M} \end{bmatrix}, \mathbf{B} = \begin{bmatrix} \sum_{l=1}^{M} \mathbf{b}_{f_1 m_l} \\ \sum_{l=1}^{M} \mathbf{b}_{f_2 m_l} \\ \vdots \\ \sum_{l=1}^{M} \mathbf{b}_{f_Q m_l} \end{bmatrix} $$

Equation (3.69) presents explicitly the velocity of all feature points in function of the motions of the manipulated points deforming the object. $\mathbf{A}$ can be considered as the interaction matrix or the deformation Jacobian by similarity with the classical visual servoing framework [Chaumette & Hutchinson, 2008; Lagneau et al., 2020a; Navarro-Alarcon et al., 2013; J. Zhu et al., 2018], while $\mathbf{b}$ is a feed-forward term that does not appear in classical shape servoing work. We will see the importance of this term in Sections 3.2.5 and 4.2.2.2, which include simulation and experimental results.

Then, to indirectly position these feature points to their desired positions $\mathbf{x}_f^* = (\mathbf{x}_{f_1}^*, \mathbf{x}_{f_2}^*, \dots, \mathbf{x}_{f_Q}^*)$ with an exponential decrease so that $\dot{\mathbf{x}}_f = -\lambda(\mathbf{x}_f - \mathbf{x}_f^*)$ where $\lambda > 0$, we deduce from (3.69) the closed-loop control law

$$ \dot{\mathbf{x}}_m = -\lambda \hat{\mathbf{A}}^+ (\mathbf{x}_f - \mathbf{x}_f^*) - \hat{\mathbf{A}}^+ \hat{\mathbf{b}} \quad (3.70) $$

where $\hat{\mathbf{A}}$ is the approximation of $\mathbf{A}$, and $\hat{\mathbf{A}}^+$ denotes the Moore-Penrose pseudo-inverse of $\hat{\mathbf{A}}$. Regarding $\hat{\mathbf{b}}$, it is the approximation of $\mathbf{b}$. In (3.70), we used $\hat{\mathbf{A}}^+$ and $\hat{\mathbf{b}}$ instead of $\mathbf{A}^+$ and $\mathbf{b}$ because we consider that the modeling proposed in (3.69) is accurate but requires exact knowledge of the model parameters. However, obtaining the exact parameter values can be challenging, leading us to use model coarse parameters. Since $\mathbf{A}$ and $\mathbf{b}$ are dependent on these model parameters, we opted for approximations, which explains our choice of $\hat{\mathbf{A}}$
and $\hat{b}$.

### 3.1.3.1 Stability analysis

In this paragraph, we explore the stability analysis and convergence of the closed-loop control law proposed in (3.70) through the application of Lyapunov analysis.

In practice, in order to bring the feature points to their desired positions, it is well-known that at least one manipulated point per feature point has to be involved [Wada et al., 2001b], which implies $M \ge Q$. However, the convergence of the system does not solely depend on the aforementioned condition, but also depends on the rank of $A$, which is influenced by the relative initial positions between the manipulated points and the feature points.

Let define a Lyapunov function $\mathcal{L} = 1/2\|x_f - x_f^*\|^2$, and by considering that the desired points are fixed and replacing (3.70) in (3.69), then we obtain:

$$ \dot{\mathcal{L}} = -\lambda(x_f - x_f^*)^T A \hat{A}^+ (x_f - x_f^*) + (x_f - x_f^*)^T (b - A \hat{A}^+ \hat{b}) \quad (3.71) $$

By denoting the second term as $\epsilon$, where $\epsilon = (b - A\hat{A}^+\hat{b})$ and assuming it remains bounded during the deformation. Then, a sufficient condition to ensure the global asymptotic stability of the system is given by:

$$ A\hat{A}^+ > 0 \quad (3.72) $$

If both matrices $A$ and $\hat{A}^+$ are of full rank $3Q$, and $\hat{A}^+$ is not too coarse, thus condition (3.72) is ensured. Note that when $A\hat{A}^+ = I_{3Q}$, which needs $\hat{A}$ to be perfectly estimated, the specified pure exponential decoupled decrease is expected if the MSM fits with the true behavior of the object. Otherwise, perturbations in this decrease will occur.

In case the system is underactuated, i.e. $M < Q$, it is possible to demonstrate that the system is locally asymptotically stable if the condition $\hat{A}^+ A > 0$ is ensured, which occurs in case $A$ and $\hat{A}^+$ are of full rank $3M$ and, once again, $\hat{A}^+$ is not too coarse. This means that the system will converge if the error between $x_f$ and $x_f^*$ is initially small [Chaumette & Hutchinson, 2008].

Furthermore, when the system is stable, then the positioning error converges to a value $\bar{e}$, with $\bar{e} = \epsilon^*/\lambda$. However, due to the system stability, it will converge to a quasi-static state where $\bar{e} \to 0$ because $b \to 0$ and $\hat{b} \to 0$. This occurs since $b$ and $\hat{b}$ encapsulate the effects of previous velocities applied to the manipulated points, which are absorbed by the MSM when it reaches a quasi-static state.
Regarding the rank of **A**, it can be analyzed by examining its terms using (3.67). Indeed, **A**<sub>fiml</sub> ≠ **0** when ρ ≥ <sup>m<sub>l</sub></sup>/<sub>r<sub>f<sub>i</sub></sub> + 1, while this condition is not ensured when a manip- ulated point is located far away from a feature point (ρ < <sup>m<sub>l</sub></sup>/<sub>r<sub>f<sub>i</sub></sub> + 1). Increasing ρ would be a solution, but a compromise has to be found to maintain a real-time performance. Furthermore, if all the manipulated points are located close to a particular feature point P<sub>f<sub>i</sub></sub> and far away from the other feature points, then the row in **A** corresponding to P<sub>f<sub>i</sub></sub> will be of full rank 3 while all other rows in **A** will only contain zero values, leading to rank(**A**) = 3.

As a conclusion, when rank(**A**) < 3Q and the initial error **x**<sub>f</sub> − **x**<sub>f</sub>\* is large, the system
may reach a local minimum for which the error **x**<sub>f</sub> − **x**<sub>f</sub>\* ≠ **0** at the end of the deformation
task. Looking at the rank of **A** is thus a precious indicator for the feasibility of indirectly
positioning feature points.

**3.2 Simulation results**

The main objective of the proposed controller (3.70) is to minimize the errors between
the current and desired positions of the feature points.

In this section, we show some simulation results using two types of objects, a 3D soft parallelepipedal object with either 1542 or 7710 nodes, each of dimension 12.5x12.5x6.25 cm (see Figure 3.6), and a 2D circular object of radius 10 cm with 256 nodes (see Figure 3.10).

The aim of these simulations is to demonstrate the validity and the robustness of the
proposed closed-loop system. Additionally, this section includes comparisons between the
proposed approach (PA) and two other model-free methods [Navarro-Alarcon et al., 2013;
J. Zhu et al., 2018].

To simulate the behavior of the objects, we use *dt* = 1.5 ms with a control scheme rate of 33 Hz, which corresponds to a Realsense RGB-D camera frequency (*ρ* = 20). For all results shown below, the same controller gain *λ* = 0.9 has been used. Simulation results are accessible via this online video: [Indirect-positioning-video].

**3.2.1 Object parameters**

**3.2.1.1 3D soft parallelepipedal object**

Two models with the same geometry but with a different resolution are used to simulate
the same 3D object (see Figure 3.6). Regarding the parameters of these models, they are