Keywords: Medical robotics, ultrasound elastography, visual servoing, haptic feedback

Optical flow RF Radio frequency ROI Region of interest SAD Sum of absolute difference SSD Sum of squared differences US

to always improve the quality of my work. I also appreciate his charisma that was never missing.

My most sincere thanks to Jonathan Vappou and Pierre Vieyres for the valuable review and comments to improve my manuscript. I would like to thank the members of my thesis committee Marie-Aude Vitrani, Pierre Janin and Eric Marchand for their insightful comments and encouragement during my defense. Special thanks to François Chaumette for accepting me in Lagadic team (now Rainbow team) and for his support during my doctorate. I also would like to express my gratitude to Paolo Robuffo Giordano and wish him the best as the head of Rainbow team. I am also very thankful to Fabien for his support and wise advices with the robotic platform.

I also appreciate the support of Hélène with all the administrative processes during my doctorate.

I warmly thank my fellow labmates in Lagadic/Ranbow team. I am grateful to Hadrien and Souriya who helped me to improve the writing of the abstract in french of this manuscript. Through these three years and several months in Rennes I had the opportunity to meet excellent mates who shared with me invaluable moments. I will always remember the past fellows: Bertrand, Mani, Le, Riccardo, Giovanni, Pierre, Lucas, Suman, Aly San and Vishnu. I want to thank you to Fabrizio, Usman, Agniva, Ide-Flore, Firas, Lesley, Jason, Marco Cognetti, Marco Agravi, Bryan, Rahf, Marc, Nöel and Quentin for all the good moments. This work has been inspired and supported by many people in my life, some of them are physically far but always in my mind. My sincere gratitude to Juan Gabriel Aviña who always encouraged me to pursue my doctorate. I am also grateful with Adan Salazar and Horacio Rostro for their support in my quest of a doctorate. Agradezco a toda mi familia i en México, especialmente a mis padres Isidro and Ruth por su constante apoyo durante toda mi vida. A mis suegros Refugio y Bertha por el apoyo que siempre me han brindado.

A mi hermano José Alberto un caluroso agradecimiento porque siempre ésta para mi en todo momento y porque su visita me recargó de energía. A mi hermana y sobrinas por su apoyo durante estos años lejos de casa. Agradecimientos especiales para mis cuñados y amigos María del Rayo, Mauricio y Felipe. Desde que los conozco, sus comentarios, amistad y cariño siempre fueron y han sido una gran fuente de energía en mi vida. A todo el resto de mi familia por todo el apoyo que me han brindado.

ii RÉSUMÉ EN FRANÇAIS La robotique médicale est apparue dans les années 1980 avec pour objectif de fournir aux médecins de nouveaux outils facilitant le traitement des patients. Un premier dispositif médical robotisé fut présenté en 1983. Il s'agit d'Arthrobot, un robot dédié à la chirurgie orthopédique que le chirurgien utilisait principalement pour des tâches fatigantes (par exemple, maintenir une partie du corps du patient à la même position pendant une longue période) et qu'il pouvait contrôler par commande vocale. Ce fut une innovation majeure dans le domaine des interventions chirurgicales assistées par la robotique. Suite au succès d'Arthrobot, d'autres applications analogues virent le jour tirant partie de la précision et de la dextérité des dispositifs robotisés. En 1985, le robot industriel PUMA 560 fut utilisé comme un outil de positionnement d'aiguille, appliqué à la biopsie du cerveau sous imagerie tomodensitométrique (TDM) [START_REF] Kwoh | A robot with improved absolute positioning accuracy for ct guided stereotactic brain surgery[END_REF]. En 1992, Think Surgical Inc introduisit ROBODOC, un assistant robotique pour l'arthroplastie de la hanche [START_REF] Taylor | An image-directed robotic system for precise orthopaedic surgery[END_REF]. A l'aide d'images TDM et durant l'intervention, le robot repérait la position de trois broches insérées dans la hanche du patient par le chirurgien, puis assurait la planification et la bonne exécution de tâches élémentaires. Ce robot fut le premier de son genre à être utilisé sur des patients. Voyant le guidage par imagerie comme un allié majeur de la robotique médicale, chercheurs et industriels proposèrent de nouvelles solutions reposant sur davantage de modalités d'imagerie médicale.

Avec le temps, l'imagerie médicale devint incontournable dans la plupart des interventions chirurgicales, celle-ci permettant au médecin de visualiser l'intérieur du corps humain sans avoir recours à une incision (imagerie non-invasive). Encore aujourd'hui, les modalités les plus notables sont la tomodensitométrie (TDM), l'imagerie par résonance magnétique (IRM) et l'imagerie échographique qui utilise les ultrasons (US). L'IRM est une technologie ayant recours à des champs magnétiques puissants pour produire une image des organes. L'acquisition d'une telle image est longue et nécessite des équipements encombrants, rendant l'IRM inadaptée pour la commande par vision des robots médicaux en temps réel. Malgré tout, des dispositifs robotisés ont été conçus iii pour être utilisés simultanément avec un scanner IRM, par exemple pour faciliter la manipulation d'aiguilles de biopsie [START_REF] Su | Real-time mri-guided needle placement robot with integrated fiber optic force sensing[END_REF]. L'imagerie échographie, quant à elle, s'appuie sur la propagation du son dans le corps humain pour construire des images des organes. Un de ses avantages majeurs par rapport aux autres modalités est son encombrement réduit, puisqu'elle nécessite uniquement l'utilisation d'une station (pouvant être portable) et d'une sonde échographique. Elle permet également d'acquérir des images en temps réel. Pour ces deux raisons, l'imagerie par ultrasons est la modalité la plus utilisée pour le guidage de robots médicaux. Elle fut pour la première fois utilisée à cette fin en 1999 [START_REF] Salcudean | A robot system for medical ultrasound[END_REF].

Afin de faciliter l'examen de l'artère carotide, un système robotique fut conçu pour déplacer automatiquement une sonde à ultrasons dans le but de maintenir la section de l'artère dans le plan de coupe de la sonde. Ce fut une innovation majeure dans le contrôle robotique par la vision (asservissement visuel). Dès lors, l'asservissement visuel par imagerie échographique suscita un vif intérêt dans le monde scientifique [START_REF] Azizian | Visual servoing in medical robotics: a survey. part II: tomographic imaging modalities -techniques and applications[END_REF].

Depuis, de nouvelles modalités d'imagerie ont vu le jour, issues de l'étude des informations obtenues à partir des modalités d'imagerie désormais classiques (US, IRM, TDM,. . . ). Parmi elles, l'élastographie permet d'acquérir des informations sur la rigidité d'un tissu, offrant alors à l'utilisateur de nouvelles données de diagnostic, tout en restant non-invasive. Ce concept fut approfondi durant les trois dernières décennies, en particulier pour la détection de tumeurs du sein [START_REF] Goddi | Breast elastography: A literature review[END_REF], de fibroses hépatiques à différents stades [START_REF] Barr | Elastography assessment of liver fibrosis: Society of radiologists in ultrasound consensus conference statement[END_REF] ou de cancer de la prostate [START_REF] Correas | Ultrasound elastography of the prostate: State of the art[END_REF]. Cependant, à ce jour la procédure permettant d'obtenir des images d'élastographie est réalisée manuellement. Elle nécessite une formation spécifique ainsi qu'une certaine expérience du praticien. Afin d'aider ce dernier, la robotique pourrait être utilisée à des fins d'assistance.

Motivations

La palpation manuelle est une technique médicale de diagnostic pratiquée depuis des siècles. Lorsqu'il a recours à cette approche, le praticien utilise le toucher pour évaluer la raideur des tissus du patient. Des changements de raideur peuvent alors être interprétés comme des signes d'une éventuelle maladie. Dans la pratique, en plus d'être non-invasive et simple à appliquer, cette méthode ne nécessite aucun équipement. En revanche, elle ne fournit que des informations qualitatives au praticien et nécessite une grande rigueur d'exécution. En effet, la raideur des tissus environnants peut influer sur le résultat et surtout le médecin n'a accès qu'aux tissus à portée de ses mains. Malgré cela, la palpation reste utilisée pour le diagnostic de nombreuses maladies. Par exemple, certaines maladies du foie, telles que la cirrhose ou l'hépatite, peuvent être détectées en repérant une fibrose. Certains types de cancers, tels que le cancer du sein, de la prostate iv ou de la thyroïde, peuvent être identifiés par la palpation. La détection précoce du cancer est cruciale pour accroître les chances de succès du traitement. Cependant, les techniques les plus couramment utilisées pour la détection du cancer, comme la biopsie ou la mastographie, peuvent perturber le confort du patient, le conduisant bien souvent à retarder les examens nécessaires. Par conséquent, il est primordial de proposer des outils de diagnostic indolores et non-invasifs.

L'élastographie ultrasonore est une méthode d'imagerie tactile utilisant les ondes ultrasonores pour mesurer la raideur des tissus. Cette approche permet de palier les défauts de la palpation manuelle, puisqu'elle fournit au praticien des informations importantes et précises. L'élastographie est le plus souvent manuelle et implique d'imposer un mouvement répétitif au tissu. Pour un tissu donné, la carte de raideur peut varier si la pression appliquée à ce dernier n'est pas régulière. Ainsi, il sera difficile, voire impossible, de reproduire des résultats d'élastographie d'un examen à un autre, en particulier s'il ne s'agit pas du même médecin ou si plusieurs interventions successives sont réalisées sur le même tissu. Ces inconvénients mettent en évidence la nécessité d'une méthode capable de fournir des informations plus fiables. Dans cette optique, la robotique est une alternative intéressante à l'intervention entièrement manuelle, de par la capacité des robots à effectuer des tâches répétitives avec un degré de précision constant et élevé. Cela permettrait, entre autres, de pouvoir produire les mêmes cartes d'élasticité d'un examen à un autre. En conclusion, la mise en oeuvre d'un système robotisé capable de mesurer la raideur d'un tissu pourrait apporter une aide considérable au diagnostic par palpation manuelle, conduisant par la suite à une meilleure prise en charge du patient. 1. Développement d'un système robotisé pour réaliser le mouvement de palpation des tissus de manière répétitive. Pour cela, la conception d'un système de commande robotique utilisant un transducteur à ultrasons conventionnel pour appliquer le mouvement est requise.

Objectifs de la thèse

2. Estimation quantitative de l'élasticité du tissu en temps réel par l'intermédiaire d'un processus élastographique rapide. [START_REF] Azizian | Visual servoing in medical robotics: a survey. part II: tomographic imaging modalities -techniques and applications[END_REF]. Utilisation de l'information précédente pour commander le robot afin d'assister le praticien durant l'intervention.

4. Télé-opération de la sonde à ultrasons pour explorer le tissu d'intérêt à distance.

Contributions

Durant cette thèse, plusieurs contributions ont été proposées dans le domaine de l'assistance médicale robotisée. Elles sont présentées ci-dessous :

• Une nouvelle méthodologie pour utiliser la carte de déformation du tissu comme entrée d'un schéma de contrôle par asservissement visuel. Cette contribution inclut l'estimation en temps réel de la carte d'élasticité des tissus et l'extraction des caractéristiques visuelles requise pour la commande d'un robot porteur d'une sonde ultrasonore. Cette méthodologie est dérivée pour les cas d'échographie 2D et 3D.

• Un système complet de palpation automatique fournissant une élastographie ultrasonore quantitative. Le schéma de contrôle, proposé dans le présent manuscrit, est composé de trois tâches robotiques hiérarchiques collaborant les unes avec les autres. Elles ont pour objectif de constamment afficher la carte d'élasticité d'une région d'intérêt donnée. Dans l'ordre des priorités les tâches proposeés sont : l'application d'un mouvement de palpation par une commande par retour d'effort, le centrage automatique d'une cible correspondant à un tissu rigide dans l'image par asservissement visuel, l'orientation automatique de la sonde ultrasonore afin d'observer la cible avec différents angles de vue. Cette contribution a été publiée dans deux articles de deux conférences internationales en robotique vi de premier rang : IROS 2016 [START_REF] Patlan-Rosales | Automatic palpation for quantitative ultrasound elastography by visual servoing and force control[END_REF] (IEEE / RSJ International Conference on Intelligent Robots and Systems) et ICRA 2017 [START_REF] Patlan-Rosales | A robotic control framework for 3-D quantitative ultrasound elastography[END_REF] (IEEE International Conference on Robotics and Automation).

• Une nouvelle méthode d'estimation de la carte des contraintes tissulaires. Cette approche est basée sur l'estimation des déplacements d'une région d'intérêt, au sein d'une image échographique. Ceux-ci sont calculés à l'aide des paramètres géométriques caractérisant un système de recalage d'image déformable.

• Une nouvelle approche d'estimation par échographie 2D de la carte des déformations d'une région d'intérêt donnée soumise à des perturbations de type mouvements physiologiques. La compensation de mouvement repose sur un asservissement visuel dense qui déplace la sonde de façon à annuler le mouvement relatif entre la sonde et le tissu mobile. Ainsi, une estimation robuste de l'élasticité des tissus peut être réalisée sur des structures en déplacement. Cette contribution a été publiée dans un article de la conférence internationale IROS 2017 [START_REF] Patlan-Rosales | Strain estimation of moving tissue based on automatic motion compensation by ultrasound visual servoing[END_REF] (IEEE / RSJ International Conference on Intelligent Robots and Systems).

• Un système haptique basé sur l'élastogramme des tissus. Ce dernier retourne à l'utilisateur la sensation de l'élasticité des tissus, tout en permettant à l'utilisateur de déplacer la position de la région d'intérêt à analyser dans l'image ultrasonore.

• Un dispositif associant la télé-opération de la sonde échographique au système haptique présenté précédemment.

Structure de la thèse

Le présent manuscrit de thèse est organisé comme suit :

Le chapitre 1 introduit les concepts de base et l'état de l'art de l'élastographie ultrasonore. Il est divisé en trois sections. La première décrit succinctement la théorie physique de la formation du faisceau d'ultrasons, puis explique la reconstruction géométrique de l'image ultrasonore. La deuxième présente les concepts élémentaires de l'élastographie, ainsi qu'un état de l'art. Ce dernier détaille plus particulièrement les contributions existantes dans le domaine de l'élastographie ultrasonore. La dernière section rappelle les principes de l'asservissement visuel qui seront utilisés dans les chapitres suivants.

Le chapitre 2 présente le système de palpation robotique que nous proposons pour réaliser une élastographie ultrasonore quantitative. Une partie traite du protocole expérimental envisagé dans le cadre de la thèse. Elle détaille notamment l'ensemble des vii équipements utilisés, mais aussi les outils de travail nécessaires à la mise en oeuvre de ce dispositif robotique. Ensuite, trois sections principales décrivent les tâches robotiques proposées pour concevoir le système de palpation robotique dédié à l'élastographie ultrasonore. Ces tâches sont au nombre de trois. La première correspond à un contrôle du mouvement d'oscillation de la sonde ultrasonore par une commande en effort. La deuxième consiste à extraire et utiliser les paramètres d'élasticité du tissu dans une commande par asservissement visuel permettant de maintenir une visibilité optimale du tissu analysé. Quant à la troisième, elle a pour objectif d'orienter automatiquement la sonde à ultrasons selon un angle fourni par l'utilisateur. Etant donné que ces trois tâches sont couplées, une approche hiérarchique est présentée afin de pouvoir les utiliser en combinaison dans le système de palpation proposé. Les résultats expérimentaux obtenus avec une sonde échographique 2D ou 3D sont décrits à la fin du chapitre.

Le chapitre 3 présente une méthode de compensation de mouvement robuste, utilisée pour estimer la carte des déformations d'un tissu en déplacement. Ce système complète et améliore celui décrit au chapitre 2. Ce chapitre comprend cinq sections.

La première décrit des travaux liés au suivi visuel de tissus déformables. La deuxième présente un aperçu des modèles de suivi visuel, testés avec des images échographiques réelles. Ensuite, une présentation détaillée de notre dispositif de suivi visuel dense est effectuée. Cette section traite également de l'utilisation du modèle de suivi pour estimer la carte des déformations du tissu en mouvement. La section 4 présente le schéma de contrôle utilisé pour compenser une perturbation, de type mouvement physiologique, par l'intermédiaire d'une sonde échographique 2D actionnée par un robot. Dans cette partie, la compensation du mouvement est basée sur un asservissement visuel dense. La dernière section introduit puis discute les résultats expérimentaux obtenus durant des expérimentations réalisées sur un fantôme simulant des tissus en mouvement.

Le chapitre 4 propose un système haptique basé sur la carte des contraintes du tissu, et restituant à l'utilisateur la sensation de l'élasticité des tissus. Ce chapitre est divisé en quatre sections. La première décrit les concepts de base du retour haptique jusqu'à aboutir au type de retour haptique utilisé dans ce chapitre. La deuxième introduit le modèle permettant de transformer la carte de contraintes en retour haptique. Elle comprend la description d'un système de télé-opération chargé d'assister l'utilisateur dans l'exploration des tissus avec une sonde à ultrasons. La combinaison de la télé-opération et du contrôle de la force oscillatoire est présentée dans la même section. La quatrième section détaille le protocole expérimental ainsi que les résultats obtenus avec le dispositif haptique. Elle approfondit également l'implémentation du système haptique. Puis, une conclusion est proposée pour clore ce chapitre.

Bibliography. 139

xiv LIST OF FIGURES

INTRODUCTION

Medical robotics emerged in the 1980s with the aim to provide the physicians with new tools to extend their ability to treat patients. In 1983, arthrobot was the first robotic assistant used for orthopedic surgery. This robot was voice-commanded by the surgeon to assist in tiring tasks (e.g, holding a limb at the same position for long periods). Arthrobot was a major breakthrough at the time, offering a wide perspective of robotic assisted procedures in medicine. The advantages of accuracy and dexterity of a robot brought more applications in medical procedures since the success of arthrobot. In 1985 the industrial robot PUMA 560 was used as a positioning device to orient a needle for biopsy of the brain [START_REF] Kwoh | A robot with improved absolute positioning accuracy for ct guided stereotactic brain surgery[END_REF]. The target was identified using computed tomography (CT) imaging.

Afterwards, many applications in medical robotics were developed, for example the RO-BODOC system (Think Surgical, Inc.) [START_REF] Taylor | An image-directed robotic system for precise orthopaedic surgery[END_REF] in 1992, which assisted a surgeon in a total hip arthroplasty procedure. This system was the first on its kind used on humans. It employs the position of three pins implanted in the hip by the surgeon, locating them in CT images for planning and performing most of the tasks involved in the orthopedic procedure. Since image-guidance offers extensive opportunities for medical robotics, more medical imaging modalities began to be employed in this field.

Currently, medical imaging has become essential in most medical procedures, providing examiners with the ability to see through the body without having to incise it (non-invasive). Along CT imaging, magnetic resonance imaging (MRI) and ultrasound imaging (US) have been the most frequently used modalities in medical robotics. MRI is a technology that uses strong magnetic fields to reconstruct images of the organs in the body. The acquisition of one image takes time and requires big sized equipment, making this modality incompatible for real-time visual control of medical robots. However, robots have been designed to be used in MRI scanner in order to assist in different tasks as for example the manipulation of a biopsy needle [START_REF] Su | Real-time mri-guided needle placement robot with integrated fiber optic force sensing[END_REF]. On the other hand, ultrasound (US) is a modality based on the propagation of the sound inside the body to generate images of the organs. The equipment used to perform ultrasound imaging is small compared to the other modalities. In addition, US has real-time acquisition capability, which INTRODUCTION makes it the most used technology for the image-guidance of medical robots. One of the first medical robots lead by ultrasound imaging was introduced in 1999 [START_REF] Salcudean | A robot system for medical ultrasound[END_REF]. This robotic system was designed to assist in the examination of the carotid artery by automatically moving an ultrasound probe such that the artery section was always visible in the ultrasound image. This system was the first one to use vision-based robot control (visual servoing) with ultrasound images. Since then, the application of ultrasound with visual servoing has gained great interest [START_REF] Azizian | Visual servoing in medical robotics: a survey. part II: tomographic imaging modalities -techniques and applications[END_REF].

New imaging modalities, such as elastography, have emerged of the analysis of the information obtained from the classical imaging modalities (ultrasound, MRI, CT, etc.).

Elastography introduces the promising concept of obtaining quantitative values of the stiffness of a tissue. The measurement of the stiffness can provide the examiner with more tools for diagnosis of diseases avoiding invasive approaches. This concept has been explored during the past three decades in medicine for the diagnosis of breast tumors [START_REF] Goddi | Breast elastography: A literature review[END_REF], liver fibrosis at different stages [START_REF] Barr | Elastography assessment of liver fibrosis: Society of radiologists in ultrasound consensus conference statement[END_REF] and prostate cancer [START_REF] Correas | Ultrasound elastography of the prostate: State of the art[END_REF]. However, the elastography process is currently performed manually, requiring high experience and training of the examiner. To overcome this issue, robotic systems can be used to assist in this medical procedure.

Motivations

Manual palpation is a medical procedure that has been used in diagnosis for centuries, in which the stiffness of the tissue of a patient is felt with the examiner's hands. It allows to recognize changes on the stiffness of the tissue, indicating a possible disease.

This practice is non invasive, simple in concept and needs no equipment. However, it requires great expertise and has significant constraints: it provides only qualitative information, it can be affected by the surrounding tissue and it is limited to the tissues within the reach of the examiner's hands. Palpation is used in the diagnosis of a wide range of diseases. Some illnesses of the liver, such as cirrhosis and hepatitis can be diagnosed by detecting fibrosis. Certain types of cancer, such as breast, prostate and thyroid cancer can be first identified by palpation. The early detection of cancer is fundamental in increasing the probabilities of successful treatment. However, the procedures most commonly performed for cancer detection, such as biopsies and mastographies, present significant drawbacks regarding the comfort of the patient, which can provoke for such patient to delay necessary examinations. Therefore, the development of tools for diagnosing which are mostly painless and non-invasive is of the utmost importance.

Ultrasound elastography is a tactile imaging method that can measure the stiffness INTRODUCTION of tissue using ultrasound waves. Consequently, it can overcome the limitations of manual palpation, providing important and precise information. Elastography is a process commonly performed by hand, which requires a repetitive motion applied to the tissue.

The generation of a stiffness map of a tissue can variate if the pressure applied on the tissue is not regular. As a consequence, the reproducibility of the results in elastography can be affected if several examinations are performed on the same tissue, specially if it is done by different examiners. This shortcoming suggests the need for an innovative method capable of providing more reliable information. On that account, a robotic system is capable of performing a repetitive task with the same pressure, which makes it a great option in the assistance for elastography process. Moreover, if the elastography is well performed, the output elasticity map of the tissue can be reproduced. Finally, the implementation of a robotic system capable of measuring the stiffness of a tissue in a patient can far extend the capabilities of diagnosis by manual palpation, creating a new tool to aid an examiner in improving the treatment of patients. The main goal of the thesis is to provide a general robotic control framework to assist an examiner in the elastography process. A general overview of the robotic framework we propose is illustrated by the block diagram presented in Figure 2 and described as follows:

Goal of the thesis

1. Development of a robotic system to perform the repetitive palpation motion on tissues, which is always needed in classic ultrasound elastography process. This goal requires the design of a robotic controller that applies motion with a conventional ultrasound transducer.

INTRODUCTION

2. Estimation of quantitative elastic information of the tissue in real time by the development of a fast elastography process.

3. Use of the elastic information of the tissue to perform robotic tasks, which can assist the examiner during the elastography process.

4. Teleoperation of the ultrasound probe to remotely explore the tissue.

Contributions

This dissertation presents several contributions to medical robotic assistance, which are listed as follows:

• A new methodology to use the strain map of the tissue as input of a visual servoing control scheme. This contribution includes the real-time estimation of the strain map, and the extraction of visual features required for the image-based control.

This methodology is developed for the cases of 2D and 3D ultrasound information.

• A complete system performing automatic palpation and providing quantitative ultrasound elastography. The proposed control process is composed of three hierarchical robotic tasks collaborating with each other. These tasks are proposed with the goal of always obtaining the strain map visibility of a region of interest. According to their priority, the three tasks are palpation motion by force control, automatic centering of a stiff tissue target by visual servoing and the orientation of the probe for tissue exploration. This contribution was published in two articles in the proceedings of the two major international conferences in robotics, IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) [START_REF] Patlan-Rosales | Automatic palpation for quantitative ultrasound elastography by visual servoing and force control[END_REF] and IEEE International conference on Robotics and Automation (ICRA) [START_REF] Patlan-Rosales | A robotic control framework for 3-D quantitative ultrasound elastography[END_REF].

• A new method to estimate the tissue strain map. This approach is based on the estimation of the motion displacements of a region of interest inside an ultrasound image. The displacements are computed through the geometric parameters involved in a deformable registration system.

• A new approach to estimate the strain map in a region of interest under motion perturbation using a 2D ultrasound probe. The motion compensation is based on a dense visual servoing approach that actuates the 2D ultrasound probe such that the perturbation motion is canceled. As a result, a robust estimation of the tissue elasticity can be performed on moving tissues. This contribution was published in an article in the proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS) [START_REF] Patlan-Rosales | Strain estimation of moving tissue based on automatic motion compensation by ultrasound visual servoing[END_REF].

INTRODUCTION

• A haptic system using the tissue elastogram that allows to the user to "feel" the elasticity of the tissue while changing the position of the region of interest in the ultrasound image.

• The combination of teleoperation of the ultrasound probe with the haptic system described in the previous contribution.

Structure of the thesis

The manuscript of this thesis is organized as follows:

Chapter 1 introduces the basic concepts and state of the art of ultrasound elastography. The chapter is divided in three main sections. First section starts with a brief description of the physics of the ultrasound beam forming, and it ends with the geometry of the ultrasound image. Second section presents the elementary concepts of elastography and a state of the art focusing particularly on ultrasound elastography. Last section recalls the visual servoing principle, which is used along the next chapters.

Chapter 2 details the robotic palpation system we propose to perform quantitative ultrasound elastography. This chapter introduces the experimental setup considered in this thesis, detailing all the equipments and the workflow used in the implementation of the robotic system. Afterwards, three main sections describe the three proposed robotic tasks in the design of the robotic palpation system for ultrasound elastography. The first task corresponds to an oscillatory force control of the ultrasound probe required for the elastography process. The second task involves the extraction and use of geometric parameters of a stiff tissue that is then automatically centered in the field of view of an ultrasound probe so the stiff tissue is always visible. The third task automatically adjusts the orientation of the ultrasound probe to a desired angle introduced by the user.

As these three tasks are coupled, a hierarchical approach is presented to combine them into the proposed palpation system. Experimental results obtained by using either a 2D or a 3D ultrasound probe are presented at the end of the chapter.

Chapter 3 presents a robust motion compensation process used to estimate the strain map of a moving tissue. This functionality complements and improves the system presented in chapter 2. This chapter is divided in five main sections. First section describes some works related to the visual tracking of deformable tissues. Second section presents an overview of the visual tracking models tested with real ultrasound images.

Then, we detail our proposed dense visual tracking system. This section also explains INTRODUCTION the use of the tracking model to estimate the strain map of the moving tissue. Afterwards, section four presents the control system used to compensate the motion using a 2D ultrasound probe actuated by a robot. In this section, the motion compensation is based on a dense visual servoing process. The last section presents and discusses the experimental results obtained from experiments performed with a phantom simulating moving tissues.

Chapter 4 proposes a haptic system based on the strain map of the tissue to provide the user with the feeling of the tissue elasticity. This chapter is divided in four sections. The first section describes the basic concepts of haptic feedback leading to the type of haptic feedback used in this chapter. The second section presents the model of the transformation from strain map to haptic force feedback. This section includes the description of a teleoperation system, which helps the user in the exploration of the tissues with an ultrasound probe. The fusion of the teleoperation with the oscillatory force control is presented in the same section. Section four shows the experimental setup and results of the haptic system. This section also describes the implementation details of the haptic system followed by the conclusion of the chapter.

Chapter 5 delivers the general conclusions of this thesis and proposes several shortterm and long-term perspectives of this work.

CHAPTER 1 ULTRASOUND ELASTOGRAPHY AND BASIC PRINCIPLES OF VISUAL SERVOING

The quest of building efficient systems to assist in medical procedures has been increasing in the last decades. Medical imaging technologies as ultrasound imaging, magnetic resonance imaging (MRI), X-ray radiography and tomography have been used by physicians to facilitate the diagnosis of illnesses and their treatment, often through complex medical procedures. However, the images produced by every one of these technologies are vastly different in terms of content, appearance and resolution and their analysis needs training and experience. When searching for malign soft tissue in a subject, medical image analysis can be exploited for visual localization, however, the elasticity properties of the studied tissue can provide more precise and valuable information. Elastography has been developed with the aim of finding the elastic parameters of a tissue to help in the detection of tumors and other malign bodies by their stiffness. This technology has been implemented using ultrasound imaging and MRI. Since the ultrasound imaging modality offers the capabilities of real-time and portability, ultrasound elastography is a promising technique to be used in building a robotic-assisted system as the one presented in this dissertation.

This chapter presents the basic principles of ultrasound imaging in section 1.1, from the beam forming to the classic b-mode image. Afterwards, a state-of-the-art related to elastography is introduced in section 1.2 that mainly focuses on techniques applied for ultrasound elastography. In section 1.3, several robotic systems involving elastography are presented. This chapter also recalls the definition of visual servoing in section 1.4

which is needed for a better understanding of the systems proposed in the next chapters.

1.1. ULTRASOUND IMAGING

Ultrasound imaging

In medical applications, the use of ultrasound imaging has become one of the most widespread along the medical community, due to its significant advantages over other technologies, such as: no radiation, low-cost, portability, to mention a few. The basic principle of US imaging consists in sending several pulses of ultrasound into the body and waiting for echoes to return. The echoes are then processed to produce an image of internal structures of the tissue. The ultrasound pulses are mechanical waves created by a vibrating object and propagated by a medium. The energy of these waves traveling through the body is then attenuated, scattered and reflected, producing

echoes. An ultrasound wave is represented as a repetitive pattern of high and low amplitudes. The distance between two peaks is known as wavelength (λ), which represents the property of the wave repetition and it is defined as

λ = c f (1.1)
where c and f are the speed and frequency of the sound. Commonly, the speed of the sound in soft tissue is assumed to be constant at 1540 m s -1 . The frequency of sound f , as mentioned before, is in the range between 2 and 40 MHz.

Generally, ultrasound waves are generated by a piezoelectric transducer, which is driven with electrical pulses sent from an ultrasound machine. The piezoelectric transducer also has the function of receiving the echoes after they have traveled through the tissue and transforming them into electrical pulses for the ultrasound machine.

The propagation speed of ultrasound in homogeneous tissue depends on two particular properties of the considered tissue, the bulk modulus B and the density ρ. It is described by the following definition: The propagation speed then depends on the material of the medium where it is propagated. For example, the speed of the sound that propagates in soft tissue is 1540 m s -1 .

c = B ρ . ( 1
When traveling in bone, the speed of the sound is around 4000 m s -1 .

Tissue is composed of different materials that influence the propagation of the ultrasound waves. This propagation is described by the Snell's law. This law relates the ultrasound wave directions of incidence, reflection and transmission at the interface between two different materials. Figure 1.1 shows two different media with sound propagation speeds c 1 and c 2 respectively. Snell's law states that if the sound wave is propagated with an incidence angle θ i = 0 with respect to the interface, then the wave is reflected in the same medium and transmitted through the second medium with the angles θ r and In practice, reflections of the sound wave do not only occur at tissue boundaries.

θ t , such that sin θ i c 1 = sin θ r c 1 = sin θ t c 2 (1.3)
Tissues are inhomogeneous thus producing local deviation of density contributing to the reflected wave. This is denominated as scatter reflection, which is represented as a collection of point scatterers retransmitting the incident sound wave in all directions. the wavelength of the sound wave induced. Otherwise, the sound wave induced is reflected If the scatterer is smaller than the wavelength of the sound wave induced, then the interference is minimum. Otherwise, the interference is important, causing occlusions.

The ultrasound wave propagation is also affected by the traveled distance. Its energy is reduced due to the scatterers and the absorption of the media. This energy loss is characterized by the attenuation coefficient α (expressed in decibels), which is dependent on the frequency of the sound wave. For soft tissue, the attenuation is usually between 0.3 and 0.6 dB/cm/MHz.

Ultrasound beam formation

The generation and detection of ultrasound waves is performed through a piezoelectric crystal that vibrates when an electric field is applied or generates an electric signal when it is subject to a mechanical vibration. The crystal is embedded in a transducer, usually an array of 128 elements, used as a transmitter and as a detector at the same time. Commonly, every piezoelectric crystal is driven with a sinusoidal electric signal that results in the emission of an ultrasound wave with a given frequency. The energy of the ultrasound wave reflected by the tissues that comes back to the crystal is then modified due to the factors previously discussed about sound propagation.

A transducer toggles from the state of transmitter to receiver every time a pulse is emitted. The ultrasound wave echo received by the transducer is captured as function of time by the ultrasound station. This detected signal is often called radio frequency (RF) signal, because its frequency range corresponds to the one of the radio waves in the electromagnetic spectrum (see Figure 1.3). An RF signal is described by the following

y(t ) = A(t )cos(ω r t + φ(t )) (1.4)
where ω r is the carrier frequency, φ is the phase and A is the amplitude of the RF signal (also known as the envelope of the RF signal). The description of the RF signal is expressed in the analytic form, however in practice it is recorded using a sampling frequency f s . The frequency f s follows the Nyquist criterion,

f s > 2 f max , (1.5)
where f max is the maximum frequency of the RF signal y(t ). For typical commercial ultrasound transducers, the sampling frequency f s is set between 20 -40 MHz. In our work, we use the beamformed RF data from an ultrasonics machine. The amplitude of the RF signal is commonly represented by a 16-bits integer. The high frequency information is removed by envelope detection (see Figure 1.4). This 1.1. ULTRASOUND IMAGING process is computed as follows,

Y env (t ) = y(t ) 2 + H (y)(t ), (1.6)
where H (y)(t ) is the Hilbert transform of the RF signal y(t ) defined as,

H (y)(t ) = - 1 π lim ς→0 ∞ ς y(t + τ) -y(t -τ) τ d τ. (1.7)
The envelope detection is also known as demodulation of the RF signal. In practice, the envelope of the RF signal is obtained in the Fourier domain as

F(H(y))(ω) = σ H (ω) • F(y)(ω), (1.8)
where F denotes the Fourier transform and σ H is a sign function:

σ H (ω) =      i , for ω < 0 0, for ω = 0 -i , for ω > 0 (1.9)
The envelope detection of the RF signal is used in the ultrasound image reconstruction.

Ultrasound image reconstruction

The envelope of the RF signal provides a filtered signal after removing the carrier frequency of the original signal. Every amplitude along the envelope is represented as a gray value of 8-bits. Then, for a transducer scanning different lines (number of scan lines SL), a brightness mode (b-mode) image is obtained as shown in Figure 1.5a.

In Figure 1.5a we can only observe specular reflections in the brightest areas. However, reflections due to the scatterers are barely visible due to the large difference of amplitude between the reflections. This issue can be easily addressed by performing a logarithmic compression of every scan line obtained from the transducer as

Y l og (t ) = A log(Y env (t )) + β (1.10)
where A is the amplification parameter and β is the linear gain parameter. This process enhances the image contrast as shown in the example presented in Figure 1.5b. Ultrasound image formation also depends on the geometry of the ultrasound transducer, more commonly known as ultrasound probe. There exist different ultrasound probe shapes, however two of the most common shapes are linear and convex probes.

The linear probes are mostly employed for vascular imaging and the convex probes are used in abdominal imaging. In the following, we explain the relation between the geometry of a 2D ultrasound probe and the image formation.

There are several parameters to consider in the reconstruction of an image from the geometry of the probe. For example, in case of a linear probe the transducer elements are co-linear to the extreme of the probe's surface and the scan lines are parallel to each other (see Figure 1.6). In the case of a convex probe, the transducer elements are positioned along the arc of the surface of the probe and the direction of the scan lines is therefore normal to this curved surface (Figure 1.7).

The geometric parameters of a linear probe are shown in Figure 1.6. Every point inside the RF array is represented by (i , j ), where i and j are the indexes of the scan line and the sample, respectively. Using those parameters, the metric coordinates (x, y) of any point located in the RF array at coordinates (i , j ) can be computed by

x = α L • (i -i 0 ) (1.11) y = α A • j (1.12)
where α A and α L are the axial and lateral resolutions, respectively. α L is the distance between two consecutive elements of the transducer, usually a measure given by the manufacturer. α A is the distance value between two adjacent samples in one scan line, and it is defined as

α A = c f s (1.13)
In the case of a convex probe, the geometry of the probe is given in Figure 1.

y = r cos(θ), (1.15)
where r is the distance of the point to the origin of F p and θ is the angle with respect to y-axis.

The RF scan lines obtained from the convex probe are stored in a rectangular array also called RF array. The origin of this array is at the top-left corner as shown in Figure 1.7 (bottom-right). The distance r is computed using the radius of the ultrasound probe r p and the distance between the origin and the j sample as,

r = r p + α A j (1.16)
where α A is the distance between two adjacent samples for every scan line i . This value is obtained as in Equation (1.13). To obtain the angle θ for the i -th scan line, first we need to know the angular distance between two consecutive scan lines α θ . This angular distance can be computed as

α θ = α L r p (1.17)
where α L is the separation between two contiguous elements of the transducer. Therefore, the angular field of view of the convex probe is

Θ = SLα θ (1.18)
where SL is the number of scan lines recorded (typical values are 128 and 192). The limits of the angular field of view correspond to θ min = -Θ 2 and θ max = Θ 2 . The angle θ for the i -th scan line is then defined as

θ = θ min + i α θ . (1.19)
With these previous relations (equations 1.14 to 1.19) we can therefore associate the (i , j ) indexes of each sample stored in the RF array to its respective metric coordinates (x, y). If the envelope detection is applied to the scan lines, then the image reconstructed is known as pre-scan image with memory indexes coordinates (i , j ) (see Figure 1.8a).

However, in a pre-scan image, the geometric structures of the tissue scanned by a convex probe are distorted since the geometry of the probe is not taken into account. The rectification of the pre-scan image is called post-scan image as shown in Figure 1.8b, which takes into account the geometry of the FOV of the ultrasound probe.

ELASTOGRAPHY: STATE-OF-THE-ART

In case of a convex probe, every pixel coordinates (u, v) in the post-scan image are then transformed to the metric Cartesian coordinates (x, y) by the relations:

u = x -(r p + SN α A ) sin(θ min ) s (1.20) v = y -r p cos(θ min ) s (1.21)
where SN is the number of samples in a scan line and s is the scaling factor of the image that represents the size of a pixel in meters. After every pixel is mapped in the post-scan image, an image interpolation is then applied to fill the missing pixel intensities between the captured scan lines.

Elastography: state-of-the-art

Many diseases cause changes in the mechanical properties of tissues. Current imaging devices such as computed tomography (CT), ultrasound (US) and magnetic resonance imaging (MRI) are not directly capable of measuring these mechanical properties. However, this information can be obtained with elastography imaging techniques that consist in applying internal or external compression on tissues and measuring the resulted strain distribution from the image. This strain distribution is related to the tissue elasticity and generates a strain image of the underlying tissues. This section provides a state-of-the-art of the elastography process, starting from the principle of elastography.

Afterwards, an overview of elastography applied on different medical imaging modalities is presented, as shown in Figure 1.9. Since this thesis concerns ultrasound elastography, the state-of-the-art will be mainly focused on techniques dedicated to the ultrasound modality that will be presented in section 1.2.5.

Figure 1.9: Organization of the state-of-the-art of elastography. The number at the bottom of every block is the subsection number where the topic is described.

Elastography principle

Mechanical properties of an organ or tissue provide essential information for diagnosis in medicine. For example, a tumor or diseased tissue can be detected by its stiffness, generally perceived by palpation, which is a physical examination used in medical diagnosis that is performed by applying pressure with a hand or fingers to the surface of the body. However, this method is limited by the accessibility of the examiner to the tissue of interest, and it only provides qualitative information that can be distorted by surrounding tissues. Elastography is one of the existing approaches able to overcome these issues. The principle of elastography is to apply external compression on tissues and to measure their resulted displacements using medical imaging in order to estimate a quantitative image of the strains (also known as elastogram or strain map). In practice, the elastogram is derived from the analysis of the pre-and post-compression states of the tissue. The solution proposed by Ophir et al. [START_REF] Ophir | Elastography: A quantitative method for imaging the elasticity of biological tissues[END_REF] uses sound waves generated by a piezoelectric transducer array (usually RF signals) and models every wave as a succession of springs (see Fig. 1.10a). If an axial force is applied over the succession of springs (see Fig. 1.10b), the length in each spring will change according to Hooke's law as follows:

F = i k i ∆l i , (1.22)
where F is the force applied to the succession of springs in Newtons (N). k i is the stiffness of every i -th spring in (N•m -1 ), and ∆l i is the deformation of each spring in meters (m).

The strain value (without units) for each i -th spring is defined as

ε i = ∆l i l i (1.23)
Fig. 1.10b shows the Hooke's law scheme, where l i is the initial length of the i -th spring, l ′ i is the length after a stress is applied and ∆l i = l ′ i -l i is the difference of lengths. The relation between the change of length ∆l i and the strain value ε i is illustrated in Figure 1.11.

The analogy of Hooke's law for the successive springs is then adapted to the echo signals and the change of length ∆l is the time-delay ∆ t between the pre-and post-compression In practice, the time-delay ∆ t is usually obtained from crosscorrelation analysis between pre-and post-compression segments of RF signal. However, the amplitudes of the RF signal change between the different states of compression of the tissue generating a false estimation of ∆ t when using cross-correlation analysis.

In Equation (1.23) we can note that an optimal estimation of ∆ t is essential to compute the best approximation of the strain value ε. Thus, other ways to estimate ∆ t are: phase zero estimation (PZE) [START_REF] Pesavento | A time-efficient and accurate strain estimation concept for ultrasonic elastography using iterative phase zero estimation[END_REF], axial velocity estimation (AVE) [START_REF] Loupas | An axial velocity estimator for ultrasound blood flow imaging, based on a full evaluation of the doppler equation by means of a twodimensional autocorrelation approach[END_REF] and optical flow (OF) [START_REF] Pan | A two-step optical flow method for strain estimation in elastography: Simulation and phantom study[END_REF]. PZE computes the displacement ∆ t based on the estimation of the zero-phase between the pre-and post compressed RF signals represented in the Fourier domain. On the other hand, AVE and OF compute the velocity between the two RF signals. AVE uses the Doppler effect to calculate the velocity while OF uses a first order Taylor approximation.

The tissue elasticity variates among all the internal materials inside the body. Young's modulus E describes the tendency of any material to deform when a stress is applied on it:

E = σ ε (1.24)
where σ is the stress applied to the material, measured in Pascals (Pa). Young's modulus is also expressed in Pa, since the strain ε has no units. The strain value can be obtained with the analogy of the Hooke's law for successive springs. However, the estimation of the Young's modulus requires knowledge of the value of the stress applied, which cannot be measured when the process of pre-and post-compression is performed manually. One alternative solution to obtain the Young's modulus is the use of shear waves, which are elastic waves generated by an external actuator (vibration), natural physiological stress (e.g. breathing) or acoustic radiation force impulses (ARFI) [START_REF] Sandrin | Shear modulus imaging with 2-D transient elastography[END_REF]. In this approach, local elasticity is estimated from the phase of the displacement, instead of the amplitude. The propagation of a shear wave is related to the shear elasticity modulus as,

ρ ∂ 2 - → u ∂t 2 = G∆ - → u (1.25)
where ρ is the density of the medium, -→ u is the displacement vector and G (in kPa) is the shear modulus of the medium [START_REF] Vappou | Magnetic resonance-and ultrasound imaging-based elasticity imaging methods: A review[END_REF]. Thus, the shear wave velocity is related to the shear modulus through:

c s = G ρ (1.26)
As shown in Equation 1.26, the estimation of the local shear wave velocity c s is required to compute the shear modulus G. The average of the shear wave velocity is calculated from the phase shift of displacement between two locations and its distance.

ELASTOGRAPHY: STATE-OF-THE-ART

The existing elastography approaches mainly concern modalities such as magnetic resonance elastography (MRE), computed tomography elastography (CTE), optical coherence tomographic elastography (OCTE) and ultrasound elastography (USE) which are described in the next sections.

Magnetic resonance elastography (MRE)

MRE is a non-invasive medical imaging technique that reconstructs the stiffness of the tissue by imaging the propagation of shear waves using MRI. MRE assesses quantitatively the mechanical properties of the tissue [START_REF] Mariappan | Magnetic resonance elastography: A review[END_REF][START_REF] Muthupillai | Magnetic resonance elastography by direct visualization of propagating acoustic strain waves[END_REF]. The technology is becoming available as an upgrade on conventional MRI scanners. MRE has proven to be beneficial as a clinical tool for the diagnosis of diseases such as hepatic fibrosis, which increases the stiffness of liver tissue [START_REF] Venkatesh | Magnetic resonance elastography of liver: technique, analysis, and clinical applications[END_REF][START_REF] Venkatesh | MR elastography of liver tumors: preliminary results[END_REF] (see Figure 1.12). The process to obtain a MRE image to 500 Hz are induced in the tissue using an external actuator. Then, the shear waves are imaged inside the body using a MRI process. Finally, the imaged shear waves are processed to generate quantitative images of the stiffness of the tissue. MRE does not have real-time capability since the MRI takes at least 1s to generate the images of shear VISUAL SERVOING waves. However, MRE offers quantitative stiffness of all the tissue scanned. One important application of MRE is the measure of the stiffness of in-vivo brain tissue as presented in [START_REF] Vappou | Assessment of in vivo and post-mortem mechanical behavior of brain tissue using magnetic resonance elastography[END_REF].

Computed tomography elastography (CTE)

CTE is a technique to obtain the stiffness of the tissue based on computed tomography (CT) images. There are two ways to implement CTE. The first approach consists in using two CT images corresponding to the pre-and post-compression states. Then, the displacement map between the two images is computed with an image registration technique, and the strain map is estimated from this displacement. A feasibility study of this approach was presented in [START_REF] Kim | X-ray elastography: A feasibility study[END_REF], where a non-rigid registration method was used to estimate the displacement maps.

The second way to implement CTE is by performing an evaluation of the stiffness of different tissues using a tactile sensor. Those tissues are imaged with CT and then segmented in the images (automatically or manually). Then, the intensities in the image and the stiffness measured with the tactile sensor are related with a curve fitting process.

Afterwards, the CTE image is reconstructed to provide quantitative stiffness of the tissue in the CT image [START_REF] Sasaki | CT elastography: A pilot study via a new endoscopic tactile sensor[END_REF].

The major drawback of CTE is the radiation exposure of the patient. However, CTE offers different information, which sometimes cannot be observed with other imaging modalities.

Optical coherence tomographic elastography (OCTE)

OCTE is a modality based on optical coherence tomography (OCT), which is a medical imaging technique that uses laser light and reaches micrometer resolution. OCT is based on low-coherence interferometry, typically employing near-infrared light. It uses long wavelength light which penetrates into the scattering medium. The information that OCT provides has high resolution. In 1998 Schmitt [START_REF] Schmitt | OCT elastography: imaging microscopic deformation andstrain of tissue[END_REF] used the OCT information to obtain strain maps of micro tissues (e.g., skin of the finger). More recently, in 2013, Sampson et al [START_REF] Sampson | Optical elastography probes mechanical properties of tissue at high resolution[END_REF] improved the differentiation of tissues pathologies, such as cancer or atherosclerosis, based on Schmitt's work.

1.2. ELASTOGRAPHY: STATE-OF-THE-ART

Ultrasound elastography (USE)

USE is a modality that has been in development for about 25 years. The principle of elastography presented in section 1.2.1 is based on ultrasound imaging. USE was developed to reconstruct an image representing the elasticity of the structures observed in the field of view of an ultrasound probe. This technique has been applied using different approaches to measure elastic parameters. These approaches can be divided in two categories according to the measurement of the displacements: quasi-static and dynamic. We show in Figure 1.13 four of the existing approaches: quasi-static, remote palpation, transient and supersonic. For all the mentioned approaches, an overview of related works is provided in the following sections.

R(n) = (s pr e ⋆ s post )(n) = M m=1 s * pr e (m)s post (m + n) (1.27)
where ⋆ is the cross-correlation operator and * denotes the conjugate of the function. n is the sample displacement (also known as lag) in a certain range (usually [-N , N ]). The output of the correlation R(n) is a signal of length (2N -1), where the amplitude at every CHAPTER 1. ULTRASOUND ELASTOGRAPHY AND BASIC PRINCIPLES OF VISUAL SERVOING n-th displacement corresponds to the similarity measured value. Then, the optimal displacement ∆ t between the two RF signal segments can be estimated by

∆ t = arg max n R(n) (1.28)
The optimal displacement ∆ t occurs when the two signals are the most similar. An illustrative example of the cross-correlation is shown in Figure 1.14. We can observe that RF signals are high frequency signals, which can cause the repetition of some segments. This repetition can produce a wrong estimation of the optimal displacement ∆ t when using N >> M . The quasi-static ultrasound elastography approach offers a solution to this issue. It consists in applying a small compression to the tissue that usually represents about 1% to 2% of the length of the tissue [START_REF] Ophir | Elastography: A quantitative method for imaging the elasticity of biological tissues[END_REF]. This quasi-static compression improves the estimation of the displacement by limiting the range of search in the cross-correlation analysis, while reducing the computational cost.

Once the displacement is estimated, the strain value of the tissue can be computed.

As explained in the principle of elastography described in section 1.2.1, the strain value can be computed through Equation (1.23). However, if the strain value is calculated for every segment of RF signal it can cause abrupt changes in the strain profile as illustrated in Figure 1.11. A better strain profile can be obtained by computing the strain value for all the samples in each of the RF signal segments. However, this process will increase the computational cost. An approach based on least-squares (LSQ) was proposed to avoid the computation of the displacement for every sample [START_REF] Kallel | A least-squares strain estimator for elastography[END_REF]. Let us consider the samples in a RF signal segment as a vector n ∈ [n i , n f ] with n i and n f as the initial and final 1.2. ELASTOGRAPHY: STATE-OF-THE-ART samples of the segment, respectively. Then, a vector of displacements ∆ t can be defined as

∆ t = an + b, (1.29)
where a and b are the slope and the offset values of the segment, respectively. The Equation (1.29) can also be rearranged as follows

       ∆ t (n i ) ∆ t (n i + 1) . . . ∆ t (n f )        =        n i 1 n i + 1 1 . . . . . . n f 1        a b , (1.30)
and can be written also in a short form as,

∆ t = A a b . (1.31)
We can observe from Equation (1.23) that the relation between the displacement ∆ t and the strain value ε can be expressed as,

∆ t = εn, (1.32)
leading to redefine a = ε and b = 0. Therefore, the Equation (1.31) can be changed to obtain ε from the displacement of the segment as

ε 0 = (A T A) -1 A T ∆ t (1.33)
Obtaining the strain values by LSQ as presented in Equation (1.33) still requires the computation of the displacement for all the samples in the segment n. However, Kallel and Ophir [START_REF] Kallel | A least-squares strain estimator for elastography[END_REF] found a reduction of the first row of A T A -1 A T as a vector g(k),

g(k) = ξ(k) 1 -k+1 VISUAL SERVOING
The methodology presented for the quasi-static ultrasound elastography is the base for the strain estimation of one RF line. Moreover, this approach can be easily adapted for 2D and 3D ultrasound information [START_REF] Rao | Correlation analysis of three-dimensional strain imaging using ultrasound two-dimensional array transducers[END_REF][START_REF] Shi | Two-dimensional multi-level strain estimation for discontinuous tissue[END_REF]. Since motion estimation is essential for optimal strain measurements, the improvement of this estimation in elastography has been studied in multiple works, like the one presented in [START_REF] Pan | A two-step optical flow method for strain estimation in elastography: Simulation and phantom study[END_REF], where a method based on two-step optical flow measure is proposed to estimate the tissue strain map. This technique has been tested in vivo in the breast of myocardial stiffness [START_REF] Konofagou | Myocardial elastography -a feasibility study in vivo[END_REF][START_REF] Varghese | Ultrasonic imaging of myocardial strain using cardiac elastography[END_REF].

Remote palpation

There are more ways to exert force on tissue than the mechanical compression presented in the quasi-static approach. Remote palpation introduces the alternative of using ultrasound energy to remotely excite the tissue [START_REF] Nightingale | On the feasibility of remote palpation using acoustic radiation force[END_REF]. In this method, acoustic radiation force impulse (ARFI) is used to locally displace the tissue in order to reveal its mechanical properties. ARFI is created by focusing a high-intensity ultrasound signal on a small region generating a local force in the direction of the propagation of the ultrasound wave.

This force is proportional to the power absorbed by the medium at the focal region:

F = 2υI c (1.35)
where υ is the absorption coefficient, I is the ultrasound intensity and c is the longitudinal wave speed in the medium. This force generates a displacement in the tissue (typically a few micrometers). This phenomenon is the result of the latency of soft tissue in responding to the excitation making the response out of phase [START_REF] Fahey | Acoustic radiation force impulse imaging of thermally-and chemically-induced lesions in soft tissues: preliminary ex vivo results[END_REF][START_REF] Nightingale | Acoustic radiation force impulse imaging: in vivo demonstration of clinical feasibility[END_REF]. The same ultrasound transducer used for imaging can generate ARFI excitation.

In the conventional method of ARFI strain imaging, the transducer excites multiple locations in the tissue at a constant depth [START_REF] Fahey | Acoustic radiation force impulse imaging of thermally-and chemically-induced lesions in soft tissues: preliminary ex vivo results[END_REF]. These locations are spread along the lateral direction. Regular ultrasound then captures the axial tissue motion before and after excitation for each lateral position to form the strain image. It is also possible to monitor temporal response of the tissue after excitation by firing a series of ultrasonic imaging signals. However, the sequence of firing limits the frame rate of ARFI imaging and the size of the region of interest. Another limitation is the imaging depth, which is limited in a focused region [START_REF] Zhai | Acoustic radiation force impulse imaging of human prostates: Initial in vivo demonstration[END_REF].

The first in-vivo results of ARFI imaging were presented in [START_REF] Nightingale | Acoustic radiation force impulse imaging: in vivo demonstration of clinical feasibility[END_REF] and demonstrated the feasibility of ARFI imaging for clinical applications. In [START_REF] Fahey | Acoustic radiation force impulse imaging of thermally-and chemically-induced lesions in soft tissues: preliminary ex vivo results[END_REF], thermally-and chemicallyinduced lesions were imaged in ex-vivo tissues. In addition, ARFI has been employed

for other applications such as assessment of breast lesions [START_REF] Sharma | Acoustic radiation force impulse imaging of in vivo breast masses[END_REF], detection of prostate cancer [START_REF] Zhai | Acoustic radiation force impulse imaging of human prostates: Initial in vivo demonstration[END_REF] and for delineation of Radiofrequency ablation (RFA) [START_REF] Fahey | Acoustic radiation force impulse imaging of myocardial radiofrequency ablation: initial in vivo results[END_REF].

1.2. ELASTOGRAPHY: STATE-OF-THE-ART

Transient elastography

Transient elastography or pulsed elastography is a technique where an ultrasound transducer probe is mounted on the axis of a vibrator (see Figure 1.15). Vibrations of mild amplitude and low frequency (arround 50 Hz) are transmitted from the vibrator to the tissues via the transducer. This induces an elastic shear wave that propagates through the tissues. In the meantime, pulse-echo ultrasound acquisitions allow the propagation of the shear wave to be followed and its velocity to be measured. The stiffness of the tissue is directly related to the velocity of the shear wave propagation: the stiffer the tissue, the faster the shear wave is propagated [START_REF] Sandrin | Transient elastography: a new noninvasive method for assessment of hepatic fibrosis[END_REF].

Supersonic ultrasound elastography

Supersonic ultrasound elastography or supersonic shear imaging (SSI) is the most recent technique for ultrasound elastography. This approach uses the same principle as remote palpation to generate mechanical vibrations using ARFI. However, supersonic ultrasound elastography generates the radiation shear waves, allowing to compute an 1.2. ELASTOGRAPHY: STATE-OF-THE-ART elasticity map of the tissue up to 5000 frames/s [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF].

The generation of the supersonic ultrasound elastography requires special equipment for fast generation of the ARFI. After the generation of the pulses with ARFI at supersonic speeds, the system switches to ultrasound imaging, requiring a fast acquisition and processing of the ultrasound information too. The processing of the ultrasound data is based on a cross-correlation analysis that estimates the axial displacement map required to compute the strain map of the tissue. SSI has been used for the assessment of liver fibrosis [START_REF] Chen | Supersonic shearwave elastography in the assessment of liver fibrosis for postoperative patients with biliary atresia[END_REF] with promising results as a noninvasive technique for evaluate liver fibrosis in children. Regarding to breast lesions, SSI has demonstrated the capacity to be considered as a tool for breast cancer diagnosis [START_REF] Athanasiou | Breast lesions: Quantitative elastography with supersonic shear imaging-preliminary results[END_REF].

Comparison of the ultrasound elastography approaches

The main approaches of ultrasound elastography can be compared as shown in Table 1.1. In this Table, the elasticity of the tissue can be estimated applying compression to the tissue in three ways: using quasi-static compression, with external vibrator or through ARFI. We also present the comparison of the advantages and limitations for every approach.

Approach

Stress type Stress source Advantages Limitations

Quasi The advantages and limitations of the ultrasound elastography approaches presented in Table 1.1 offer us a perspective to select the most convenient method for our robotic system. In order to avoid special equipment (besides the robot), we opt for the quasi-static approach, which is compatible with all ultrasound probes. The main limitation of this approach is the necessity of an operator holding the ultrasound probe to perform quasi-static compression. However, this limitation can be overcome by having a robot performing this task with precision.

CHAPTER 1. ULTRASOUND ELASTOGRAPHY AND BASIC PRINCIPLES OF VISUAL SERVOING

The following section presents a brief overview of medical robotic systems involving ultrasound elastography.

Robotic-assisted systems for elastography

Very few investigations have been undertaken regarding the use of ultrasound elastography in robot-assisted procedures. These works are related to the field of minimal invasive surgery (laparoscopy). For example, a snake-like robot was presented in [START_REF] Sen | Enabling technologies for natural orifice transluminal endoscopic surgery (N.O.T.E.S) using robotically guided elasticity imaging[END_REF], where a micro ultrasound probe attached at the distal part of the robot was used to find hard lesions by palpation motion. This system controls an 11 degrees of freedom (DOF) robot in order to perform three types of motion: coarse positioning, fine positioning and palpation motion. This system was tested on a prostate phantom with some stiff regions, and the performance of the palpation system is displayed in Figure 1.18. The da Vinci surgical robot (Intuitive Surgical Inc.) has been used to obtain elastic information of a tissue of interest by controlling the motion of a laparoscopic 2D ultrasound probe [START_REF] Billings | System for robot-assisted real-time laparoscopic ultrasound elastography[END_REF]. This robot-assisted system applies a palpation motion that is mixed with the teleoperated motion of a laparoscopic ultrasound probe. Therefore, it allows to obtain the elastogram of the tissue while the surgeon teleoperates the ultrasound probe with the da Vinci robot (see Figure 1.19). In a similar framework [START_REF] Schneider | Remote ultrasound palpation for robotic interventions using absolute elastography[END_REF], a mechanical vibrator placed on the skin of the patient was used to replace the palpation through the controlled motion of an ultrasound probe.

In ultrasound remote palpation, a robotic system was built to control the contact forces between an ultrasound probe and the tissue. Then, ARFI were applied to obtain a measure of the elasticity of the tissue [START_REF] Bell | Force-controlled ultrasound robot for consistent tissue pre-loading: Implications for acoustic radiation force elasticity imaging[END_REF]. The block diagram of this approach is shown in Figure 1.20, where a proportional-integral-derivative (PID) controller was designed to Force control scheme of the robotic system presented in [START_REF] Bell | Force-controlled ultrasound robot for consistent tissue pre-loading: Implications for acoustic radiation force elasticity imaging[END_REF] for applying the ARFI required for the tissue elasticity measurement.

The works mentioned in this Section introduced robotic systems to assist in the elastography process. However, none of these robotic systems used the information of the CHAPTER 1. ULTRASOUND ELASTOGRAPHY AND BASIC PRINCIPLES OF VISUAL SERVOING elastogram as part of the robotic controller, which is one of the main contributions of this thesis that will be presented in Chapter 2. The following Section briefly recalls the principle of visual servoing since it will be considered in this thesis work.

Visual servoing principle

A brief overview on visual servoing is presented here to explain the basic concepts of this control technique. The reader can find a more detailed explanation of visual servoing in the two parts tutorial presented by François Chaumette and Seth Hutchinson [START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF][START_REF] Chaumette | Visual servo control, part II: Advanced approaches[END_REF]. Two main configurations related to the visual sensor position are used in visual servoing: eye-in-hand and eye-to-hand configurations (see Figure 1.22). In the eye-in-hand configuration, the visual sensor is placed on the robot, so its motion is guided by the robot. Alternatively, in the eye-to-hand configuration, the visual sensor is in a fixed remote location, and observes the moving robot's end-effector interacting with the scene.

Interaction matrix

The configuration in any robotic mechanism depends on the position of its joints q(t ) ∈ R p at time t , with p as the number of joints (see Figure 1.23). The number of joints p is also known as the number of degrees-of-freedom (DOF) of the robot. Since we will consider a robotic arm holding an ultrasound probe, the following expressions will be recalled for the eye-in-hand visual servoing configuration. Therefore, the pose of the visual sensor r(q, t ) can be linked to the joints position using the forward kinematic model of the robot. In the design of the control scheme, the variation of the robot's end-effector ṙ(t ) is related to the time variation of the features ṡ(t ). This relation is defined where v is the velocity screw vector of the end-effector obtained from the time variation of r. ∂s ∂t is the variation of the features in the environment through time. This means that in a static environment ∂s ∂t = 0. The remaining term, ∂s ∂r is called the interaction matrix of k × 6, also defined as,

L s = ∂s ∂r .
(1.38)

The interaction matrix links the velocity screw vector of the visual sensor, v, to the variation of the visual features ṡ(t ) as,

ṡ = L s v. (1.39)
We consider here the case where the feature variation is only due to the robot displacement meaning that ∂s ∂t = 0.

Control law

The relation between time variation of the error ė and the sensor velocity v can be computed using equations (1.36) and (1.39) as follows,

ė = L s v. (1.40)
Since the goal of visual servoing is to minimize e, the variation of the error ė is usually set as an exponential decrease of the error ė = -λe,

where λ > 0 is the gain of the control law. This gain can be set as a constant value or as a variable (e.g, adaptive gain) dependent on the current error value [START_REF] Kermorgant | Dealing with constraints in sensor-based robot control[END_REF], such as,

λ( e ) = (λ 0 -λ ∞ ) e - λ ′ 0 λ 0 -λ∞ e + λ ∞ , (1.42)
where λ 0 = λ(0) and λ ∞ are the gains for the smallest and highest values of e respectively, and λ ′ 0 is the gain slope at e = 0. Now, we can compute the velocity control law to be applied to the robot by using equations (1.40) and (1.41) to obtain,

v = -λL + s e, (1.43)
1.5. CONCLUSION where L + s ∈ R 6×k is the Moore-Penrose pseudoinverse of L s defined by

L + s = L ⊤ s L s -1 L ⊤ s when L s is full rank. However, if L s is square (k = 6
) and the det(L s ) = 0, then it is possible to invert L s , giving the velocity control law as v = -λL -1 s e.

We should notice that the interaction matrix L s cannot be known perfectly on a real system and an approximated value Ls is usually considered. To ensure an asymptotic stability of the system using Ls , the condition Ls L s > 0 must be valid as demonstrated in [START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF]. Therefore, the control law becomes v = -λ L+ s e (1.44)

Conclusion

This chapter has introduced the principles of ultrasound imaging, ultrasound elastography and visual servoing. First, we explained the process to reconstruct an ultrasound image from the acquired information of the ultrasound propagation through the tissue.

This introduction allows the reader to be introduced to the basic concepts of medical ultrasound imaging which are important for the understanding of this thesis. Mainly, the focus of this chapter is elastography, presented in Section 1.2.1, where the most widely used approaches for this process were introduced. Magnetic resonance elastography (MRE), computed tomography elastography (CTE) and optical computed tomography elastography (OCTE) require expensive and large equipment to be used. On the other hand, ultrasound elastography (USE) requires small and less expensive equipment which is already present in most medical facilities. The different methods used in USE were also presented in Section 1.2.5, where the principle and state-of-the-art of every approach was provided. The comparison of these techniques was presented in Section 1.2.5.5, which led us to select the classic quasi-static approach for our robotic framework that will be detailed in the following chapters. The quasi-static approach was chosen due to its compatibility with most ultrasound systems, requiring no additional devices, nor special ultrasound transducers. All the concepts presented in this chapter are widely used in the next chapters. Chapter 2 will present a novel approach to build a robotic-assisted elastography system. This approach not only helps in the generation of the tissue elasticity map, but it also exploits the elastic information in a visual servoing task.

CHAPTER 2 AUTOMATIC PALPATION FOR ULTRASOUND ELASTOGRAPHY

In the previous chapter we presented the background of elastography and its first usage in robotic applications. Strain information has been used to localize malign tissue based on its elasticity, which is not possible with common b-mode ultrasound images. This chapter focuses on the development of a quasi-static elastography approach that ensures the compatibility with widespread conventional medical ultrasound systems. The major contribution presented in this chapter is the development of a robotic assistant palpation system that autonomously provides tissue elasticity information based on a quasi-static elastography estimation process. Our assistant robotic system is composed by a robotic arm that holds and controls an ultrasound probe. Its goal consists in continuously applying a palpation motion on the tissues and maintaining the visibility of a stiff tissue of interest in the ultrasound image. To achieve these assistance functionalities, we propose to perform three hierarchical robotic tasks that collaborate together. The first task consists in automatically applying a periodical compression motion to the tissues with the ultrasound probe in order to obtain the pre-and post-compression states of the tissues. The secondary task is based on a visual servoing control scheme that uses directly the strain information as visual feature to automatically maintain a selected tissue of interest in the field of view of the ultrasound probe. To the best of our knowledge, it is the first time the elastography strain information is used as input of a robot controller.

The third proposed task is the automatic orientation of the ultrasound probe, allowing the user to explore the surrounding area of the target tissue.

This chapter is structured as follows. In Section 2.1, the experimental setup used for 2.1. EXPERIMENTAL SETUP the experiments presented in this chapter and Chapter 3 is first detailed. Then, in Section 2.2, the process to obtain the elastogram by applying periodical soft tissue deformation by a force control is described. Section 2.3 develops the visual servoing approach we propose to automatically align the tissue of interest with the center of the FOV of the ultrasound probe. The automatic orientation of the ultrasound probe is presented in Section 2.4, and the fusion of the three control tasks is detailed in Section 2.5. Experimental results obtained with the robotic system are presented in Section 2.1 for the use of 2D and 3D ultrasound probes interacting with different kinds of phantoms.

Experimental setup

This section presents the setup and the equipment used for the experiments. First, we explain the workflow presented in Figure 2.1, and then we present all the components used to build the setup. In Figure 2.1, we illustrate the connections between the elements of our experimental setup. We have a robot equipped with a force/torque sensor and an ultrasound probe attached to the end-effector. Force sensor data is sent to a workstation, where all the algorithms are implemented. The robot is connected bidirectionally to the workstation to send its status and to receive motion velocity commands.

Ultrasound equipment

We used the ultrasound station SonixTOUCH (BK Ultrasound, MA) as shown in Figure 2.3a. This diagnostic ultrasound system is packed with an ultrasound research The probe used for the experiments is the 4DC7-3/40 (see Figure 2.3b) which is a convex 3D probe. This probe has a motor to orientate a 2D curvilinear transducer array.

The motor, with a radius of 2.72 cm, makes the probe suitable to acquire 2D and 3D ultrasonic data. The curvilinear array of the transducer has a frequency range between 7 and 3 MHz, a focal depth range between 5 and 24 cm, and an image field of view of 78°.

The motor sweeping maximum angle range is 75°.

Phantoms

Automatic palpation

We define as "autonomous palpation" the robotic task and image processing needed to automatically compute a strain map. It is performed in real-time by applying periodic compression motion to the tissue with an ultrasound probe attached to the end-effector of the 6-DOF Viper manipulator robot. The Figure 2.5 shows the location of the Cartesian frames we considered in this study, where the force sensor and the robot end-effector are positioned at frames F f and F e , respectively. We define two frames attached to the mechanical part holding the ultrasound probe that is plugged to the robot end-effector:

first, the frame at the gravity center of mass F g , second, the frame at the probe first contact point F cp . We define one frame F U S at the center of the image acquired by the ultrasound probe. All these frames are defined using the metric system (meter and radian) and we consider in the following of this manuscript that both intrinsic and extrinsic parameters of the ultrasound probe have been calibrated using a method like the ones presented in [START_REF] Lasso | Plus: Opensource toolkit for ultrasound-guided intervention systems[END_REF].

Force control

In order to obtain the pre-and post-compression states of the tissue, we propose to apply a varying force along the axial direction (y cp ) of the ultrasound probe at the contact frame F cp by using a force control scheme. We denote the velocity screw vector of the probe at the frame F cp , as v = [v x v y v z ω x ω y ω z ] ⊤ . The first three components of v correspond to the translational velocities, and the last three elements to the angular velocities.

The 6-axes force/torque sensor provides a force tensor measurement H f f in the force sensor frame F f . To measure the interaction force between the probe and the tissue expressed in the contact frame F cp , we have to consider the probe mass m p in order to compensate the gravity force tensor H g g = [0 0 9.81m p 0 0 0] ⊤ defined at F g . The force tensor applied to the tissue can then be expressed in the frame F cp as follows:

H cp cp = F cp f H f f -F f g H g g (2.1)
where F f g and F cp f are force twist transformation matrices from the gravity frame F g to the frame F f and from the frame F f to the frame F cp , respectively. The force twist transformation matrix is used to transform the force/torque vector expressed at a frame F b into a frame F a and it is defined by a 6×6 matrix: Since our goal is to control only the force component along the y-axis (axial direction) of the probe, we define the feature vector to be regulated as

F a b = R a b 0 3×3 t a b × R a b R a b (2.2)
s f = [0 1 0 0 0 0] H cp cp .
(2.3)

In order to apply a continuous compression motion, we propose to implement the following desired varying force that is based on a sinusoidal function:

F d (k) = ∆ F 2 sin (4k -T )π 2T + 1 + F 0 , (2.4)
where k is the discrete time and ∆ F is the amplitude of the sinusoidal function as shown in Figure 2.6. T is the period of the desired force signal expressed in sample time and F 0 is the initial desired force value. In order to apply this varying force along the y-axis, we define the desired force as

s * f = F d (k)
and the force error to minimize as e f = s f -s * f . An exponential decrease of e f is achieved by imposing the desired error variation of the error such as ė * f = -λ f e f with λ f being the force control gain. To generate a velocity control law that minimizes e f , we need to express the interaction matrix L f which relates the variation of the force feature to the probe velocity tensor such as ṡ f = L f v. In this work we consider an approximation of the interaction matrix L f = [0 K 0 0 0 0], where K is a coarse estimation of the contact stiffness between the probe and the tissue. The force control law is then obtained by applying the following velocity to the ultrasound probe:

v f = L + f ( ė * f + ṡ * f ), (2.5)
where the operator "+" represents the Moore-Penrose pseudo-inverse defined as

L + f = (L ⊤ f L f ) -1 L ⊤ f when L f is full rank. ṡ * f (k) = ds * f (k) dk
is the differential of the desired force variation. If we analyze Equation (2.5), the term ṡ * f (k) can be neglected in this work, since our goal is just to obtain a sinusoidal variation of the force. Simulated outputs of the force controller are presented in Figure 2.7 to show the temporal evolution of the resulting force with and without considering the term ṡ * f (k) in the control law. Since a perfect phase between the desired and measured force is not necessary for the palpation motion task, we simplified the control law as:

v f = L + f ė * f . (2.6)
The velocity v f is applied in the contact point frame F cp . However, in practice the velocity must be expressed at the robot's end-effector frame F e to be applied. Thus, the velocity at F e is formulated as follows, where V e cp is a 6×6 velocity twist transformation matrix defined in Equation (2.8).

v e = V e cp v f (2.7)
V e cp = R e cp t e cp × R e cp 0 3×3 R e cp
(2.8)

The velocity v e can finally be applied using the robot's velocity kinematics model as follows, q = J + e v e (2.9)

where q is the join velocity control vector and J e is the Jacobian of the robot estimated through its kinematics. The remaining of the manuscript assumes that Equations (2.7)-(2.9) are always applied after expressing any control law at the frame F cp .

Elastogram estimation

Force control gives us the mechanical compression required for the elastography. The process that estimates the elastogram is illustrated in The elastogram is generated using a method based on motion estimation and strain filtering. First, we detail the motion estimation process for the 2D and 3D cases, and then the filter used to estimate the strain map.

2D Motion estimation

The essential part of the elastography method is the displacement estimation of the RF signals from pre-to post-compression states. This can be computed using motion estimation. We propose to use a subpixel motion estimation approach [START_REF] Chan | Subpixel motion estimation without interpolation[END_REF] with the purpose of achieving real time elastography imaging capability. Motion estimation is divided in two steps: integer displacement estimation and sub-displacement estimation. Integer displacement estimation is obtained with the block matching algorithm (BMA). Figure 2.9 shows the parameters used for this approach. Let us define two arrays of RF signals, the RF ROI in pre-compression as f (i , j ) and post-compression as g (i , j ), where i is the scan-line index and j is the sample index of the RF scan line. BMA 2.2. AUTOMATIC PALPATION divides the RF frame in blocks of M × N size (B 1 and B 2 for f (i , j ) and g (i , j ) respectively). Then, the displacement for each block in the current frame f is estimated with respect to the next one, g . The search of the best match is performed over a region of size (2N -1)×(2M -1) (as the one shown in Figure 2.9 in orange). The search region size can be changed to optimize the computational cost.

The best match of the i -th block over the search region can be found using common similarity measures (MSE, SAD, ZNCC, etc.). In our case, the sum of absolute differences (SAD) was selected due to its low computational cost. Therefore, the minimization is computed as follow,

(u 0 , v 0 ) = arg min 1 M N M -1 m=0 N -1 n=0 |B 1 (m + u 0 , n + v 0 ) -B 2 (m, n)| (2.10)
where B 1 ∈ f and B 2 ∈ g are blocks (matrices) of M ×N size. The u 0 and v 0 are the integer displacements corresponding to lateral and axial displacements, respectively.

Unfortunately, BMA only estimates integer displacements, and we need to estimate the sub-displacements to obtain an accurate displacement map. There are approaches for sub-displacement estimation that are based on parabolic interpolation, cosine interpolation, optical flow (OF) and splines, to name a few. In our work, we use OF as base for sub-displacement estimation. The OF between the blocks B 1 and B 2 is estimated by solving the next linear system,

  m,n ∂B 1 ∂i 2 m,n ∂B 1 ∂i ∂B 1 ∂ j m,n ∂B 1 ∂i ∂B 1 ∂ j m,n ∂B 1 ∂ j 2   δ u δ v = m,n (B 2 -B 1 ) ∂B 1 ∂i m,n (B 2 -B 1 ) ∂B 1 ∂ j (2.11)
Once we have obtained the sub-displacements (δ u , δ v ), we can compute a more accurate displacements as, The overlapping can be tunned to achieve good results. For example, a real-time displacement map was obtained using a zero-overlapping in [START_REF] Zhou | A motion estimation refinement framework for real-time tissue axial strain estimation with freehand ultrasound[END_REF]. This produced good results when a 2% compression between the pre-and post compress states was applied. In our case, we consider 25% of overlapping to obtain a finer displacement map while allowing real-time processing capability. The displacement maps, U 0 (i , j ) and V 0 (i , j ) for lateral and axial displacements respectively, are the outputs of this process (see Figure 2.11). Since the maps are of different size than the RF ROI due to the overlapping, we apply a bilinear interpolation to obtain the real sized U (i , j ) and V (i , j ) of the considered ROI. After obtaining the displacements for every block, their axial and lateral components are stored in two arrays V 0 (i , j ) and U 0 (i , j ), respectively.

u = (i c -u 0 ) + δ u (2.12) v = ( j c -v 0 ) + δ v (2.

AUTOMATIC PALPATION

Strain filtering

Strain values depend directly on the motion estimation. The base of the strain tensors is the partial derivative of the directional displacement map with respect to the axis of the tensor component. In the specific case of elastography, the compression is performed in axial direction. Then, the axial strain values can be computed as,

ε j j = ∂v ∂ j . (2.14)
where ∂v ∂ j indicates the variation of the axial displacement v through j -direction.

Additionally, to obtain a better quality elastogram, we use the LSQ strain filter proposed by Kallel and Ophir [START_REF] Kallel | A least-squares strain estimator for elastography[END_REF]:

h(n) = ξ(n) 1 -n+1 2 1 2 . . . n 1 1 . . . 1 (2.15)
where

ξ(n) = 12 n(n 2 -1)
and n is the number of samples in the interval ∆ j as shown in Figure 2.12. As demonstrated in [START_REF] Kallel | A least-squares strain estimator for elastography[END_REF], the convolution of this filter h(n) with the axial component of the motion estimation V (i , j ) can generate a smooth strain map ε(i , j ).

The strain map ε(i , j ) provides the elastic information inside the ROI. However, the strain is a measure depending on the constant stress applied to the tissue. This stress is performed with the ultrasound probe by applying the force controller presented in Section 2.2.1.

Evaluation of the estimated elastogram using a ground truth from finite element model (FEM) simulation

The process to obtain the elastogram should be evaluated based on a ground truth model. Therefore, we propose to evaluate the strain map estimation based on finite element analysis (FEA). The physical phenomenon of strain is expressed using partial differential equations. Solving these equations for any shape using analytical formulation is challenging. However, FEA is a numerical methodology to approximate the solution to these partial differential equations. The principle of FEA is to divide a rigid body into finite elements using a mesh and compute a solution for every element. The accuracy of the solution improves with the increment of the number of elements, but the computational cost also increases. FEA is used in many engineering applications to design and test mechanical structures under several boundary conditions. There is a wide variety of software for FEA, from open-source to paid license. In our case, we use COMSOL Multiphysics 5.0 (COMSOL, Inc.) to obtain the ground truth of the strain map. To estimate the elastogram as in section 2.2.2, we need the pre-and postcompressed RF frames. We obtain those RF frames using an ultrasound simulator called

Field II [START_REF] Jensen | Computer phantoms for simulating ultrasound B-mode and CFM images[END_REF]. This software can simulate the data acquired by a virtual ultrasound system based on an accurate ultrasound wave propagation model. We first defined a virtual stiff tissue target by randomly positioning scatterers in a geometrical mesh that represents a virtual organ at the pre-compression state as illustrated in Figure 2.13. Then we simulated tissue deformation by applying displacements computed from the FEA to the Once the models are defined for the FEM, we need to apply the compression. In our work, we fixed the edge of the base and define a Neumann boundary condition, which specifies the normal derivative of the function on a surface. A 2% compression was applied and we obtained the results reported in Figure 2.16 with the image of axial displacements presented in the first row and the strain image in the second row.

The value of the displacements d ∈ R 2 obtained by FEM are then applied to the scatterers position. We define the scatterers position at pre-compression frame as p s , and the scatterers position at post-compression as p ′ s . We use each frame as input in FIELD II to obtain its RF frames. Figure 2.17 shows the output of the three models described previously, where the pre-and post-compression frames are displayed in the first and second row, respectively. The output is shown in b-mode.

The output of FIELD II gives the RF data, which is the input to the processes described in Section 2.2.2.1 and Section 2.2.2.2. These processes were implemented in Qt with C++, and the library used to compute the FFT (Fast Fourier Transform) in the

Extension for 3D elastogram estimation

We have defined how to compute the elastogram in 2D. Now, in this section, we use an ultrasound probe able to obtain volumetric data. For this case, the ultrasound probe motor is enabled to move, with a specific angle step, as opposed to the 2D case where the motor position is fixed. In this case, one volume of RF signals can be seen as a set of N f RF frames in 2D. These frames are acquired during one directional sweep of the motor as shown in Figure 2.19.

If we estimate the 3D elastogram with the same principle as the 2D case, then we need to wait for the acquisition of two complete RF volumes, one for each state: preand post-compression. However, the computational cost of the motion estimation in process. First, we store the RF volume for one state V r (pre-or post-compress). Afterwards, the 2D elastograms are computed online for every motor position in the current motor sweep, where we also grab the frames in the RF volume V c . At the end of the current motor sweep we obtain the complete 3D elastogram, defined as,

V s (i , j , k) = ε k (i , j ) (2.16)
where i and j are the indexes of the ultrasound scanline and the sample in the scanline, respectively. ε k is the strain map in 2-D estimated at the k-frame pair of V r and V c in the interval [k 0 , k n-1 ] with k 0 and k n-1 as the indexes of the initial and final frame pair in the volumes V r and V c , respectively.

Automatic centering of a stiff tissue in ultrasound image

In the last section, we presented how to estimate the elastogram by automatically moving an ultrasound probe with a robot. Now, in this section, we define a robotic task to automatically align the center of the probe with the stiffest tissue in a ROI. This process can assist the examiner by always maintaining the visibility of the target tissue during a medical procedure. We propose to use visual servo control for this robotic task by considering visual features extracted from the estimated strain image directly as inputs of the control scheme.

In order to automatically center a stiff object at the middle of the full image by visual servoing, we propose to isolate the biggest rigid region in the elastogram and use its barycenter coordinates as the visual features.

The method we propose to extract these features consists first in segmenting the biggest stiff region from the elastogram and then computing the coordinates of its centroid as described next.

Stiff tissue segmentation

First, we propose to generate an image I g (i , j ) by filtering the strain map with a Gaussian function as:

I g (i , j ) = e -ε(i , j ) 2 e ε 2 max
(2.17) 53

AUTOMATIC CENTERING OF A STIFF TISSUE IN ULTRASOUND IMAGE

where ε max = max(|ε(i , j )|) and I g (i , j ) ∈ [0, 1]. The aim of this filter is to enhance the intensity of the rigid objects and to decrease the intensity of the rest of the area.

The segmentation is obtained by applying a binarization on the image I g (i , j ) such as,

I u (i , j ) = 1 if I g (i , j ) ≥ Γ 0 otherwise (2.18)
where Γ is the threshold value. In practice, the value of Γ = 0.5 was used since it provides optimal results.

A similar approach can be used for the segmentation of a 3D strain volume. In this case, we generate first a filtered volume V g (i , j , k) ∈ [0, 1] as,

V g (i , j , k) = e -V s (i , j ,k) 2 e V 2 smax , (2.19)
where

V s max = max(|V s |).
The segmentation of the volume V g is then computed as,

V u (i , j , k) = 1 if V g (i , j , k) ≥ Γ 0 otherwise (2.20)

Centroid estimation

Once we have isolated all the stiff regions, we need to know which is the biggest one by labeling all the regions. For this process, we can use a connected components algorithm [START_REF] Samet | Efficient component labeling of images of arbitrary dimension represented by linear bintrees[END_REF], which is a graph theory-based approach to label different regions. In Figure 2.21, we show an example of the labeling of different regions. We represent a binary image with a tree (e.g., a quadtree for a 2D case) to search the adjacent nodes and classify them in one region. Every node is defined as a pixel in the 2D case, or as a voxel in the 3D case.

The biggest region is the region having the bigger area (or volume in the 3D case) of the labeled regions, and its centroid is computed from the image moments as proposed in [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]. The general definition of the image-based moments in 2D is given by where the integration is done using all the elements in the RF plane Ω. The area A and the centroid (i c , j c ) are computed as follows,

M mn = Ω I u (i , j )i m j n d i d j , (2.21)
A = M 00 (2.22) i c = M 10 M 00 (2.23) j c = M 01 M 00 (2.24)
In a similar way for the 3D case, the moments in 3D are defined by,

M mno = Ψ V u (i , j , k)i m j n k o d xd yd z, (2.25)
and the integration is done using all the elements in the RF volume Ψ. The volume ν and the centroid (i c , j c , k c ) are estimated as,

ν = M 000 (2.26) i c = M 100 M 000 (2.27) j c = M 010 M 000 (2.28) k c = M 001 M 000 (2.29)
The area (volume for 3D case) of every labeled region is used to find the biggest one.

Then, the centroid of the biggest region is obtained through the image moments as previously defined. Figure 2.22 shows the different steps of the centroid computation in a 2D case.

c x = -α L i r 0 + i c - w p 2 (2.30) c y = α A j r 0 + j c - h p 2 (2.31)
where w p in RF lines and h p in RF samples are the width and height of the RF frame.

The coordinate (i r 0 , j r 0 ) is the initial point (top-left corner) of the ROI with respect to the RF frame. The constants α L and α A are the scale values to express RF in metric units for RF lines and samples, respectively. We should mention that α L and α A were defined in Section 1.1.2. Due to the geometry of the convex probe, the centroid (c x , c y ) is distorted, and it requires scan conversion to the real metric centroid (x c , y c ). This scan conversion is estimated through Equations (1.14) and (1.15).

To center the ultrasound probe with the stiffest tissue centroid, we need to minimize

x c to zero. In a 2D case, to achieve this robotic task we are controlling in-plane motions. This means that we can control only 3-DOF of the ultrasound probe, v x , v y and ω z . Then, we define the measure of the lateral component of the centroid as our measures in the control loop, s t = x c , and our desired centroid position as s * t = 0 corresponding to the horizontal component of the center of the ultrasound image. The relation of the variation of s t with respect to the ultrasound probe velocity v is defined as,

ṡt = L t v (2.32)
where L t is the interaction matrix that relates the probe velocity and the measure variation. The value of L t can be found through the Varignon's formula using the relation between the centroid variation with respect to the probe's velocity as,

   ṡt ẏc 0    =    ẋc ẏc 0    = -    v x v y v z    -    ω x ω y ω z    ×    x c y c 0    . ( 2

.33)

As we mentioned before, we want to control only 3-DOF (v x , v y and ω z ). Then, the remaining velocities are zero, and we compute the interaction matrix

L t = [-1 0 0 0 0 y c ].
The goal of this control task is to minimize the error e t = s t -s * t with an exponential decrease of e t . Therefore, the desired variation of the error is defined as ė * t = -λ t e t with λ t being the centering control gain. Then, the control law for automatic horizontal centering is defined as, 57

AUTOMATIC CENTERING OF A STIFF TISSUE IN ULTRASOUND IMAGE

v t = v = L + t ė * t (2.34)

Automatic centering in 3D

Figure 2.24: Axes for the frame F p (frame of the ultrasound probe). Lateral and top views of the probe in contact with a virtual human torso are shown in left and right images, respectively. A target, in green, is placed inside of the torso to relate the two views.

Similar to automatic centering in 2D, the goal is to center the stiffest object on the plane X -Z of the full volume (see Figure 2.24). First, we need to convert the value of the centroid in RF units to the metric coordinates. To do this, as we use a convex ultrasound probe, we perform a scan conversion of each point inside the RF volume to the Cartesian coordinates (see Figure 2.25), s(i , j , k) → p(X , Y , Z ), in order to obtain the metric location with respect to a Cartesian frame. We define the scan conversion using the ultrasound probe parameters as described in [START_REF] Lee | Intensity-based visual servoing for non-rigid motion compensation of soft tissue structures due to physiological motion using 4D ultrasound[END_REF]. In our case, RF data is considered instead of prescan images. We recall the scan conversion formulation as, The quasi-spherical coordinates are computed in function of the RF coordinates as,

r = v s f s j + r p (2.38) φ = -0.5α l (N f -1) + α l i (2.39) θ = -0.5η(N f -1 -2k) (2.40)
where v s is the speed of the sound (1540 m/s), f s is the sampled frequency (40 MHz), α l is the angle between neighboring scanlines and η is the angle of the field of view (FOV) of the motor in the ultrasound probe for a motor angular step. All these parameters are given by the ultrasound system and probe specifications.

Using the scan conversion defined above, we can compute the metric value of any point in the RF volume with respect to the probe frame F p . We calculate the centroid inside of a VOI, which means that the centroid coordinate in the full RF volume is expressed as,

i cm = i r 0 + i c + w p 2 (2.41) j cm = j r 0 + j c + h p 2 (2.42) k cm = k r 0 + k c + d p 2 (2.43)
where (i r 0 , j r 0 , k r 0 ) is the initial point of the VOI (top-left-front corner) in RF units. w p , h p and d p are the dimensions of the VOI in RF units. Then, using equations (2.35) to (2.40) for (i cm , j cm , k cm ), we obtain (X c , Y c , Z c ) in metric units.

Keeping a target in the FOV of the volume of analysis is a task which requires the displacement of the ultrasound probe on the X -Z plane (see Figure 2.24). The centroid of the target is defined as (X c , Y c , Z c ), but in this case we only use the values of X c and Z c . This means that the centroid's coordinates with respect to F p are the same with respect to the frame F cp . Then, similar to the 2D case, we define a visual feature as s t =

PROBE ORIENTATION

[X c Z c ] ⊤ and the desired feature vector to reach the centering of the object of interest in the probe FOV is directly s * t = 0 2×1 .

The error is defined as e t = s t -s * t . Similar to the 2D case, an exponential decrease of the error can be obtained by defining the desired error variation as ė * t = -λ t e t with λ t being the target-probe centering control gain.

Using the Varignon's formula, we determine the relation between the probe velocity v and the variation of the retained features as,

Ẋc Żc = -1 0 0 0 -Z c Y c 0 0 -1 -Y c X c 0 v. (2.44)
The Equation (2.44) can be written as ṡt = L t v, where L t is the interaction matrix related to s t . Then, the control law for the target-probe centering can be expressed as

v t = v = L + t ė * t (2.45)

Probe orientation

Probe orientation is the third proposed task in our approach. This task offers to the user the capability to explore the surrounding area of the target tissue with an automatic orientation of the ultrasound probe. In the following, we detail this robotic task in both 2D and 3D cases.

2D Probe orientation

CHAPTER 2. AUTOMATIC PALPATION FOR ULTRASOUND ELASTOGRAPHY

The aim of this task is to automatically orient the probe to a desired angle s * θ (in the image plane) from the current angle of the probe s θ = θθ i ni t . θ i ni t is the angle of the initial probe orientation and θ is the angle measured during the probe orientation control (see Figure 2.26). Both angles are obtained by the odometry measures of the robot. The variation of the angle feature s θ due to the probe velocity is defined as:

ṡθ = L θ v (2.46)
where L θ = 0 0 0 0 0 -1 is the interaction matrix related to the variation of s θ .

The angle error is defined as e θ = s θ -s * θ , and the desired angle error variation as ė * θ = -λ θ e θ with λ θ being the probe orientation control gain. Therefore, the control law for the orientation of the probe is defined as, The error to be minimized is defined as e θ = s θ -s * θ , and the desired exponential error decrease can be achieved by the desired error variation expressed as ė * θ = -λ θ e θ where λ θ is the orientation control gain.

v θ = L + θ ė * θ (2.47)

3D Probe orientation

As in the previous 2D case, we determine the interaction matrix that relates the feature vector variation, ṡθ , with the probe's velocity v as,

L θ =    0 3×3 -1 0 0 0 -1 0 0 0 -1    (2.48)
The control law for the probe orientation is then provided by,

v θ = L + θ ė * θ (2.49)

Control fusion

The system that we propose requires the three control tasks presented in this chapter.

However, if we analyze the interaction matrices, we can observe a coupling between the automatic centering and the probe orientation controls. This means that these two tasks are disturbing each other. We can deal with this through the redundancy control framework [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], where a hierarchical method for the i -th control task ( ėi , L i ) is proposed as,

v 0 =0 v i =v i -1 + (L i P i -1 ) + ( ėi -L i v i -1 ) (2.

50)

where P i -1 is the projection operator onto the null-space of (L 1 , . . . , L i -1 ), and it is defined as, P 0 =I

P i =P i -1 -L + i L i (2.51)
This formulation allows us to establish the control tasks priorities, giving to the ith task a lower priority with respect to the previous i -1 task so it does not disturb it.

Considering this hierarchical approach, we assign to the palpation task by force control e f the highest priority, since it is needed for the elastogram estimation process. Then, the automatic horizontal centering e t and the automatic probe orientation e θ are set as the second and third priorities, respectively. We can express these tasks using the redundancy control framework. First we designated the task e f with the highest priority such that,

v 1 = v f = L + f ė f . (2.52)
Then, the projector of v 1 on the next task is defined as,

P 1 = I 6 -L + f L f =           
1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

0 0 0 0 0 1            . ( 2

.53)

We can observe in P 1 that the first task, force control, constrains the next tasks on the translation along the y-axis. Then, projecting the matrix L t onto the null space of L f (ker(L f )) we obtain,

L t P 1 = L t . (2.54)
This means that the first and second tasks are decoupled, and the product L t v f = 0.

Then, the secondary task can be defined as,

v 2 = L + t ėt , (2.55)
where L + t is the Moore-Penrose pseudoinverse of L t that we analytically calculated for the 2D and 3D cases as,

2D case: L

+ t = C 1 -1 0 0 0 0 Y c ⊤ , (2.56)

3D case: L

+ t = C 1 C 2 -(X 2 c +Y 2 c +1) 0 -(X c Z c ) -(X c Y c Z c ) -Z c (Y 2 c +1) Y 2 c (X 2 c +Y 2 c +1) -(X c Z c ) 0 -(Y 2 c +Z 2 c +1) -Y c (Y 2 c +Z 2 c +1) X c (Y 2 c +1) X c Y c Z c ⊤ , (2.57)
where

C 1 = 1 1+Y 2 c and C 2 = 1 X 2 c +Y 2 c +Z 2 c +1 .
The projector of this task on ker(L f ) and ker(L t ) can be defined for both cases as,

P 2 = P 1 -L + t L t , (2.58)
producing the projectors that constrain the motion for the third task on the x-and ytranslations for the 2D case, and for the 3D case x-, y-and z-translations.

The third task, probe orientation, can be expressed as, 2.5. CONTROL FUSION

v 3 = (L θ P 2 ) + ( ėθ -L θ (v f + v 2 )). (2.59)
Finally, we obtain the control law that fuses the three hierarchical tasks with the following expression:

v = v 1 + v 2 + v 3 .
(2.60)

The behavior of this control law can be compared with another approach, where there are not task priorities, which can be defined by a simple interaction matrices stacking as:

2D case: v =    -1 0 0 0 0 Y c 0 K 0 0 0 0 0 0 0 0 0 -1    +    ėt ė f ėθ    (2.61) 3D case: v =            -1 0 0 0 -Z c Y c 0 K 0 0 0 0 0 0 -1 -Y c X c 0 0 0 0 -1 0 0 0 0 0 0 -1 0 0 0 0 0 0 -1            +       ėt x ė f ėt z ėθ       (2.62)
In the 3D case, ėt x and ėt z are the first and second row of ėt , respectively. We can observe Equation (2.61) (2D case) that the first and third rows of the interaction matrices stacking, representing the second and third tasks, respectively, are disturbing each other.

Similarly for the 3D case in Equation (2.62), we can notice this disturbance between the same tasks, the second task and third tasks represented in the first row and third to sixth rows, respectively. This issue represents a coupling between the second and third task which can cause instability in the control law when using equations (2.61) and (2.62).

However, we can see that the approach presented in Equation (2.60) can deal with the coupling of the second and third tasks offering us stability in the control system.

In order to demonstrate the difference between the control laws presented with and without the hierarchical approach, we performed a simulation with the initial and desired features shown in Table 2.3. This table presents the parameters of the controller for the 2D and 3D cases. The initial values are set for all the parameters. However, we changed at t =2s the desired probe orientation in order to observe the behavior of the control laws due to the modification of the third task. The evolution of the probe velocities is shown in Figure 2.28.

Modality

v x v y v z ω x ω y ω z Figure 2
.28: Plots of the probe control velocities. First and second rows are for the 2D and 3D cases, respectively. The performance of the three tasks are displayed using a global interaction matrix at the left, and using the redundancy control framework at the right. The gray strip highlights the time segment where we can observe the difference between the two approaches.

EXPERIMENTAL RESULTS

We can observe in Figure 2.28 the differences between the two approaches that occur during the time periods that are highlighted with the gray strips. As expected, the hierarchical control framework allows to obtain a perfect decoupling of the three tasks.

On the other hand, the use of the interaction matrices stacking, presented in equations (2.61) and (2.62), shows the coupling between the second and third tasks as we expected.

For this reason, we will retain the hierarchical approach in the next of our work.

Experimental results

We present here the experimental results obtained by applying the control framework that combines the three robotic tasks. We consider both the 2D and 3D cases that correspond respectively to the use of a 2D or 3D ultrasound probe. In the first part of this section, we present the results and performance obtained for the 2D case by considering in the setup the abdominal phantom that was presented in Figure 2.4a. Then, in the second part of this section, we present the results obtained with the implementation of our control approach to the 3D case by the use of a 3D US probe interacting with a homemade phantom.

For all experiments, the acquisition of RF data was implemented using a server-client TCP/IP communication in a local network. We used as server the SonixTouch ultrasound scanner, and as client a Linux workstation (Intel Xeon CPU @2.1 GHz) that performs all the imaging process, control law computation and communication with the robot. The RF data from the server is sent to the client at the rate of 24 FPS (frames per second).

In the 3D case, each 3D image (volume) is composed of 31 RF 2D frames which result in a volume acquisition rate of 0.77 VPS (volumes per second). On the client side, we developed a multi-thread software application in C++ based on ViSP [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF], VTK [START_REF] Schroeder | Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition[END_REF] and Qt [START_REF] Eng | Qt5 C++ GUI Programming Cookbook[END_REF] libraries. It provides a graphical user interface (GUI) to activate and supervise the proposed functionalities of the elastography robotic system. The ViSP library was used to perform all the computation process related to the elastogram estimation, control law, and robot communication. The elastogram display and application GUI were implemented using VTK and QT. The details of the implementation are presented for every case.

Results of the 2D case implementation

The multi-thread software application presented in Figure 2.29 was developed for the implementation of the 2D case. The RF frame is created as a shared pointer 1 , which is shared with three threads (RFtoBMode, Acquisition and Elastography). The RF frame is updated every time the Acquisition thread acquires the RF data from the client-server communication, previously described. Then, the RF frame is converted to b-mode image to be displayed in the main thread containing a GUI. The same RF frame is stored by the Elastography thread such that the first RF frame is defined as the pre-compress and the next one as the post-compress state. After the second frame is sent into the Elastography thread, the incoming RF frames are always defined as post-compress state, and the previous post-compress RF frame is shifted to the pre-compress state. In this thread the elastogram is computed after a ROI is selected by the user in the GUI. Then, the elastogram in the ROI is overlaid in the b-mode image. The robot is activated in the GUI to perform the palpation motion, and the automatic centering of the stiffest tissue in the center of the FOV of the US probe is also enabled through the GUI. The required centroid for this task is computed in the main thread using the current elastogram. The process time of the elastography algorithm was 20 ms corresponding to 50 FPS, over a ROI of 50% of the RF frame size. Therefore, it is compatible with the time constraint of a real-time robotic control scheme. 1 Address of the shared memory that allows data exchange between the different threads 67 2.6. EXPERIMENTAL RESULTS

Experimental results using a training abdominal phantom

For the next experiments we used the setup presented in section 2.1 with the 3D convex ultrasound probe in 2D imaging mode. The force control law was performed with a frequency higher than the other tasks (200 Hz). The visual control was performed with the same period as the image capture (24 fps). The experiments were performed over the ABDFAN ultrasound examination training model (see section 2.1) simulating the abdomen of a patient. The phantom manufacturer specifications note the presence of lesions and tumors. Experiments were performed selecting a ROI including hepatic lesions and pancreatic tumors.

In the experiments, the probe was initially positioned above the phantom, without contact, and oriented with an initial angle θ i ni t and initial force F 0 = 0. Then, to demonstrate the efficiency of the general control law (2.60), we set a desired sinusoidal force signal with F 0 = 5 N and ∆ F = 2 N. We selected ∆ F = 2 N, because it is the minimum force variation required to compute an elastogram, which was found by applying different forces in the finite element analysis. The lower relative error compared with the ground truth computed by the FEA was 5.3%, which proves that our elastogram is well estimated. The automatic horizontal centering of the ROI is activated once the user selects this area in the graphical interface. We performed several experiments and present here the results obtained from one of the experiments.

A set of five desired angles for the probe orientation task is considered:

θ 0 = θ i ni t - 10 • , θ 1 = θ 0 +5 • , θ 2 = θ 1 +5 • , θ 3 = θ 2

Elastogram quality improvement

To improve the elastogram quality, we propose to align and average the different elastograms obtained for all probe orientations. This alignment is performed by a warping function that consists in applying a translation based on the centroid relative position and the image relative rotation between each elastogram of the object of interest (blue region in the images of the third row in the Figure 2.31, where dark blue is the lower strain, and dark red is the highest strain). Once we obtain the warped elastograms Further, in elastography, the concept of contrast-noise-ratio (CNR e ) is expressed as,

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o)
CNR e = 2(µ s -µ b ) σ 2 s + σ 2 b (2.63)
where µ s and µ b represent the mean value of the stiff and background tissues, and σ 2 s , σ 2 b denote the strain variance, respectively. The value of the CNR e allows to make a decision on accepting or rejecting the presence of a lesion as presented in [START_REF] Varghese | An analysis of elastographic contrast-to-noise ratio[END_REF], and a higher level of CNR e suggests better ability to detect the lesion. Therefore, to evaluate the elastogram quality in our experiment, we compute the CNR e in the strain images for each probe orientation and for their mean as shown in Table 2.4. The highest CNR e is obtained for the image of the mean of the elastograms as expected.

Results of the 3D case implementation

For the implementation of the 3D process, we developed a multi-thread software application as described in Figure 2.34. A shared pointer related to the memory of the RF volume is continuously updated by the acquisition thread (frame by frame). This shared pointer is read by the RFtoBMode thread (process in charge of converting the RF volume to b-mode volume) and the Elastography thread (process to compute the 3D elastogram in a VOI) once a volume is completed. The Display object, in the main thread, contains the functions to display three orthogonal planes (sagittal, axial, coronal) of the volume (see Figure 2.33) using VTK library (Visualization Toolkit [START_REF] Schroeder | Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition[END_REF]). This object also allows the user to select the VOI by displacing the planes to the desired position. We perform a set of experiments on a homemade gelatin phantom containing two duck gizzards (see Figure 2.35), and we present the evolution of one experiment, which is the base for every experiment performed.

We set the values of F mi n = 1.5 N, F max = 2.5 N and λ f = 0.002. F mi n and F max were estimated empirically for the gelatin phantom, and they can be adapted for any other kinds of tissue. In this implementation, for the second and third tasks, we used adaptive gains (see Equation (1.42)). The parameters in the adaptive gains in the control law of Equation 2.60 are set as λ t (0

) = 0.1, λ t (∞) = 0.03 , λt (0) = 0.3, λ θ (0) = 1.5, λ θ (∞) = 0.2 and λθ (0) = 2.3.
The experiment begins with an initial probe position where a stiff object of interest is located in the 3D US probe FOV (red point in Figure 2. 35-left). Then, the automatic palpation task with the robot is activated to perform the compression of the tissue. Next, four points delimiting the VOI are chosen using the developed GUI and displayed by small yellow spheres as shown in Figure 2.33. Once the VOI is selected, the 3D elastogram is estimated for every pair of RF volumes. The centroid of the 3D elastogram is computed as we previously described in section 2.3, and it is sent to the automatic centering control task. The probe orientation is always active, and the user can change the desired orientation of the probe any time through the GUI. We show in Figure 2.36 the plots of the evolution of the probe velocities and error of the three tasks for one experiment. We can observe at the beginning of the experiment that the only active velocity is the v y (force control). At time t ≃ 23s the system is paused to select the VOI, and at time t ≃ 72s the process continues. Then, the center of mass of the biggest stiff tissue is computed and the velocities v x and v z applied by the visual

Conclusion

We have proposed a new approach for automatic palpation in this chapter. We based our methodology on three hierarchical tasks. The main task, the compression motion based on the force variation, is required every moment for the elastography system. Then, the secondary task, ultrasound probe centering with a target tissue, performs the visual servoing-based approach to automatically center the probe's FOV with the biggest stiff tissue in a ROI or VOI. The third task, the probe orientation, is used to explore the tissue targeted with the secondary task from different orientations of the probe. In the 2D case, our system was also used to improve the quality of the elastogram. It is based on the elastograms captured at different probe orientations using the third task. This third task will be revisited in chapter 4, where the probe orientation will be controlled with a haptic interface. The control system of the 2D case was designed to control the 3-DOF of the ultrasound probe corresponding to the in-plane motions. We have also presented an extended approach for the use of a 3D US probe allowing our system to perform realtime 3D quantitative ultrasound elastography. It is based on a control scheme similar to the 2D case but controlling the 6-DOF of a 3D motorized ultrasound probe in order to consider the out-of-plane motions. The experimental results demonstrated the feasibility of the proposed concept. The experiments were performed on static tissues which is not the case with the human body. However, next chapter will present an extended dense approach to deal with the physiological motions and also with the control of the out-of-plane motions with a 2D probe.

CHAPTER 2. AUTOMATIC PALPATION FOR ULTRASOUND ELASTOGRAPHY

RELATED WORKS

to the tracking of deformable tissues. Then, the general principle of the visual dense tracking approach is presented in Section 3.2. In the same section, deformable tracking models and their relation with the strain map are presented. Afterwards, in Section 3.3, the motion compensation of the ultrasound probe is introduced. A dense visual servoing based on ultrasound b-mode images is elaborated to compensate the in-plane and outof-plane motions of a targeted tissue of interest. Experimental results obtained on a moving abdominal phantom are presented and discussed in Section 3.4.

Related works

Image tracking has gained interest over the past 20 years in the medical imageprocessing field. Image tracking is the process of aligning two or more images. Such images can originate from a single imaging modality or from different modalities; they can be taken from different patients to study the same organ, tissue or structure; or they can be obtained from an acquisition through time, where temporal structural changes are analyzed. Image tracking can extract valuable information that can be spread on two or more images. Defining the transformation model that best aligns the structures or tissues of interest present in the images is of the utmost importance. Deformable transformations are capable of managing significant changes of biological structures.

Accordingly, deformable image tracking is a fundamental task in medical image processing. Currently, there is a wide variety of techniques developed for medical imaging tracking as presented in [START_REF] Alam | Evaluation of medical image registration techniques based on nature and domain of the transformation[END_REF]. We briefly describe few works related to the tracking method of deformable tissues that will be presented in the next section.

In heart surgery, the surface of a beating heart was tracked in stereoscopic images using a thin-plate spline (TPS) deformation model [START_REF] Richa | Three-dimensional motion tracking for beating heart surgery using a thin-plate spline deformable model[END_REF] (see Figure 3.1). Similar methodology was presented to track the deformation applied to a soft tissue CHAPTER 3. ROBUST MOTION COMPENSATION phantom from 3D ultrasound images [START_REF] Lee | Intensity-based visual servoing for non-rigid motion compensation of soft tissue structures due to physiological motion using 4D ultrasound[END_REF]. In this work, a 3D non-rigid tracking algorithm based on a thin-plate spline 3D deformation model was considered and a visual servoing scheme was designed to automatically move the ultrasound probe to compensate the rigid motion components of the tissue. Recently, an approach based on a massspring-damper model combined with the dense information contained in a sequence of ultrasound 3D images was proposed to track the deformation of liver tissues [START_REF] Royer | Realtime target tracking of soft tissues in 3d ultrasound images based on robust visual information and mechanical simulation[END_REF].

However, none of these works were applied to obtain the strain information of a moving tissue which is the goal of our work presented in this chapter.

Dense visual tracking

In Chapter 2, a geometrical-feature-based visual tracking approach was presented, where the centroid of the stiffest tissue in a selected ROI was extracted from the elastogram using a segmentation algorithm. Unlike this previous approach, we propose here to use directly the appearance of the b-mode image to perform a non-rigid visual tracking of the ROI containing the deformable tissue of interest. In addition, the elastogram is estimated exploiting the output of the visual tracking process. The use of intensities, colors and textures instead of geometrical features has been proposed in several works. Visual tracking based on image pixel intensities was introduced by Lucas and Kanade [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] as an approach for image matching (registration) in stereo vision. Through the years, image registration has been extended in several approaches. For example in [START_REF] Irani | A unified approach to moving object detection in 2D and 3D scenes[END_REF], a technique for moving object detection was presented. A two-view approach for moving objects detection was introduced in [START_REF] Schindler | Two-view multibody structure-and-motion with outliers through model selection[END_REF]. The common factor of these works that are referenced as dense registration techniques is the use of image templates as visual features.

One of these approaches is the photometric visual tracking technique that considers the intensity of all pixels in the image registration process [START_REF] Hager | Efficient region tracking with parametric models of geometry and illumination[END_REF]. However, these dense registration approaches are computationally expensive due to the large quantity of pixels to match. Nevertheless, an efficient optimization for direct image registration has been presented in [START_REF] Benhimane | Real-time image-based tracking of planes using efficient second-order minimization[END_REF][START_REF] Dame | Second order optimization of mutual information for real-time image registration[END_REF]. Then, the process to match the pixel position of I c with I t is defined as:

p = arg min p N -1 k=0 E (I t (x k ), I c (w (x k , p))) (3.1)
where w (x, p) warps the image point coordinates x using the transformation parameters p and E is a similarity function.

Similarity measures

The similarity function E is the cost function to be optimized by finding the parameters p and it is an essential part of the visual tracking process. E measures the similarity or dissimilarity between two images. We present two similarity metrics that we consider in this thesis.

Sum of squared difference (SSD)

The sum of squared differences (SSD) measures the difference between the pixel intensities of a template image I t and a current image I c as,

SSD(I c , I t ) = N -1 k=0 I t (x k ) -I c (w (x k , p)) 2 (3.2)
This function is simple and computationally efficient, which makes it a widely used function in image registration [START_REF] Shi | Good features to track[END_REF]. However, SSD function lacks of robustness to intensity changes and occlusions in the image. These are common occurrences in ultrasound images when the probe is moving as illustrated in Figure 3.2.

Sum of conditional variance (SCV)

To deal with issues as the ones presented in Figure 3.2, a more robust similarity measure is required. In our work we propose to use the sum of conditional variance (SCV) [START_REF] Richa | Visual tracking using the sum of conditional variance[END_REF],

which is robust to global illumination changes and that is expressed as follows:

SC V (p) = N -1 k=0 Ît (x k ) -I c (w (x k , p)) 2 (3.3)
where Ît is the image intensity adaptation of the template relative to the image intensity conditions in the current warped image I c (w (x, p)). The image template adaptation is performed through the expectation operator E as Ît = E(I c (w (x, p))|I t (x)). Thus, the adaptation of every gray level for the reference I t is,

Ît (x) = L-1 i =0 i p I t I c (i , j ) p I t ( j ) (3.4)
where L is the maximum gray level of the template image I t and current image I c . p I t is the probability density function of I t and p I t I c is the joint probability density function of I t and I c . These functions are computed as follows:

p I t I c (i , j ) = p I t I c (I c (w (x, p)) = i , I t (x) = j ) (3.5) = 1 N N -1 k=0 δ(I c (w (x k , p)) -i )δ(I t (x k ) -j ) p I t (i ) = L-1 j =0 p I t I c (i , j ) (3.6)
where δ is a Dirac delta function such as δ(u) = 1 ⇔ u = 0.

Optimization

The selection of the SCV as a similarity function replaces E in Equation (3.1) as follows:

p = arg min p N -1 k=0 Ît (x k ) -I c (w (x k , p)) 2 (3.7)
The variation of the intensity values of Ît (x) with respect to the coordinates x is nonlinear. Therefore, the Equation (3.7) is a nonlinear optimization, and we can solve it using some iterative strategies. Nonlinear optimization iteratively updates the values of the parameters p until convergence. A good performance of the optimization depends on

DENSE VISUAL TRACKING

Strategy Warp parameters increment Update rule

Forward additional [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] ∆p l = arg min

∆p l N -1 k=0 I c l (w (x k , p l + ∆p l )) -Ît (x k ) 2 p l +1 = p l + ∆p l Direct com- positional ∆p l = arg min ∆p l N -1 k=0 I c l (w (w (x k , ∆p l ), p l )) -Ît (x k ) 2 w (x, p l +1 ) = w (w (x, ∆p l ), p l ) Inverse com- positional ∆p l = arg min ∆p l N -1 k=0 I t (w (x k , ∆p l )) -Îc l (w (x k , p l )) 2 w (x, p l +1 ) = w (w -1 (x, ∆p l ), p l )
Table 3.1: Nonlinear optimization strategies. ∆p l is the increment of the parameters at the iteration l . the initialization of parameters. Table 3.1 presents three of the most used strategies to solve the Equation (3.7). Table 3.2 summarizes the advantages and drawbacks of the different optimization strategies presented in Table 3.1. The iterative process of every strategy ends when SC V is minimum or when l (iteration index) has reached N i t (maximum iteration number). In our work, we select the inverse compositional approach due to the advantages of efficiency with respect to the other strategies. This strategy helps us to converge in less number of iterations. Since the nonlinear optimization is solved using an increment CHAPTER 3. ROBUST MOTION COMPENSATION ∆p l , we perform the first order Taylor expansion of I t (w (x k , ∆p l )) as:

Strategy

I t (w (x k , ∆p l )) = I t (w (x k , 0)) + ∇I t ∂w ∂∆p l ∆p l , (3.8)
where ∇I t ∈ R 1×2 is the image gradient of I t defined as ∇I t = ∂I t ∂x ∂I t ∂y . We assume that w (x k , 0) is the identity warp [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF], such that w (x k , 0) = x k . Therefore, the SC V (∆p l ) can be defined as:

SC V (∆p l ) = N -1 k=0 I t (x k ) + ∇I t ∂w ∂∆p l ∆p l -Îc l (w (x k , p l )) 2 (3.9)
The goal of every l -iteration is to minimize the value of SC V (∆p l ), which can be achieved by nullifying the gradient of SC V with respect to ∆p l as:

∂SC V (∆p l ) ∂∆p l = 2 N -1 k=0 ∇I t ∂w ∂∆p l ⊤ I t (x k ) + ∇I t ∂w ∂∆p l ∆p l -Îc l (w (x k , p l )) = 0 (3.10)
We can obtain the change of the parameters for every iteration from Equation (3.10) as follows,

∆p l = -J(∆p l ) + Îc l (w (x, p l )) -I t (x) , (3.11)
where J(∆p l ) ∈ R N ×N p , defined in Equation (3.12), is the Jacobian matrix with N p as the number of parameters in the warp function w . Îc l (w (x, p l )) and I t (x) are row vectors of N -elements containing every pixel in the images Îc l and I t , respectively.

J(∆p

l ) = ∇I t ∂w ∂∆p l (3.12)
We have defined how to solve iteratively a dense visual tracking problem as presented in the Equation (3.7). This solution is invariant to global illumination changes due to the robustness of the SCV similarity metric. Next, we define the warp function w and how to select it according to the complexity of the image registration problem.

Warp transformation

The warp function

R 2 → R 2 : x ′ = w (x,

Rigid transformations

Geometric rigid transformations preserve the distance between points. These transformations are based on Euclidean geometry and include rotations and translations.

A rigid transformation in 2D can be defined as a homogeneous transformation matrix

T ∈ SE (2) 1 : T(ψ, t x , t y ) = R 2×2 (ψ) t 2×1 (t x , t y ) 0 1×2 1 (3.13)
where R ∈ SO(2) and t ∈ R2 are the rotation matrix (2 × 2 size) and the translation vector (2 × 1 size) respectively. These two elements of the rigid transformation matrix are formally defined as:

R(ψ) = cos(ψ) -sin(ψ) sin(ψ) cos(ψ) (3.14) t(t x , t y ) = t x t y (3.15)
The rigid transformation, defined in Equation (3.13), is applied to any pixel with coordinates x = [x y] ⊤ to obtain the new pixel coordinates:

x ′ 1 = T(ψ, t x , t y ) x 1 (3.16)
This leads us to obtain the warp function w to compute the transformed coordinates

x ′ = [x ′ y ′ ] ⊤ as, x ′ = w (x, p) = x cos(ψ) -y sin(ψ) + t x
x sin(ψ) + y cos(ψ) + t y (3.17) where p = [ψ t x t y ] ⊤ are the parameters of the warp function.

CHAPTER 3. ROBUST MOTION COMPENSATION

Once we have defined the function w , we can estimate the Jacobian defined in Equation (3.12) as:

J(∆p l ) = ∇I t -x sin(ψ) -y cos(ψ) 1 0
x cos(ψ) -y sin(ψ) 0 1 (3.18) where J(∆p l ) ∈ R N ×3 is the Jacobian that estimates the parameters variation in the iterative process to solve the image registration problem (Equation (3.11)) with a rigid transformation.

In order to evaluate the performance of every transformation model presented in this chapter, we implemented the visual tracking algorithm in C++ in a Linux notebook (Intel i7 CPU @2.1 GHz). For this evaluation, the convergence conditions are the maximum number of iterations N i t set to 50 or the minimum SC V value set to 1×10 -6 . Then, we acquired a sequence of 500 ultrasound b-mode images during the application of the palpation motion task, presented in Section 2.2, and the presence of lateral in-plane motion also applied with the ultrasound probe. Due to palpation motion, deformations are produced along the image sequence. We tested the rigid transformation with this image sequence using a ROI as shown in Figure 3.4a delineated in green color. This ROI was tracked through the image sequence until the last image as shown in Figure 3.4. The performance of the dense rigid tracking was evaluated with the image absolute error,

e abs = e ⊤ e N , (3.19)
where e = Îc (w (x, p)) -I t (x). The algorithm converged after 29 iterations and the value of e abs reached 9.70 × 10 -7 . It is well known that the palpation motion introduces deformations that are non-rigid, therefore, it is likely that a rigid transformation will not be sufficient for this visual tracking problem.

Non-rigid transformations

The compressions required in the elastography generate deformations in the ultrasound images which cannot be approximated using rigid transformations. Therefore, it is necessary to introduce non-rigid transformations as warp functions. These transformations are commonly used when the image template is distorted or deformed. Moreover, nonrigid transformations are classified as linear and non-linear. Linear transformations are used for image distortion, which can come from the change in the normal of a planar object in an image sequence. Image distortions also appear when a planar object is viewed from a different point of view in relation to the camera perspective. The most wellknown linear transformations are the affine and the projective transformations. However, the projective transformation is not suitable for the ultrasound images, since the image is not reconstructed from a perspective model as images from a camera. Nonlinear transformations are used when deformation is applied to the image template in an image sequence. For example, in an image sequence where an elastic material is deformed by an external force, a non-linear transformation should be defined to track the image template through the sequence.

Next, we describe the most common linear transformation, the affine transform. After, we detail the free-form deformation (FFD) and the TPS as non-linear transformations for image deformation. For every transformation, the warping function and the Jacobian (required for the visual tracking process) are formulated. In order to present a comparative of the image registration process, we use the same ultrasound image sequence presented for the rigid tracking case. This image sequence was acquired with our experimental setup by applying deformation on the phantom with the ultrasound probe.

Affine

The main characteristic of an affine transformation is that it preserves parallels lines in the image after being applied. This transformation combines the rigid transformation motions, along with scale and shear in a set of six parameters p. Basically, four parameters of p modify rotation, scale and shear of a pixel coordinate x = [x y] ⊤ . The remaining two parameters are directly related to the translation of x. An affine transformation can be define as a matrix as follows,

A(p) = p 0 p 1 p 2 p 3 p 4 p 5 (3.20)
where p = [p 0 p 1 p 2 p 3 p 4 p 5 ] ⊤ is the vector containing the parameters of the affine transformation with N p = 6.

The affine transformation A is applied to any pixel coordinate x to obtain the transformed coordinate x ′ . If we use the augmented version of A expressed in homogeneous coordinates, we can then compute the transformed coordinate as:

x ′ 1 =    p 0 p 1 p 2 p 3 p 4 p 5 0 0 1    x 1 . (3.21)
The warp function w that maps x → x ′ using the affine transformation is then expressed as,

x ′ = w (x, p) = xp 0 + y p 1 + p 2 xp 3 + y p 4 + p 5 (3.22)
in order to be adapted to the image registration process. The resulted Jacobian of the Equation (3.12) can be then expressed as:

J(∆p l ) = ∇I t x y 1 0 0 0 0 0 0 x y 1 (3.23)
In the dense affine registration process, the Jacobian J(∆p l ) ∈ R N ×6 allows to estimate the parameters variation in the iterative process to solve the image registration (Equation (3.11)).

We show in Figure 3.5 an example of the performance of the dense registration process using the affine model. For comparative purposes of the warp functions performance, we use the same image sequence as the one used for rigid registration. In this case, the image absolute error after 31 iterations was 9.23 × 10 -7 which is better than the error obtained with the rigid transformation.

Free-Form deformation

Free-Form deformation (FFD) is a common technique in computer graphics and animation design. The main concept of FFD relies on the use of hierarchical transformations to deform an object [START_REF] Barr | Global and local deformations of solid primitives[END_REF]. The transformations include twisting, bending, tapering and stretching of the object. A most generalized approach for FFD was presented in [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF] allowing to apply global and local deformations to surfaces of any degree (e.g., plane, quadric or parametric). This generalization was performed by using Bernstein polynomials to design the spline functions. In the context of image registration, FFD has been used to match breast MRI images where deformations are present [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast mr images[END_REF]. More recently, image matching with FFD has been also improved in terms of computational cost [START_REF] Brunet | Feature-driven direct non-rigid image registration[END_REF] using an intuitive feature-driven framework.

FFD is one of the most common transformation models in medical imaging, where a rectangular grid of N p = N p x × N p y control points are placed on the template image. The displacement of the control points deforms the image using the products of univariate splines. The deformation of the image using FFD is obtained by applying the warping function:

w (x, p) = N p y j =1 N px i =1 p k B i (x)B j (y), (3.24)
where p k ∈ R 1×2 is the k-th control point with index number k = ( j -1)N p x + i and B i is the basis function of the cubic B-splines:

B i (x) =                  B 1 (x) = x3 6 if x ∈ [k i , k i +1 ] B 2 (x) = -3 x3 +3 x2 +3 x+1 6 if x ∈ [k i +1 , k i +2 ] B 3 (x) = 3 x3 -6 x2 +4 6 if x ∈ [k i +2 , k i +3 ] B 4 (x) = -x3 +3 x2 -3 x+1 6 if x ∈ [k i +3 , k i +4 ] 0 otherwise (3.25) with x = x-k l δ |x ∈ [k l , k l +1 ], δ =∥ k l +1 -k l ∥ . k l is the l-th control point.
A generalization of the Equation (3.24) can be expressed as

w (x, p) = ψ ⊤ P, (3.26)
where ψ ∈ R N p and P ∈ R N p ×2 are the vectors of the basis functions and control points respectively:

ψ ⊤ = [B 1 (x)B 1 (y) . . . B N px (x)B 0 (y) . . . B N px (x)B N p y (y)] (3.27)
P = [p 1 . . . p N p ] (3.28)
One big advantage of the Equation (3.26) is that the vector of the basis functions can be precomputed for the pixel coordinate x. This helps us to improve the computational cost of this warping function. Therefore, the deformation of the current image depends on the variation of the parameters p. Moreover, the computational cost also increases with a large number of control points.

As in previous transformations, we define the Jacobian of the warp function required for the image registration process. This one is expressed as,

J F F D = ∂w ∂∆p = ψ ⊤ 0 1×N p 0 1×N p ψ ⊤ (3.29)
The Jacobian J F F D ∈ R 2×2N p helps us to solve iteratively the image registration process of the Equation (3.11) using the FFD warp function. We show in Figure 3.6 the performance of the dense visual tracking system using FFD transformation applied in the same image sequence as the previous transformations. In this case a 5×5 grid of control points was used. The image absolute error after 50 iterations was 4.69×10 -7 , which is a better result than the previous transformation models. This improvement is due to the FFD transformation fitting better the image template to the deformation presented in the image.

DENSE VISUAL TRACKING

Thin-plate splines

This nonlinear transformation was suggested for image registration [START_REF] Goshtasby | Registration of images with geometric distortions[END_REF], and it is based on the analogy of how a metallic thin plate is deflected by normal forces at discrete points. Thin-plate splines (TPS) have also been used for dense image tracking in [START_REF] Delabarre | Dense non-rigid visual tracking with a robust similarity function[END_REF] and [START_REF] Brunet | Feature-driven direct non-rigid image registration[END_REF]. TPS warping function is a combination of an affine transformation and deformation parameters (control points) as follows:

w (x, p) = a 0 a 1 a 3 a 4 x + a 2 a 5 + N c -1 k=0 κ k x κ k y φ(d (x, c k )) (3.30)
where N c is the number of control points c, κ k x and κ k y are the weights of each k control point along the x and y axes respectively. These weights represent the force amplitude applied at the control point position. φ is the thin-plate kernel defined as,

φ(x) = x 2 l og (x) 2 and (3.31) d (x
, y) is the euclidean distance between the points x and y.

The parameter vector of the warping function p of dimension 2N c +6 is expressed as:

p ⊤ = a 0 a 1 a 2 a 3 a 4 a 5 κ ⊤ x κ ⊤ y (3.32)
where the first six parameters are the parameters of the affine transformation. κ x and κ y are vectors (N c elements) containing the weights κ k x and κ k y respectively:

κ ⊤ x = κ 0 x . . . κ N c -1 x (3.33) κ ⊤ y = κ 0 y . . . κ N c -1 y (3.34)
In the registration process, the positions of the control points c are initially distributed in an equidistant grid inside the image. Then, the values of the forces applied at every control point are changed to match a deformed image with a template image. We illustrate in Figure 3.7 how the TPS warp function can be adapted to deform the template image.

Image registration with TPS is commonly performed by changing the parameters p to adjust the current image, in an image sequence, with the image template. The optimization requires the Jacobian of w , as we shown in the previous transformations, which can be obtained as: ∂w ∂∆p = J A J κ (3.35) where:

J A = x y 1 0 0 0 0 0 0 x y 1 ∈ R 2×6 , J κ = φ 0 • • • 0 0 • • • 0 φ ∈ R 2×2N c , φ = φ(d (x, c 0 )) • • • φ(d (x, c N c -1 )) ,
x and y are the pixel coordinates in the image I t .

The Jacobian J(∆p l ) ∈ R N ×(6+2N c ) can be obtained through Equation (3.12). This Jacobian is employed to solve the image registration system of the Equation (3.11) using the TPS transformation. We show in Figure 3.8 the performance of this process using the same example as the one tested with the previous transformations. This figure shows the results of the dense tracking process using the TPS transformation with 5×5 control points. The image absolute error after 21 iterations was 1.33 × 10 -7 . This result is better than the result obtained by the FFD transformation. In addition, TPS reaches the convergence in the image registration process in less number of iterations than the FFD image registration.

The performance of the registration process using the different warp functions provides a perspective to select the best function for our application. Due to the mechanism of the elastography using the compression of the tissue, the ultrasound image tends to present deformation. The non-linear warp functions perform better than the other functions under deformations. Therefore, from our performance comparison study of the different tested approaches, summarized in Table 3.3, we propose to choose the TPS as the warp function in our visual tracking system due to its faster convergence and minimum absolute error. The computational time is also a major factor considered in our decision, since our system requires real-time visual tracking capability. Table 3.3: Performance evaluation of the transformations in the visual tracking system.

Strain estimation based on optical flow

Strain estimation is a process that depends on the motion estimation of the elements contained in the ROI to generate an elastogram. From the TPS registration, we can obtain the displacement maps U (x, y) and V (x, y) (lateral and axial directions respectively).

Let us define x ′ = w (x, p) as the corresponding coordinates of x after the tissue deformation. For every x we have a displacement vector D(x) = x ′ -x. Then, we can obtain U (x) and V (x) as the lateral and axial components of the displacement vector D(x) (see Figure 3.9).

Assuming that we have at least a grid of 3 × 3 control points, then the elastogram ε(x, y) can be computed as in Section 2.2.2.2 by convolving a least-squared (LSQ) strain filter with the axial displacement map V (x, y). Therefore, the elastogram ε(x, y) can be computed using the information from the deformable registration.

Motion compensation

We have developed a method to track a template image (in our case the ROI) when the tissue is deformed. Since the elastography process requires the compression of the tissue, deformable registration is essential to track the ROI. However, some of the physiological motions cause the ROI to go outside of the image plane when using a 2D probe.

This leads to a failure of the image tracking. Therefore, to solve this issue, we propose to use a control system that uses prior information of parallel images planes to the plane containing the ROI. This control system helps us to fully control the 6-DOF of the ultrasound probe in order to always maintain the visibility of the ROI. Preserving the position of the ROI stable even when the tissue is moving is essential for the right estimation of an elastogram. Moving tissue can cause motions not only inplane but also out-of-plane as emphasized in [START_REF] Krupa | Real-time tissue tracking with b-mode ultrasound using speckle and visual servoing[END_REF] (see Figure 3.11). This tissue motion can be compensated by controlling the motion of a 2D US probe using the intensitybased ultrasound visual servoing presented in [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF][START_REF] Nadeau | Intensity-based ultrasound visual servoing: Modeling and validation with 2-D and 3-D probes[END_REF].

6-DOF motion compensation by dense visual servoing

We propose to use a similar visual servoing approach to the one presented in [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF] to control a robotic arm holding the probe in order to automatically compensate the relative non-axial motions between the probe and the moving tissue of interest to analyze. Nonaxial motions correspond to one lateral and one rotational motions in the US plane and one lateral and two rotational motions out of the US plane (see Figure 3.11).

In [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF], a visual servoing method that uses the intensity information of the pixels inside a ROI has demonstrated the feasibility to control the 6 DOF of a 2D ultrasound probe for compensating both in-plane and out-of-plane rigid motions. However, in this previous work the tissue was assumed rigid without considering deformation due to internal physiological motion or the presence of mechanical compression. To deal with the soft tissue deformations, we propose to improve the method of [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF] by using the non-rigid motion estimation algorithm presented in Section 3.2.5.3

We briefly recall the principle of the ultrasound dense visual servoing approach [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF].

The aim is to control the probe velocity v expressed in the frame F cp (see Figure 2.5).

The visual features vector, s, used in this control scheme contains directly the intensities of the pixels inside a ROI such as:

s I = (I 1,1 , . . . , I M ,N ) (3.36)
where I u,v is the intensity in gray level for the 2D pixel coordinates (u, v) in the US image.

The interaction matrix L I u,v ∈ R 1×6 that relates the variation of the pixel intensity to the probe velocity v, such that I u,v = L I u,v v, is given by:

L I u,v = ∇I x ∇I y ∇I z y∇I z -x∇I z (x∇I y -y∇I x ) (3.37)
where ∇I u,v = [∇I x ∇I y ∇I z ] corresponds to the 3D image gradient associated with the pixel (u, v). The three components,

∇I x = ∂I u,v ∂x , ∇I y = ∂I u,v ∂y and ∇I z = ∂I u,v
∂z are obtained with 3D derivative filters (see Figure 3.12), as performed in [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF], applied to a thin volume composed of 5 parallel slices captured by moving the probe during an initial procedure before launching the visual servoing. The values of x and y are the metric coordinates of the pixel (u, v) in the image obtained from the intrinsic parameters of the probe:

x y = s x (u -u cp ) s y (v -v cp ) (3.38)
as e s = s w I -s * I , and we establish the desired visual error variation as ė * s = -λ s e s with λ s being the visual control gain. Unlike [START_REF] Nadeau | Intensity-based direct visual servoing of an ultrasound probe[END_REF], here s w I = I c (w (x, p)) is the current image warped with the TPS warping function using the current parameters p. This major improvement allows the visual servoing approach to be robust to the presence of the non-rigid motion induced by the tissue deformation. The desired pixel intensities vector corresponds to the intensities value of the pixels contained in the ROI of the initial image I t , s * I = (I t 1,1 , . . . , I t m,n ). Then, the control law applied to the probe for performing the automatic motion compensation is provided by:

v s = L + s ė * s .
(3.40)

Control fusion

In order to fuse the automatic motion compensation by visual servoing and the force control (presented in Chapter 2 Section 2.2.1), we can define a control law for the probe velocity v using the redundancy control framework (presented in Chapter 2 Section 2.5).

We set the force control law as the highest priority task, remaining as in Equation (2.6).

Then, the secondary task that corresponds to the visual servoing can be expressed as:

v s = (L s P f ) + (ė * s -L s v f ) (3.41)
where P f = I-L + f L f is the projector operator onto the null space of L f . I is the identity matrix of size 6.

Finally, the general control law that allows to control the 6-DOF of the 2D probe is given from (2.6) and (3.41) as:

v = v f + v s (3.42)
This control fusion allows to control the 6-DOF of the ultrasound probe in order to automatically compensate the motions and keep the ROI always visible. In the next section we present the experimental results of the proposed approach that makes possible the estimation of the elastogram of a moving tissue.

Experimental results

We present the results obtained with the same experimental setup proposed in Chapter 2. The images from the scanner were sent to the workstation at a rate of 40 FPS

EXPERIMENTAL RESULTS

(frames per second). Force control was performed at higher frequency (200 Hz). We developed a C++ software with a graphical user interface (GUI), and we used ViSP [START_REF] Marchand | ViSP for visual servoing: a generic software platform with a wide class of robot control skills[END_REF] for the communication with the robot. The experiments were performed on the ABDFAN ultrasound phantom. Now, we describe the complete process for one experiment. Initially, the probe was positioned above the phantom without contact. Then, through the GUI, we enabled the force control without oscillation (F 0 =5 N and ∆ F =0 N), and we can see that the measured force value reaches the desired force in the strip (light gray background) of the plot in the that the cyst is very visible as opposed to the elastogram presented in Figure 3.15b where a perturbation motion was applied to the phantom without activating the automatic compensation by visual servoing. It is clear that in this case the elastogram estimation is perturbed and can not provide any useful information. Figure 3.15c shows the obtained elastogram when the phantom is moving and the automatic motion compensation by visual servoing is activated. This last test demonstrates the efficiency of our approach since the cyst is very visible and similar to the case where the phantom was motionless.

Conclusion

Physiological motions are always present in real tissue making it necessary to consider motion compensation in the design of our robotic-assisted system for elastography. We have encountered several challenges to obtain the elastogram of a moving tissue. First, since the elastogram computation depends on the axial motion estimation between the pre-and post-compression states, a large lateral motion causes wrong measurements in the elasticity of the tissue. The process presented in Chapter 2 which includes a block matching algorithm (BMA) and an optical flow (OF) algorithm for motion estimation can deal with small lateral motions introduced to the tissue. However, it does not consider large motion perturbations in the tissue when using a 2D probe, neither it can deal with out-of-plane motions. On the other hand, the system using a 3D probe, also presented in Chapter 2, takes into account in-and out-of-plane motions in the robotic centering task. However, the slow acquisition rate of the volumetric information makes the system fail when fast motions are induced to the tissue. These reasons lead us to the design of a robust approach for motion compensation.

In this chapter, we use the 2D ultrasound probe to reach a faster acquisition rate than with the 3D probe, in order to be reactive to fast motion introduced to the tissue.

Contrary to the approach presented in Chapter 2, we consider the b-mode image for our visual tracking process instead of the RF signals. This makes the process two times faster, since the visual tracking does not require two acquisitions of RF data. The visual tracking system presented here was tested with several dense approaches considering an image template as our ROI to estimate the elastogram. Rigid and non-rigid transformations were implemented to evaluate the performance of our visual tracking system.

This also demonstrated that the sum of conditional variance (SCV) is robust to intensity changes usually observed in ultrasound images. The thin-plate splines (TPS) transformation was selected due to its fast performance with respect to the other transformations evaluated (rotation-translation, affine and free-form deformation transformations). TPS transformation uses control points placed inside of the ROI and then computes the displacement of the points due to the deformation of the tissue reflected in the ultrasound image. The use of the control point displacements that are estimated by the proposed non-rigid dense visual tracking allows us to avoid the axial motion estimation from the pre-and post-compressed RF signal. However, the palpation motion task is still required to generate a slight deformation of the tissue along the axial direction. In addition, our approach uses the b-mode image ROI in a dense visual servoing approach that automatically moves the ultrasound probe to compensate the in-plane and out-ofplane motions of a moving tissue using the 2D ultrasound probe. This process considers the information of five parallel images captured at the initial position of the probe in order to estimate the 3D image gradient that is needed in the control law for compensating both the in-plane and out-of-plane tissue motions. Preliminary ex-vivo results have demonstrated the feasibility to estimate the strain map of a moving tissue.

CHAPTER 4 TELEOPERATION AND HAPTIC FEEDBACK BASED ON ELASTOGRAM

The term teleoperation refers to a process performed at a distance. It comprises a robotic system where a human operator controls a remote robot. Teleoperation has been used in applications where close manual operation is hazardous or where access is limited (e.g., nuclear waste manipulation, underwater exploration). The first teleoperation system was designed to handle nuclear material in 1940 [START_REF] Vertut | Teleoperation and robotics[END_REF]. Since the 1990s the use of teleoperation for medical purposes began to appear along the concept of computed-assisted surgery. The first teleoperated system related to this field was the ZEUS surgical robot (Computer Motion, Inc.) developed in 1995 [START_REF] Kumar | Telesurgery[END_REF]. This system comprises three robotic arms mounted on a table, one holding an endoscope which provides a view of the internal operating field, and the others holding surgical instruments. The robotic arms are controlled by the surgeon through a console. Currently, ZEUS is discontinued and the da Vinci surgical system (Intuitive Surgical, Inc.) is the most widely used robotic system for telesurgery in hospitals. The da Vinci surgical system provides to the surgeon a console that renders force feedback and 3D vision of the internal operating field. The current version of the da Vinci system is equipped with four robotic arms with one holding the camera and the others actuating the surgical instruments. Force feedback is provided to the surgeon via the joysticks of the console when the instruments are in contact with the tissues. Usually, the force feedback is associated with the name of haptic feedback and its rendering is performed at a higher frequency than the visual feedback in order to provide a responsive interface to the user [START_REF] Cavusoglu | Multirate simulation for high fidelity haptic interaction with deformable objects in virtual environments[END_REF]. Such haptic feedback allows the surgeon to increase his perception of the scene and it is therefore a functionality of great interest to facilitate the execution of the intervention. For the same reason, we propose in this chapter to provide a haptic feedback functionality to our robotic palpation system.

HAPTIC FEEDBACK

We develop hereafter an approach that will provide a force feedback to the examiner from the estimated strain map. Indeed, our objective is to give the examiner the abilities of feeling the rigidity of a tissue while visualizing its elastogram with the use of a haptic device. First, the basic concepts of haptic feedback are introduced in Section 4.1.

In Section 4.2, the design of the system that translates the elastogram into force feedback is presented. Afterwards, in the same section, a method to remotely control the ultrasound probe held by the robot using a haptic device is described. In Section 4.3, results obtained from experiments performed on the abdominal phantom are presented.

Finally, Section 4.4 concludes this chapter.

Haptic feedback

The sense of touch is used to perceive and explore the world, helping us to identify objects or warning us when touching something dangerous. There are two types of force feedback when making physical contact with an object, kinesthetic and tactile. The kinesthetic feedback is the perception of the internal status of the body. On the other hand, tactile feedback is the response that allows us to feel the material or texture of any object. The combination of these two feedbacks facilitates the adjustment in the configuration of the body to interact with the object. For example, the configuration adopted with a hand to hold something lightweight such as a cotton ball is not the same as the configuration to hold a ceramic mug.

In robotics, the haptic feedback is commonly associated with the two types of force feedback used to explore an object by hand. A wide variety of devices have been designed to emulate the kinesthetic and tactile force feedbacks with the aim of feeling virtual objects. Those devices can be easily identified by their structure: the kinesthetic haptic devices are usually grounded while the tactile ones are wearable (see Figure 4.1).

Currently, the most common haptic devices are the kinesthetic kind, which can be classified by their configuration: manipulandum, grasp, and exoskeleton (see Figure 4.2). The manipulandum configuration involves all grounded devices with 3 to 6 DOF. Grasp configuration concerns the devices simulating grasping interaction at the user's hand. The exoskeleton configuration is associated with devices adapted to the user's body, providing forces at the joints.

In this thesis, we use the Virtuose 6D (Haption S.A.) shown in Figure 4.3. This device is a kinesthetic haptic device with a manipulandum configuration (referred on the next sections as haptic device) with 6 DOF. It also has the capability of applying a 6-DOF force feedback to the user, three translational forces (maximum force of 31N) and three The impedance type devices are the most common force feedback devices (see Figure 4.4). The input of the impedance type device is the motion applied by the user to the kinesthetic haptic device and the output is a force feedback. On the other hand, the admittance type devices have as input the force applied by the user to the kinesthetic haptic device and, as output, the user feels motion (see Figure 4.5). Admittance type devices are not as common as the impedance ones. In our case, the Virtuose 6D haptic device can be configured as both types (impedance or admittance), but we adopted the impedance configuration in order to be compatible with most of the similar haptic devices. ELASTOGRAM

In the next section, we present the design of our robotic system to provide haptic force feedback estimated from an elastogram acquired in real-time. The system allows the teleoperation of a 6-DOF robot, actuating an ultrasound probe (as shown in previous chapters) with a haptic device that is also used to provide force feedback rendered by the elastogram.

Haptic force feedback from elastogram

The use of a haptic device allows us to move a virtual probe within a virtual representation of the environment. The location of the virtual probe can be modified by physically manipulating the haptic device, generating a force feedback. To compute the output force, the algorithm transforms the motion into a force value. For example, if the virtual probe is positioned at the middle of an empty container or box, the force feedback would be null. However, if the position of the virtual probe is displaced to the location of any wall of the container, then the force would be estimated based on the strength of its material. This example provides an idea on how to translate the location of the virtual probe into haptic force feedback.

We consider the virtual probe as a virtual region located at the center of the ROI in the ultrasound image. The virtual probe moves according to the haptic device motion, implying that the ROI is also moving. Next, we present the development of a new approach to estimate the force feedback based on the location of the virtual probe. We also present a teleoperation system to position the ultrasound probe according to the haptic's motion. Figure 4.6 shows a short block diagram of the robotic system proposed that will be developed in the following sections. Two operational modes are presented in this block diagram: impedance haptic control and teleoperation control. The haptic control mode uses the elastogram in a ROI of the ultrasound image to estimate the force feedback that will be applied to the haptic device as we will describe in Section 4.2.1. On the other hand, the teleoperation control mode applies the motion introduced to the handler of the haptic device into the ultrasound probe as we will describe in Section 4.2.3.

Force estimation from elastogram

The diagram of the proposed process to estimate the force based on the elastogram is shown in

(i , j ) = g m (i , j ) max(G m ) , (4.1)
where max(G m ) is a constant value representing the maximum component value of G m , and

g m (i , j ) = e -i 2 2σ 2 x + j 2 2σ 2 y 2πσ x σ y , (4.2)
where σ x and σ y are the standard deviation for lateral and axial directions, respectively.

In our case, these values are typically set as σ x = N 4 and σ y = M 4 aiming to obtain a Gaussian distribution inside of a rounded area. The center of the rounded area is located at the center of G m .

Following the diagram of Figure 4.7, the filtering of the elastogram is performed using the following expression

E f = E • G m , (4.3)
where • is the Hadamard product operator and E f is the resulting elastogram after filtering.

Next, the average scalar strain value of E f , µ ε is computed and then used in the process to estimate the displacement ∆x of a virtual spring after being compressed with a force F (see Figure 4.7). Based on the classic definition of mechanical strain, we can obtain the value of the displacement of the spring as,

∆ x = µ ε L (4.4)
where L is the original length of the spring. The force F of the spring according to Hooke's law is defined as,

F = -k∆ x (4.5)
where k is the stiffness value of the spring given by

k = E πr 2 d (4.6)
where E is the Young's modulus of the soft tissue. The Young's modulus of healthy tissue has values between 2-4kPa. d is the thickness of the compressed tissue and r is the The translation from strain to force is then achieved through Equation (4.5). Now, in the following sections we define how the motion applied by the user on the haptic device is used to move the virtual probe which corresponds to the ROI in the ultrasound image.

Afterwards, we detail the process we developed for teleoperating the ultrasound probe.

We also show how moving the virtual probe gives us the force feedback in function of the tissue elasticity by using the approach previously explained.

Impedance force feedback system

This section describes the impedance force feedback scheme we implemented on the haptic device. The relation between the elastogram and the output force has been established in the previous section. Now, we explain how to relate this output force to the user input motion.

In this study, we consider the use of a 2D ultrasound probe instead of a 3D one in order to obtain 2D elastograms in real time. Therefore, the motion of the 2D US probe will be limited to pure in-plane translations (lateral and axial translations).

h0 M h = ( b M h0 ) -1b M h . (4.8)
where the operator -1 represents the inversion of a homogeneous matrix. The value of the homogeneous matrix b M h0 is measured at the initialization and remains constant, and b M h represents the current measure of the handler pose. As previously mentioned, we only need the relative translation, which is extracted from the relative homogeneous matrix h0 M h as h0 t h ∈ R 3 . Therefore, the relative motion of the handler ∆ h in the x-y plane is obtained from the x and y components of the relative translation h0 t h . The value

∆ d = S I R h ∆ h , (4.9)
where I R h ∈ SO( 2) is a 2 × 2 rotation matrix between the x-y plane of F h and the image frame F I . S ∈ R 2×2 is a 2 × 2 diagonal matrix containing the scale values for lateral and axial direction. These values are computed from the calibration of the ultrasound probe to convert pixels to meters (s x and s y ). Therefore, the matrix S is defined as The displacement of the point (u c , y c ) generates a new elastogram which is translated to a force feedback F using the method described in Section 4.2.1 through Equation (4.5). This completes the impedance haptic system that generates force feedback every time the user applies motion to the handler of the haptic device. As we previously explained, this approach is implemented for the 2D FOV of the ultrasound probe.

S = 1 s x 0 0 1 s y . ( 4
In addition, the adaptation of the probe teleoperation in the robot's workspace to explore the tissue while moving the handler of the haptic feedback can provide more freedom to the user. However, the teleoperation of the probe and the force feedback system can not work together since they would need to share the same DOF of the haptic device.

To deal with this issue, we design two operation modes as shown in Figure 4.6, where the user can switch between teleoperation of the ultrasound probe and virtual probe for haptic sensing through the haptic device handler's buttons. Next section presents the teleoperation of the US probe, and the experimental results of teleoperation and haptic feedback.

Robotic teleoperation

In this section, the teleoperation of the ultrasound probe is presented. This process is based on a master-slave system, where the master is a unit manually operated by a user, in our case the haptic device. The slave robot holding the US probe is duplicating the motion applied on the master device. The master-slave communication is essential to perform delicate procedures. Figure 4.12 shows the proposed system for the masterslave teleoperation of the ultrasound probe.

In the system presented in Figure 4.12, the haptic state contains the information of the kinematics of the haptic device. This information provides the current pose of the handler of the haptic device. An initial pose of the handler is represented by the homogeneous matrix b M h0 once the system is launched. Afterwards, the current pose described by b M h is used to compute the relative pose of b P h with respect to the initial one b P h0 as defined in Equation (4.8). The relative pose represented by h0 M h is measured with respect to the frame F h and to duplicate the displacement with the ultrasound probe we need to define it with the same direction as the ultrasound contact frame F cp (see is obtained using the rotation matrix cp R h between the frames F h and F cp : As we want to reach the desired displacement with the ultrasound probe, we define the current pose as r M cp , and its relative pose cp0 M cp with respect to the initial pose r M cp0 as, cp0 M cp = ( r M cp0 ) -1r M cp (4.13)

M ∆ = cp R h 0 3×1 0 1×2 1 h M h0 (4.11)
We can extract the relative probe displacement from cp0 M cp using the same principle as the one used to extract the displacement component of M ∆ . In this case the six displacement elements are defined inside the measured vector s ∆ :

s ∆ = ∆t p x ∆t p y ∆t p z ∆θp x ∆θp y ∆θp z ⊤ (4.14)
The error between the measured and desired displacement is defined as,

e ∆ = s ∆ -s * ∆ , (4.15)
and the desired exponential error variation as,

ė∆ = -λ ∆ e ∆ , (4.16)
where λ ∆ is the gain of the error variation.

The variation in the measured displacement ṡ∆ is related to the ultrasound probe velocity v with the following expression,

ṡ∆ = L ∆ v, (4.17)
where L ∆ = I 6 is the interaction matrix that relates the variation of the measured displacement ṡ∆ with the probe velocity v and I 6 is a 6 × 6 elements identity matrix.

The variation of the measured displacement ṡ∆ is directly related to the variation of the desired exponential error ė∆ . Therefore, if we replace ṡ∆ with ė∆ in the Equation (4.17), then we obtain a relation to compute the velocity v corresponding to the desired error variation:

v = L + ∆ ė∆ (4.18)
The Equation (4.18) is the velocity control law for the 6-DOF at the frame F cp . Since this control law is designed for the full motion of the probe at frame F cp , we need to limit its movement in the y-axis for security reasons, such that the force control in the y-axis has full priority. Indeed, force control is needed for the elastography process, but it also brings safety when combined with the teleoperation task. The force control law The content of the Elastogram shared pointer is changed by the elastography thread at 24 FPS. The elastography process computes the elastogram with the approach presented in Section 2.2.2. Then, the elastogram is converted into force feedback as presented in Section 4.2.1 in the haptic system thread. The haptic system measures the handler current pose with an update rate frequency of 100Hz, and the new relative displacement ∆ d is sent to the elastography process to change the position of the ROI. Inside the haptic system thread, a mechanism can be activated on demand by the user to switch between the impedance force feedback and the teleoperation control mode. This mechanism uses the two buttons located on the handler to switch between the two control modes. When the teleoperation process is activated, the Robot Control thread computes and applies the velocity to the robot by using the Equation (4.21). The Robot Control thread is also in charge of applying the continuous oscillation in the y-axis force needed to obtain the pre-and post-compressed states of the tissues required for the elastogram estimation. The thread RFtoBMode performs the conversion of the RF frame to a B-mode image and sends it to the main thread for displaying purposes. In this implementation, most of the threads are running at different frequencies. However, the synchronization between the threads is achieved by the Qt connection signal-slot mechanism. Now, we present an experiment that consists in two parts. First, the initial state of the system is in teleoperation mode, and the user can explore the tissue by moving the handler of the haptic device. Figure 4.16 shows in the first row different configurations of the haptic device when the user applied manual motion on the handler during the teleoperation control mode. The second row of Figure 4.16 presents the resulting pose of the robot holding the probe and the third row provides the observed ultrasound image for each configuration. Figure 4.17 shows the temporal evolution of the measured force and the control velocities of the robot applied at the contact frame F cp . The experiment begins with the ultrasound probe above the phantom without making contact, as shown in to variate in order to replicate the motion introduced by the user in the handler of the haptic device. After, in Figure 4.17d we highlighted two moments when the user was applying continuous motion to the haptic device (strips in light-blue color) and when the user stops the motion (strips in light-green color). We can observe the fast convergence (∼ 0.5s) of the teleoperation system in both cases. In these plots, one can observe the variation in the velocities and the errors due to the different motions introduced to the handler of the haptic device. The parameters of the force controller F 0 and ∆ F were set to 5N and 3N, respectively. We can observe that the measured force follows correctly the desired oscillation reference, since the force control task that has the highest priority is not disturbed by the teleoperation task.

The second part of the experiment presents the results of the impedance haptic system when the user switches from the teleoperation mode to the haptic force feedback mode. After selecting a ROI where the elastogram is estimated in real time, the user can feel the force computed from the elastogram while moving the handler of the haptic device. We present the results of one experiment where the user moves the virtual probe position, corresponding to the ROI, inside the ultrasound image with the haptic device.

The user feels the force feedback while moving the haptic device. We can observe in

Conclusion

The sense of touch in palpation examination is essential to differentiate tissue stiffness.

In addition to the capability of ultrasound elastography to provide quantitative elastic information, the interaction with the examiner to feel the elasticity can help to locate ELASTOGRAM anomalies while moving a ROI in the ultrasound image. Two assistant modes were proposed in this chapter to aid the examiner to perform ultrasound elastography of the tissue and to simultaneously feel the elasticity of the tissue with a haptic device. This system with the teleoperation of a 2D ultrasound probe offers the capabilities to remotely perform ultrasound elastography on a patient or simply to confirm the tissue elasticity displayed in the elastogram. We have demonstrated experimentally a good performance for both teleoperation and impedance haptic control modes. The force feedback using the elastogram was statistically evaluated to determine its reproducibility. However, the estimation of the force feedback from the elastogram assumes a specific stiffness value depending on the Young's modulus of the soft tissue and this may variate between different kinds of soft tissues. Despite of the coarse assumption value of the Young's modulus, the force feedback feeling obtained from the experimental results is promising, and offers the possibility to perform a study with expert physicians to validate this force feedback functionality assistance.

CHAPTER 5 CONCLUSION AND PERSPECTIVES

Conclusion

The research work presented in this thesis focused on the design of a robotic system to assist the examiner in the ultrasound elastography process. The assistance of the robotic system could facilitate the diagnosis of diseases done by the physician. In addition, the system introduced in this work has demonstrated the ability to perform fatiguing tasks that the examiner usually performs. In Chapter 1, we introduced the principles of ultrasound imaging, the state-of-the-art of the ultrasound elastography and the definition of visual servoing. In Chapter 2, we presented a robotic-assisted system for quantitative ultrasound elastography combining three hierarchical robotic tasks. As far as we know, the system presented in Chapter 2 is the first robotic-assisted system that uses direct parameters from the elastography in a robotic task. This contribution showed its capacity to be used with 2D and 3D ultrasound probes. A different methodology to maintain the stiff tissue of interest in the field of view of the ultrasound probe was designed in Chapter 3. This approach took into account motion perturbations in the tissue (e.g., physiological motions), which were not addressed in the system presented in Chapter 2. The approach of Chapter 3 considered a 2D ultrasound probe using a hybrid ultrasound visual servoing process to compensate this perturbation motion with a 6-DOF robot. The automatic motion compensation with the robot led to another contribution, the estimation of the elasticity map of a moving tissue. In Chapter 4, another robotic application with elastography was presented. This last application was the development of a haptic system based on the elasticity map given by the robotic palpation system. In addition, the haptic device was used to teleoperate the ultrasound probe in order to explore the tissue. In the following, we present the general conclusions of the thesis. Afterwards, several ideas The first Chapter presented the three concepts most used in this work. First off, the main concepts of ultrasound imaging modality were introduced. The definition of the image reconstruction from ultrasound was described. Then, the state-of-the-art of the elastography described the most common approaches to obtain the strain map of the tissue. The entire classic approach of elastography was detailed to introduce the idea of the process in which a robot can assist. Afterwards, the few robotic systems related with ultrasound elastography were described. However, these previous robotic systems do not exploit the elastic information in their robotic controllers. This chapter also recalled the basic principle of visual servoing.

Chapter 2 presented one of the major contributions of this thesis, a robotic-assisted system for quantitative ultrasound elastography. The proposed system allows an automatic and real-time generation of the elasticity map of the observed tissues. It is based on the development and implementation of three robotic tasks: palpation motion, automatic centering of stiff tissue, which is the target of interest and orientation of the ultrasound probe. The palpation motion was performed by a force controller that uses the force measured with a force/torque sensor placed between the end-effector of the robot and the ultrasound probe. This controller was designed to continuously apply a periodical pressure on the tissues by controlling the probe velocity component along its axial direction. An algorithm was then proposed based on BMA (Block Matching Algorithm) and optical flow to estimate the elasticity image (elastogram) from pre-and post-compression RF arrays acquired during this palpation task. The second task was employed to automatically center the ultrasound probe with the stiffest tissue in a ROI.

The location of the stiffest tissue was extracted from the elastogram image generated by the palpation task. Then, the location of the stiffest tissue was used as input of a visual servoing process to laterally center it in the FOV of the ultrasound probe. The third robotic task that we proposed was an automatic re-orientation of the probe that can be activated on demand by the examiner to observe the stiff tissue of interest under different view angles. The three tasks were combined by formulating a hierarchical control.

The palpation motion has the highest priority, followed by the automatic centering of the stiffest tissue and lastly, the orientation of the probe has the lower priority. Experimental results obtained on different phantoms demonstrated the feasibility of our proposed assistant robotic system for 2D and 3D elastography. The experiments performed with the 2D probe demonstrated convergence of the robotic system when a ROI was selected in the ultrasound image's FOV. As the use of a 2D probe limited the motion control of the probe in the observation plane, we also proposed an adaptation of our system for the use of a 3D probe. This adaptation gave to the system the possibility of performing 3D elastography and centering the probe on a stiff tissue target that does not initially CHAPTER 5. CONCLUSION AND PERSPECTIVES intersect with the central US observation plane. However, the acquisition time of an ultrasound volume with a motorized probe is important (around 1s per volume) and this limitation caused the system to be slow compared with the one developed with a 2D probe. Additionally, a study of the use of this robotic system to improve the quality of the 2D elastography was also included in this chapter. This study consisted in capturing elastograms at different probe orientations and then averaging them to compute the resulting elastogram, which has better quality than an individual elastogram at one orientation. The quality assessment in the elastogram was performed by measuring the CNRe (contrast-noise ratio elastogram).

The limitations of the robotic system presented in Chapter 2 led us to design an alternative approach to assist in the elastogram estimation. This approach has to be robust to any motion perturbation on the tissue. This system was presented in Chapter 3, where a 2D ultrasound probe was used to perform the compensation of perturbing motions. The robotic system comprised two tasks, palpation motion and motion compensation. The first task, palpation motion, remained the same as in the previous approach. However, the motion compensation using a 2D ultrasound probe was the key issue in the design of the new system. As explained before in the previous approach, a 2D ultrasound probe was limited to in-plane motions. However, prior knowledge of close parallel planes to the initial 2D ultrasound image provides enough information to estimate the motion out-ofplane through 3D image gradients, as presented in [START_REF] Nadeau | Intensity-based ultrasound visual servoing: Modeling and validation with 2-D and 3-D probes[END_REF]. This approach can compensate motion based on a dense image tracking of a ROI. However, the existing approach only considered a rigid model transformation to track the ROI, which in our case is not always sufficient due to the deformation of the tissue. Therefore, we proposed to adapt this dense image tracking method by including a deformable model to track the ROI in the image in order to compensate for the deformations of the ROI due to the compression applied by the palpation motion with the ultrasound probe. This improves the efficiency of the robotic motion compensation task. In addition, the information of the deformable tracking of the ROI was employed to compute the displacement map and then used for the estimation of the elastogram in the ROI. The experimental results demonstrate that the system is capable of compensating a wide range of external motions applied to the tissue, keeping the ROI in the FOV of the ultrasound probe. This improved the system presented in Chapter 2, providing a better alternative for a robotic-assisted system for elastography.

The possibility of measuring the elasticity of the tissue with the assistance of a robotic system helps the examiners to detect stiff tissues. However, the cooperation between the robotic system and the examiner is indispensable to allow a better diagnosis of diseases based on the elasticity of the tissue. Feeling the elasticity of the tissue is a fundamental task in palpation examination, which complements the elastography process.

PERSPECTIVES

Therefore, in Chapter 4, we presented a haptic system to apply force feedback based on the motion of a virtual probe and the estimated elastogram of the tissue. This approach was designed as a master-slave system, where the master was a haptic device and the robot performing palpation motion with the ultrasound probe was the slave. This system included a teleoperation of the probe with a haptic device to explore the tissue while the ultrasound probe was always in contact with the tissue using the palpation motion task.

The approach was experimentally tested over the abdominal phantom used in the other two robotic systems previously described, and the results of the haptic force feedback were analyzed. The results showed that the haptic system is an excellent approach that can help doctors to just not only measure the elasticity, but to be able to feel it through a haptic device.

Perspectives

There are some important perspectives of this work which are classified in two categories accordingly to the time of development: short-term and long-term perspectives.

Short-term perspectives

Many studies can be performed with the robotic framework presented in this work. For example, the use of a 3D ultrasound probe with faster acquisition of ultrasound data can be explored to expand the approach presented in Chapter 2. The exploration and comparison of different ultrasound elastography approaches with the robotic system can also provide another perspective in determining which technique performs better with the robotic system. A multi-modality image registration process can also improve the system presented in Chapter 2. The information of the 3D stiffness map of the tissue can be obtained beforehand with MRE (Magnetic resonance elastography) to apply a multi-modality image registration with the 2D stiffness map estimated in real-time with the ultrasound. This approach can be applied in laparoscopy, guiding the surgeon in the location of stiff tissues. The multi-modality image registration for elastography can also be used to avoid the estimation of the tissue stiffness and instead using the corresponding stiffness value computed by the MRE to calculate the force feedback in the haptic system presented in Chapter 4.

The occlusions on the ultrasound propagation due to bones or other artifacts are cause of low quality in the image. The quality in the acquisition of the ultrasound information is fundamental to compute the elasticity of the tissue. Recently, a new robotic 128 CHAPTER 5. CONCLUSION AND PERSPECTIVES approach to improve the ultrasound imaging quality was presented in [START_REF] Chatelain | Confidence-driven control of an ultrasound probe[END_REF]. This work is based on the optimization of the confidence map, which is an image representation of the ultrasound quality. Therefore, the adaptation of this approach to the robotic systems presented in this thesis could be an advantageous perspective, leading to increase the quality of the elastogram of the tissue.

In terms of applications, guidance on the insertion of a needle for a biopsy using elastography can be an interesting approach. This application can be obtained with the combination of the robotic needle insertion approach presented in [START_REF] Chevrie | Online prediction of needle shape deformation in moving soft tissues from visual feedback[END_REF] and the approach presented in Chapter 2 to locate the stiffest tissue in a region of interest inside the field of view of an ultrasound probe. The real-time computation of the location of the stiffest tissue, usually corresponding to a tumor, would be the goal to achieve a robotic needle insertion system.

Long-term perspectives

The estimation of tissue deformation can be improved by an online registration between the real-time ultrasound information and prior knowledge of the geometrical structure of the tissue. This kind of approach has been applied in radiotherapy of the neck [START_REF] Coevoet | Registration by interactive inverse simulation: application for adaptive radiotherapy[END_REF] where weight loss during several weeks of therapy modifies the volume of anatomical structures. This approach uses a finite element method (FEM) obtained from the first CT scans. Afterwards, an interactive registration of few points is performed to obtain the deformation and biomechanic parameters of the soft tissue helping in the limitation of radiations.

Currently, soft robots have been studied for minimal invasive surgery [START_REF] Cianchetti | Soft Robots in Surgery[END_REF]. This kind of robots provides higher dexterity and flexibility than the classic robots. The mobility is one of the most interesting elements in a soft robot, it can help in performing elastography of regions occluded by bones in real-time. The design of a soft robot to perform elastography is a perspective that can help surgeons to localize and remove tumors in a minimal invasive surgery.

Finally, the three robotic-assisted systems presented in this dissertation are just the beginning in the research of a robotic tool that someday can be employed in hospitals to facilitate the diagnosis of diseases or the planning of surgery. The approaches for a robotic system to assist in elastography presented in this work were only evaluated on phantoms. However, the main goal is to take these approaches to real medical applications. This is a long process that involves evaluations of the systems with in-vivo tissues, which requires the approvement of hospitals and clinicians. Despite the accuracy of an 129 5.2. PERSPECTIVES industrial robot, the volume and safety of this kind of robots are limitations that should be considered for a certified medical robotic system [START_REF] Taylor | Medical robotics in computer-integrated surgery[END_REF]. Therefore, the implementation of the proposed robotic assistance tasks on a smaller and safer robot is a necessary step.

Figure 1 :

 1 Figure 1: Cadre robotique général pour l'élastographie ultrasonore quantitative.

Le chapitre 5 3 . 1 4 . 1 5 Conclusion and perspectives 125 5 . 1

 5314151 présente la conclusion générale de la thèse et propose des perspectives à court et à long termes de ces travaux. ix x Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Dense visual tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Similarity measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Warp transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Rigid transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.5 Non-rigid transformations . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.6 Strain estimation based on optical flow . . . . . . . . . . . . . . . . . . 3.3 Motion compensation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 6-DOF motion compensation by dense visual servoing . . . . . . . . . 3.3.2 Control fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii CONTENTS 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Teleoperation and haptic feedback based on elastogram 103 Haptic feedback . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Haptic force feedback from elastogram . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Force estimation from elastogram . . . . . . . . . . . . . . . . . . . . . 4.2.2 Impedance force feedback system . . . . . . . . . . . . . . . . . . . . . 4.2.3 Robotic teleoperation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Short-term perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Long-term perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 :

 2 Figure 2: General robotic framework for quantitative ultrasound elastography.

Figure 1 . 1 :

 11 Figure 1.1: Sound propagation through two different media. Left illustration depicts the change in the wavelength of the sound for every media. Right sketch shows the orientation of the sound when is induced (θ i ), transmitted (θ t ) and reflected (θ r ) between the two media.

Figure 1 .Figure 1 . 2 :

 112 Figure 1.2 shows the representation of the reflections of the sound wave of incidence produced by one scatterer point.

Figure 1 . 3 :

 13 Figure 1.3: RF signal recorded from a real transducer.

Figure 1 . 4 :

 14 Figure 1.4: RF envelope detection. Red curve represents the envelope of the RF signal in green.

Figure 1 . 5 :

 15 Figure 1.5: B-mode image reconstructed from the RF signals envelopes. (a) b-mode image reconstructed only from the envelope detection for a set of scan lines. (b) b-mode image obtained with the logarithmic compression of the envelope detection for the same set of scanlines used in (a).

Figure 1 . 6 :

 16 Figure 1.6: Geometry of a linear ultrasound probe.

  7. The elements of the transducer are positioned along an arc of a circumference with center at F p . The metric coordinates of every point (x, y) inside the field of view (FOV) of the probe are related to the polar coordinates (r, θ) by x = r sin(θ), (1.14)

Figure 1 . 7 :

 17 Figure 1.7: Geometry of a convex ultrasound probe.

Figure 1 . 8 :

 18 Figure 1.8: Pre-scan and post-scan b-mode images. (a) pre-scan and (b) resulted postscan b-mode images obtained with a convex ultrasound probe.

1. 2 .Figure 1 . 10 :

 2110 Figure 1.10: (a) RF signal analogy with a succession of springs S 1 , S 2 and S 3 . (b) the Hooke's law scheme.

Figure 1 . 11 :

 111 Figure 1.11: Strain profile for three consecutive springs. Springs 1 and 3 have the same length l and stiffness value k in the pre-compression state.

Figure 1 . 12 :

 112 Figure 1.12: Results where the MRE demonstrates increasing liver stiffness values with increasing stage of fibrosis (Figure taken from [89]). The top row shows share wave images from four patients with biopsy-proven hepatic fibrosis from stage 1 to 4. The lower row shows corresponding elastograms for these patients.

Figure 1 . 13 : 1 . 2 . 5 . 1

 1131251 Figure 1.13: Classification of ultrasound elastography approaches. The principle of each category is described in the section indicated at the bottom of the box.

Figure 1 . 14 :

 114 Figure 1.14: Cross-correlation between the segments of RF signals s pr e and s post . The absolute cross-correlation is normalized to one and the lag where the amplitude is maximum is indicated.

Figure 1 . 15 :Figure 1 . 16 :

 115116 Figure 1.15:Transient elastography process with an ultrasound transducer. First, the shear wave is produced by the vibrator and then, the ultrasound is emmited. This process produces the elastic information of the tissue scanned by the ultrasound transducer.

Figure 1 . 17 :

 117 Figure 1.17: Liver stiffness measurement using FibroScan (Figure taken from [25]).

Figure 1 . 18 :

 118 Figure 1.18: Results of the snake robot performing palpation of a prostate phantom presented in [76]. (Left) shows the experiment setup and (right) presents the stiffness of the artificial prostate.

1. 3 .Figure 1 . 19 :

 3119 Figure 1.19: View of the da Vinci console displaying the real-time images including the elastogram (Figure taken from [10]).

Figure 1 .

 1 Figure1.20: Force control scheme of the robotic system presented in[START_REF] Bell | Force-controlled ultrasound robot for consistent tissue pre-loading: Implications for acoustic radiation force elasticity imaging[END_REF] for applying the ARFI required for the tissue elasticity measurement.

Figure 1 . 21 :

 121 Figure 1.21: Visual servoing closed-loop.

Figure 1 . 22 :

 122 Figure 1.22: The two configurations of the visual sensor location in visual servoing.

Figure 1 . 23 :

 123 Figure 1.23: Robot joints q for a 6-DOF robot, and the pose r of the end-effector.

Figure 2 . 1 :

 21 Figure 2.1: Scheme and workflow of the experimental setup.

CHAPTER 2 .Figure 2 . 2 :

 222 Figure 2.2: Viper s850 robot.

2. 1

 1 . EXPERIMENTAL SETUP (a) SonixTouch ultrasound system. (b) Ultrasound transducer 4DC7-3/40.

Figure 2 . 3 :

 23 Figure 2.3: Ultrasound equipment used for the experiments.

  (a) Abdominal phantom ABDFAN US-1B. (b) Two-layers gelatin phantom containing two duck gizzards.

Figure 2 . 4 :

 24 Figure 2.4: Phantoms used in the experiments.

Figure 2 . 5 :

 25 Figure 2.5: Cartesian reference frames attached to the robotic arm.

b

  as the rotation matrix and translation vector, respectively. t a b × is the three-by-three skew-symmetric matrix representation of t a b .

Figure 2 . 6 :

 26 Figure 2.6: Desired sinusoidal force applied by the force controller.

Figure 2 . 7 :

 27 Figure 2.7:Measured force when the force control law is applied. Green curve is the desired force variation. The red and blue curves correspond to the force measured when applying Equations (2.5) and (2.6) respectively.

Figure 2 . 8 .

 28 Two RF frames are grabbed before and after the tissue compression and represent respectively the pre-and post-compression states. A region of interest (ROI) or a volume of interest (VOI), in case of the use of a 2D or 3D ultrasound probe, is selected from the RF data to estimate the elastogram and display it on the b-mode image.

Figure 2 . 8 :

 28 Figure 2.8: Elastography process. From left to right: mechanical compression with the ultrasound probe over a static soft tissue; RF frames acquired for the pre-and postcompression states; b-mode image with an elastogram overlaid on a ROI.

Figure 2 . 9 :

 29 Figure 2.9: Displacement estimation. Parameters in the motion estimation process using the RF data observed in the ROIs of the pre-and post-compressed frames.

  [START_REF] Cavusoglu | Multirate simulation for high fidelity haptic interaction with deformable objects in virtual environments[END_REF] where (i c , j c ) is the center point of the search region as shown in Figure2.10. At this stage, we have shown the estimation of the motion for one block. The same process is repeated for the next blocks, but we need to define their position. We define the position changing as a shifting process, where the next block position can represent a block overlapping with respect to one or more previous blocks (seeFigure 2.10). The percentage of blocks overlapping is related to the resolution of the axial and lateral displacements

Figure 2 . 10 :

 210 Figure 2.10: Parameters involved in the motion estimation algorithm using the information of the post-compression RF frame (element access through g (i , j )). The grid inside the search region represents the integer step between the RF elements.

Figure 2 . 11 :

 211 Figure 2.11: Displacement maps from motion estimation. After obtaining the displacements for every block, their axial and lateral components are stored in two arrays V 0 (i , j ) and U 0 (i , j ), respectively.

Figure 2 . 12 :

 212 Figure 2.12: LSQ strain estimation. (a) strain estimation parameters for the sample at (i s , j s ) in the displacement map, and the respective (b) strain map.

Figure 2 . 13 :

 213 Figure 2.13: Generation of the ground truth using a FEM.

Figure 2 .

 2 Figure 2.14 shows few scatterers as red points in a frame acquisition by the simulated ultrasound machine (similar to a real ultrasound machine). The same figure also shows how the scatterers have been displaced due the mechanical compression. Note that in order to work, FIELD II requires a number of scatterers of order 10 3 points in a random Gaussian distribution.

Figure 2 . 14 :

 214 Figure 2.14: Scatterers position in the process of compression.

Figure 2 . 15 :

 215 Figure 2.15: Three input models used to simulate the compression by FEM.

Figure 2 . 16 :

 216 Figure 2.16: Output due the compression by FEM.

Figure 2 . 17 :Figure 2 . 18 :

 217218 Figure 2.17: Output FIELD II using B-Mode.

Figure 2 . 19 :

 219 Figure 2.19: RF volume for one motor sweep.

Figure 2 . 20 :

 220 Figure 2.20: 3D elastogram reconstruction based on the 2D process.

Figure 2 .

 2 Figure 2.20 shows the process to compute one elastogram in 3D based on the 2D

Figure 2 . 21 :

 221 Figure 2.21: Connected components process to label different regions for 2D and 3D images.

Figure 2 . 25 :Figure 2 .

 2252 Figure 2.25: 3D scan conversion. We show the parameters needed to convert any point s in RF units to a point p in metric units.

Figure 2 . 26 :

 226 Figure 2.26: Probe orientation.

Figure 2 . 27 :

 227 Figure 2.27: 3D probe orientation.

Figure 2 . 29 :

 229 Figure 2.29: Block diagram of the 2D case implemented in the multi-thread software application.

  +5 • and θ 4 = θ 3 +5 • . The curves of the error evolution for the three tasks are shown in Figure 2.30.We can see that the force error ranges between ±1N due to the sinusoidal desired force variation. Once the ROI is selected (at time 20 s), the centroid of the elastogram is computed, and it is sent to the automatic centering control task. The object centering error decreases towards zero but still exhibits a low remaining oscillation of ±3 mm due to the elastogram noise. However, the ROI is horizontally maintained close to the image center by the visual servoing task even when the user successively changed the probe desired angles at times 21, 90, 126, 167 and 205 s, keeping automatically the object of interest in the field of view.

Figure 2 .Figure 2 . 30 :

 2230 Figure 2.30: Phantom experiment. Evolution of the system during the experiment. (a) Force error evolution, with F 0 = 5N and ∆ F = 2N. (b) Probe orientation error evolution. (c) Horizontal target centering error evolution. (d) Velocities applied to the 3-DOF involved in the control law.

Figure 2 . 31 :

 231 Figure 2.31: Experiment with different probe orientations: (a)-(e) the probe oriented at different angles, from left to right the angles are θ 0 = -10 • , θ 1 = -5 • , θ 2 = 0 • , θ 3 = 5 • and θ 4 = 10 • . (f)-(j) b-mode image where the target is centered with the image and the elastogram ROI overlaid for each probe orientation. (k)-(o) elastograms obtained for each probe orientation shown in the ROI of the images.

42 Table 2 . 4 :

 4224 71.45 65.87 87.52 87.02 127.Comparison of the CNR e (in dB) of the estimated elastography images at the different probe orientations of the experiment and their mean.

Figure 2 . 32 :

 232 Figure 2.32: Mean of aligned elastograms obtained from 5 probe orientations.

Figure 2 . 33 :

 233 Figure 2.33: Display of the three orthogonal planes and the VOI with the 3D elastogram.

Figure 2 . 34 :

 234 Figure 2.34: Short diagram of the implemented multi-thread software application.

Figure 2 . 35 :

 235 Figure 2.35: Experiment with a gelatin phantom containing two duck gizzards.

2. 7 .

 7 CONCLUSIONservoing task start to variate. The convergence is reached at t ≃ 150s when the probe has been automatically aligned with the stiff tissue of interest (duck gizzard) by following an exponential decrease of the visual error e t as expected. A change of the orientation of the probe is introduced by the user at time t ≃ 160s and we can see at this point the variation of ω x , ω y and ω z that allow reaching the desired orientation (s * θ = [ π 36 -π 36 π 18 ] rad). At this time all the tasks are activated which lead us to keep the target in the center of the ultrasound image.

  Angular velocities (e) Probe position at the VOI selection (t =72s) (f) Probe position at t =80s (g) Probe position at t =110s (h) Probe position at t =150s (i) Probe orientation at t =172s

Figure 2 . 36 :

 236 Figure 2.36: Evolution through time of the 3D reference experiment. (a)-(d) evolution of the errors and velocities during the experiment. (e)-(h) evolution of the pose of the ultrasound probe during the centering of the object in the FOV with the 3D ultrasound images at the top. The stiffest object (green color) is also overlaid in the 3D ultrasound images. (i) probe orientation control of the ultrasound probe.

Figure 3 . 1 :

 31 Figure 3.1: Surface tracking of a beating heart [61]. Left and right images of the stereo camera and the TPS surface approximation.

  The main advantage of dense image registration techniques is the suppression of an image feature extraction step, therefore avoiding errors due to false detection or bad segmentation in the image. Let us define I t ∈ R m×n and I c ∈ R m×n as the initial and current image templates, respectively. The templates are composed of N = m × n pixels.

Figure 3 . 2 :

 32 Figure 3.2: Ultrasound images acquired after applying small motion to the probe.(a) reference image. (b)-(c) images acquired after small motion of the probe. (c) an oclussion (shadow) at the right side of the image is present due to the lack of echogenic gel between the probe and the tissue surface.

  p) is a function that maps a point or set of points x to a new location x ′ by applying a geometric transformation, as illustrated in Figure 3.3. The parameters p in the warp function depend on the type of geometric transformation to use. These transformations are classified in two groups, rigid and non-rigid transformations. The most commonly used transformations in visual tracking are briefly explained in the next section.

Figure 3 . 3 :

 33 Figure 3.3: Arbitrary warp transformation applied to a set of points x.

Figure 3 . 4 :

 34 Figure 3.4: Example of performance in dense rigid tracking. (a) initial image (i = 0) of an image sequence (green rectangle is the region to be tracked). (b) last image (i = 500) in the sequence. (c) template image from (a) at top and warped image from (b) at bottom. (d) error image computed from Îc (w (x, p)) -I t (x).

Figure 3 . 5 :

 35 Figure 3.5: Example of the performance in dense affine tracking. (a) initial image (i = 0) of an image sequence (green rectangle is the region to be tracked). (b) last image (i = 500) in the sequence. (c) template image from (a) at top and warped image from (b) at bottom. (d) error image computed from Îc (w (x, p)) -I t (x).

Figure 3 . 6 :

 36 Figure 3.6: Example of the performance in dense FFD tracking. (a) initial image (i = 0) of an image sequence (green rectangle is the region to be tracked). (b) last image (i = 500) in the sequence. (c) template image from (a) at top and warped image from (b) at bottom. (d) error image computed from Îc (w (x, p)) -I t (x).

Figure 3 . 7 :

 37 Figure 3.7: TPS image transformation. Image deformation using the TPS warp function. Red dots represent the location of the control points.

Figure 3 . 8 :

 38 Figure 3.8: Example of the performance in dense TPS tracking. (a) initial image (i = 0) of an image sequence (green rectangle is the region to be tracked). (b) last image (i = 500) in the sequence. (c) template image from (a) at top and warped image from (b) at bottom. (d) error image computed from Îc (w (x, p)) -I t (x).

Figure 3 . 9 :

 39 Figure 3.9: Displacement map obtained from the deformation of the ROI. x and x ′ represent respectively the coordinates of one pixel in I t and its new coordinates in the current image I c after tissue deformation. D(x) is the vector displacement from x to x ′ .

Figure 3 .

 3 Figure 3.10 shows the workflow of our proposed method. The main problem produced by a moving object is the generation of noise in the elastogram, due to the introduction of non-axial motion. To counter this, our system combines ultrasound dense visual servoing, force control and the non-rigid motion estimation to compute the elastogram of moving tissue presented in last section. In the next subsections we describe the elements of the proposed control system featured in Figure3.10.

Figure 3 . 10 :

 310 Figure 3.10: Proposed methodology to obtain the strain map of a moving tissue. This diagram shows the steps to estimate the elastogram from two ROIs in the ultrasound images I t and I c . These images are the reference and current images respectively. I w represents the image I c modified with a non-rigid transformation (w ) that is performed from the motion estimation in order to reduce the absolute difference with I t . [U, V] are the lateral and axial displacement maps computed between I c and I t after the motion estimation.

Figure 3 . 11 :

 311 Figure 3.11: Possible motions with a 2D US probe.

Figure 3 . 13 :

 313 Figure 3.13: Evolution through time of the reference experiment. (a) to (d) show the evolutions of the velocities, measured force and visual error during the experiment.

Figure 3 .

 3 Figure 3.13c.

Figure 3 . 14 :

 314 Figure 3.14: Perturbations introduced to the phantom. (a) to (d) show some of the states during the experiment. (e) to (h) show the b-mode images observed at these states. (i) to (l) show the strain maps for every state.

Figure 3 .Figure 3 .Figure 3 . 15 :

 33315 Figure 3.15 shows three elastograms obtained under different conditions. InFigure3.15a, the strain map has been estimated using only oscillatory force control without the automatic motion compensation when the phantom was motionless. We can see

Figure 4 . 1 :

 41 Figure 4.1: Types of haptic force feedback.

Figure 4 . 2 :

 42 Figure 4.2: Configurations of kinesthetic haptic device.

Figure 4 . 3 :

 43 Figure 4.3: Virtuose 6D (Haption S.A.). Haptic device used in this thesis work.

Figure 4 . 7 .Figure 4 . 6 :

 4746 Figure 4.6: Short block diagram of the system proposed in this chapter.

Figure 4 . 7 :

 47 Figure 4.7: Force estimation based on strain information.

4. 2 .r = σ 2 x + σ 2 y( 4 . 7 )Figure 4 . 8 :

 2224748 Figure 4.8: Thickness of the compressed tissue with an ultrasound probe.

Figure 4 . 9 illustrates

 49 the principle that consists in moving the ROI to follow the displacement of the user measured by the handler of the haptic device. This figure also shows the Cartesian frames F b and F h corresponding to the base and the handler of the haptic device, respectively, and F I corresponding to the ultrasound image. The point (u r , v r ) is the origin point of the ROI with respect to the ultrasound image frame F I , and the point (u c , v c ) is the center of the ROI. If the user applies motion at the handler of the haptic device F h , then the point (u c , v c ) is shifted with a displacement ∆ d proportional to the displacement of the handler.

Figure 4 . 9 :

 49 Figure 4.9: Moving elastogram inside of the ultrasound image.

Figure 4 . 10 :

 410 Figure 4.10: Handler displacement into the ultrasound image.

. 10 )

 10 The displacement ∆ d of the point (u c , y c ) in the ultrasound image affects directly the position of the ROI where the elastogram is estimated. This produces a displacement of ∆ d in all the elements inside the ROI, and the motion of the point (u c , y c ) must be bounded in the region R I ∈ ( w r 2 (w I -w r2 )] × ( h r 2 (h I -h r 2 )] (seeFigure 4.11). Setting boundaries for the motion of the point (u c , y c ) is necessary to ensure the estimation of the elastogram in a ROI of size w r × h r .

Figure 4 . 11 :

 411 Figure 4.11: Motion boundary R I of the point (u c , y c ) in the ultrasound image.

Figure 4 .

 4 Figure 4.13). Therefore, the new relative pose expressed by the homogeneous matrix M ∆

Figure 4 . 12 :

 412 Figure 4.12: Main components of the master-slave teleoperation system and the connexions between them.

Figure 4 . 13 :

 413 Figure 4.13: Cartesian frames used for teleoperation of the probe.Once the relative pose is expressed in the frame F cp , we can extract the six displacement components, three angular displacements (∆θ x , ∆θ y , ∆θ z ) and three linear displacements (∆t x , ∆t y , ∆t z ). The linear displacements are computed from the translational vector inside M ∆ . On the other hand, the angular displacements are obtained through the rotation matrix R ∆ contained in M ∆ and expressed as three angular components (∆θ x , ∆θ y , ∆θ z ) with Rodrigues' formula. These six displacement components are the desired values that the ultrasound probe should reach with respect to the initial probe position relative to the robot's base frame F r , r M cp0 " and we can enclose them in a desired feature vector as,

Figure 4 . 14 :

 414 Figure 4.14: Experimental setup of the proposed haptic system.

Figure 4 . 15 :

 415 Figure 4.15: Proposed multithread workflow for the implementation. Directional connections between the threads are depicted with red lines.

Figure 4 .

 4 Figure 4.16d. Then, the palpation motion task is activated, initiating the force control to reach contact with the phantom as shown in Figure 4.16e.Figure 4.17a shows at t =∼ 2.5s

Figure 4 .

 4 17a shows at t =∼ 2.5s the beginning of the force variation needed for the palpation motion. The teleoperation of the US probe with the haptic device starts at t =∼ 13s as indicated with the black arrow in the plot illustrating the evolution of the teleoperation errors (Figure4.17d). At the same time, we can also observe inFigures 4.17b

Figure 4 . 16 :

 416 Figure 4.16: Teleoperation system states. (a) initial pose of the handler and (d) the corresponding pose of the ultrasound probe with the (h) resulting ultrasound image before the contact with the phantom. (b) rotation around z-axis of the handler and (f) the corresponding pose of the ultrasound probe with the (j) resulting ultrasound image. (c) arbitrary pose of the handler and (g) the corresponding pose of the ultrasound probe with the (k) resulting ultrasound image.

Figure 4 . 17 :

 417 Figure 4.17: Measured force, velocities and errors in the teleoperation system.

Figure 4 .

 4 Figure 4.18 the plot that represents the haptic force feedback applied to the handler of

Figure 4 . 18 :

 418 Figure 4.18: Result of the force feedback of the impedance system. First row shows the different states while moving the handler of the haptic device. The motion of the ROI containig the elastogram is shown in the second row for the different states of the handler motion. The temporal evolution of the force feedback is ploted at the bottom and the position of the states are indicated with red arrows.

Figure 4 .

 4 [START_REF] Chevrie | Online prediction of needle shape deformation in moving soft tissues from visual feedback[END_REF] shows the force feedback average after the 50 repetitions of the green path, where the black line is the force feedback average and green area represents the interquartile range (IQR). The small IQR for all the positions shown in the plot of Figure4.20 describes the standard deviation (SD). Based on the observed small variation (maximum SD of 0.21N) of the force feedback after 50 repetitions of the green path illustrated in Figure4.19, we can conclude that our force feedback measurement in the ROI is highly reproducible.

Figure 4 . 19 :Figure 4 . 20 :

 419420 Figure 4.19: Repetitive motion path of the virtual probe. The four images correspond to the four corners of the square path. The green path was used to measure the standard deviation of the force feedback.

5. 1 .

 1 CONCLUSIONfor perspectives of this work are provided.

  .[START_REF] Chaumette | Visual servo control, part II: Advanced approaches[END_REF] Results of the robotic system presented in[START_REF] Sen | Enabling technologies for natural orifice transluminal endoscopic surgery (N.O.T.E.S) using robotically guided elasticity imaging[END_REF] . . . . . . . . . . . . . . . . . . 1.19 Results of the robotic system presented in [10] . . . . . . . . . . . . . . . . . . .13 Generation of the ground truth using a FEM. . . . . . . . . . . . . . . . . . . . 2.14 Scatterers position in the process of compression. . . . . . . . . . . . . . . . 2.15 Three input models used to simulate the compression by FEM. . . . . . . . . 2.16 Output due the compression by FEM. . . . . . . . . . . . . . . . . . . . . . . . Experiment with different probe orientations in 2D . . . . . . . . . . . . . . . 2.32 Mean of aligned elastograms obtained from 5 probe orientations. . . . . . . 2.33 Display of the three orthogonal planes and the VOI with the 3D elastogram. 2.34 Short diagram of the implemented multi-thread software application. . . . 2.35 Experiment with a gelatin phantom containing two duck gizzards. . . . . . 2.36 Evolution through time of the reference experiment (3D case) . . . . . . . . Strain maps obtained under different conditions. . . . . . . . . . . . . . . . . 4.1 Types of haptic force feedback . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Configurations of kinesthetic haptic device . . . . . . . . . . . . . . . . . . . . 4.3 Virtuose 6D (Haption S.A.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Workflow of the impedance feedback . . . . . . . . . . . . . . . . . . . . . . . 4.5 Workflow of the admittance feedback . . . . . . . . . . . . . . . . . . . . . . .

	LIST OF FIGURES	LIST OF FIGURES LIST OF FIGURES
	xvii LIST OF TABLES 4.20 Force feedback average after 50 repetitions of a path . . . . . . . . . . . . . . 122 2.31 xviii xx 3.15 xix 1.1 Comparison between some of the ultrasound elastography approaches. . .

1 Cadre robotique général pour l'élastographie ultrasonore quantitative. . . . v 2 General robotic framework for quantitative ultrasound elastography. . . . . 3 1.1 Sound propagation through two different media . . . . . . . . . . . . . . . . 9 1.2 Scatterer point reflections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3 RF signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.4 RF envelope detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.5 B-mode image reconstructed from the RF signals envelopes . . . . . . . . . 13 1.6 Geometry of a linear ultrasound probe . . . . . . . . . . . . . . . . . . . . . . 14 1.7 Geometry of a convex ultrasound probe . . . . . . . . . . . . . . . . . . . . . . 15 1.8 Pre-scan and post-scan b-mode images . . . . . . . . . . . . . . . . . . . . . . 16 1.9 Organization of the state-of-the-art of elastography . . . . . . . . . . . . . . . 17 1.10 RF signal analogy with a succession of springs and Hooke's law scheme . . . 18 1.11 Strain profile for three consecutive springs . . . . . . . . . . . . . . . . . . . . 18 1.12 Magnetic resonance elastography of liver fibrosis [89] . . . . . . . . . . . . . 20 1.13 Classification of ultrasound elastography approaches . . . . . . . . . . . . . 22 xv LIST OF FIGURES 1.14 Cross-correlation between two segments of RF signals . . . . . . . . . . . . . 1.15 Transient elastography process with an ultrasound transducer . . . . . . . . 1.16 Experimental setup for transient ultasound elastography [12] . . . . . . . . . 1.17 Liver stiffness measurement using FibroScan (Figure taken from [25]). . . . 11.20 Force control scheme of the robotic system presented in [7]. . . . . . . . . . 1.21 Visual servoing closed-loop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.22 The two configurations of the visual sensor location in visual servoing. . . . 1.23 Robot joints q for a 6-DOF robot, and the pose r of the end-effector. . . . . . 2.1 Scheme and workflow of the experimental setup. . . . . . . . . . . . . . . . . 2.2 Viper s850 robot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Ultrasound equipment used for the experiments. . . . . . . . . . . . . . . . . 2.4 Phantoms used in the experiments. . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Cartesian reference frames attached to the robotic arm. . . . . . . . . . . . . 2.6 Desired sinusoidal force applied by the force controller. . . . . . . . . . . . . 2.7 Measured force when the force control law is applied . . . . . . . . . . . . . . 2.8 Elastography process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Displacement estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Parameters involved in the motion estimation algorithm . . . . . . . . . . . . xvi LIST OF FIGURES 2.11 Displacement maps from motion estimation . . . . . . . . . . . . . . . . . . . 2.12 LSQ strain estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.17 Output FIELD II using B-Mode. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Results using simulation data from FEM+FIELD II. . . . . . . . . . . . . . . . 2.19 RF volume for one motor sweep. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 3D elastogram reconstruction based on the 2D process. . . . . . . . . . . . . 2.21 Connected components process to label different regions. . . . . . . . . . . . 2.22 Centroid estimation of the biggest region . . . . . . . . . . . . . . . . . . . . . 2.23 Cartesian frames in the ROI for automatic centering process. . . . . . . . . . 2.24 Axes for the frame F p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.25 3D scan conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.26 Probe orientation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.27 3D probe orientation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.28 Plots of the probe control velocities . . . . . . . . . . . . . . . . . . . . . . . . 2.29 Block diagram of the 2D case implemented in the multi-thread software application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.30 Phantom experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Surface tracking of a beating heart [61] . . . . . . . . . . . . . . . . . . . . . . 3.2 Ultrasound images acquired after applying small motion to the probe . . . . 3.3 Warp transformation applied to a set of points x . . . . . . . . . . . . . . . . . 3.4 Example of performance in dense rigid tracking. . . . . . . . . . . . . . . . . 3.5 Example of the performance in dense affine tracking. . . . . . . . . . . . . . 3.6 Example of the performance in dense FFD tracking. . . . . . . . . . . . . . . 3.7 TPS image deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Example of the performance in dense TPS tracking. . . . . . . . . . . . . . . . 3.9 Displacement map obtained from the deformation of the ROI . . . . . . . . 3.10 Methodology to obtain the strain map of a moving tissue. . . . . . . . . . . . 3.11 Possible motions with a 2D US probe. . . . . . . . . . . . . . . . . . . . . . . . 3.12 Spatial derivative filters (case of three parallel slices). . . . . . . . . . . . . . . 3.13 Evolution through time of the reference experiment. . . . . . . . . . . . . . . 3.14 Perturbations introduced to the phantom . . . . . . . . . . . . . . . . . . . . 4.6 Block diagram of the system presented in Chapter 4 . . . . . . . . . . . . . . 4.7 Force estimation based on strain information . . . . . . . . . . . . . . . . . . 4.8 Thickness of the compressed tissue . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Moving elastogram inside of the ultrasound image . . . . . . . . . . . . . . . 4.10 Handler displacement into the ultrasound image . . . . . . . . . . . . . . . . 4.11 Motion boundary of the point (u c , y c ) . . . . . . . . . . . . . . . . . . . . . . . 4.12 Teleoperation system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Cartesian frames used for teleoperation of the probe . . . . . . . . . . . . . . 4.14 Experimental setup of the proposed haptic system . . . . . . . . . . . . . . . 4.15 Proposed multithread workflow for the implementation . . . . . . . . . . . . 4.16 Teleoperation system states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Measured force, velocities and errors in the teleoperation system. . . . . . . 4.18 Result of the force feedback of the impedance system . . . . . . . . . . . . . 4.19 Repetitive motion path of the virtual probe . . . . . . . . . . . . . . . . . . . . 2.1 Force/torque sensor range and resolution. . . . . . . . . . . . . . . . . . . . . 2.2 Parameters used in FIELD II. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Parameters for the control simulation test. . . . . . . . . . . . . . . . . . . . . 2.4 Quality comparison probe orientation and average. . . . . . . . . . . . . . . . 3.1 Nonlinear optimization strategies . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Features of nonlinear optimization strategies. . . . . . . . . . . . . . . . . . . 3.3 Performance evaluation of the transformations in the visual tracking system. xxi LIST OF TABLES xxii

Table 1 .

 1

	-static	Mechanical	Manual compression	Compatibility with all ultrasound probes	Operator dependent
	Remote palpation Transient	Shear wave	ARFI Transient force	Can assess deeper located tissue Accuracy for liver fibrosis staging	Special or additional equipment
	Supersonic		ARFI	Fast acquisition

1: Comparison between some of the ultrasound elastography approaches.

Table 2

 2

	Axis	Sensing range Resolution
	F x , F y	65 N	12.5 mN
	F z	200 N	25 mN
	T x , T y and T z	5 N•m	0.75 µN•m

.1. As shown in Figure

2

.2, the robot is holding an ultrasound probe plugged to its end-effector.

Table 2 .

 2

1: Force/torque sensor range and resolution.

Table 2 . 2 :

 22 2.2. AUTOMATIC PALPATIONscatterers in order to obtain the post-compressed state of the virtual organ. FIELD II was then used to simulate the RF data of the pre-and post-compression states of the virtual tissue. In practice, FIELD II has to be initialized with the parameters of an ultrasound probe and in our case we set it with the parameters of our real probe (4DC7-3/40) in 2D mode, shown in Table2.2. Parameters used in FIELD II.

	Parameter	Value
	Transducer center frequency [Hz]	3.5 × 10 6
	Sampling frequency [Hz]	4.0 × 10 7
	Speed of sound [m/s]	1540
	Wave length [m]	1540/3.5 × 10 6
	Width of element[m]	4.25 × 10 -4
	Kerf [m]	5.5 × 10 -5
	Number of elements in the transducer	128

Table 2 . 3 :

 23 Parameters for the control simulation test.

		Parameter		Initial value	Value at t =2s
	2D 3D 2D and 3D	s t s * θ s t s * θ F 0 ∆F	[ π 6	9 cm π 6 rad [9 3] ⊤ cm π 32 π 9 ] ⊤ rad [ 5π 18 3 N 2 N	-3 rad π -7π 36 -5π 36 ] ⊤ rad --

Table 3 . 2 :

 32 Features of nonlinear optimization strategies.

		Convergence efficiency	Advantages	Disadvantages
	Forward additional	Low	Simple approach	Jacobian computed every at iteration
	Direct com-positional	Medium	Intuitive approach and faster convergence than	Jacobian computed every at iteration
			forward additional ap-
			proach
	Inverse com-positional	High	Jacobian computed be-fore iterations	Requires good knowledge about the warp function w

(a) (b)

1 2 . . . k 1 1 . . . 1 (1.34)whereξ(k) = 12 k(k 2 -1)is the variance of the estimated displacement and k is the number of samples in the RF signal segment. Then, by convolving the vector g(k) with the displacement ∆ t , the strain values are estimated. This filter was designed such that ∆ t contains only two displacements (at n i and n f ) reducing the computational time considerably in the strain estimation. The LSQ strain estimator has been used to obtain a real-time ultrasound elastography[START_REF] Turgay | Identifying the mechanical properties of tissue by ultrasound strain imaging[END_REF].

The special Euclidean group SE (2) = SO(2) ⋉ R

with SO(2) as the special orthogonal group.

ACKNOWLEDGEMENTS I would like to express my sincere gratitude to my supervisor Alexandre Krupa for his valuable support during my Ph.D. research. He gave me autonomy and valuable advices

Automatic centering in 2D

The centroid of the stiffest region of the elastogram, extracted using the prior process, will be used to horizontally center the rigid object in the full image. As the centroid coordinates are expressed in the ROI's frame (in RF frame units), it is necessary to express them in the ultrasound contact probe frame F cp (metric units).

ROBUST MOTION COMPENSATION

In the previous chapter, a palpation assistant system for ultrasound elastography was proposed to automatically move a 2D or 3D ultrasound probe to maintain the visibility of a stiff tissue of interest at the center of the FOV. This approach requires a processing step to segment and compute the centroid of the stiff tissue which is then automatically centered in the image by visual servoing. The system was designed for motionless tissues, and it gives good results under small perturbations. However, the previous implementation could not provide satisfactory output in case the tissue is moving since perturbation motion generates large noises in the estimated elastogram. This is also an important issue for a clinician when performing manual ultrasound elastography by moving the probe with his hand. Therefore, the main objective of this chapter is to present a method to estimate the strain map of a moving tissue, which was not possible with the process previously described. We propose a new robotic solution that exploits the intensity information of the 2D b-mode images. It is based on a tissue deformation tracking algorithm and an automatic 6-DOF compensation of the perturbation motion by ultrasound visual servoing.

In this chapter, a method for non-rigid motion estimation in 2D ultrasound images is proposed to estimate the displacement map required to compute the tissue strain map.

For this, the intensity changes in the ultrasound images due to the force applied by the probe are considered. Moreover, to estimate the strain map of a moving tissue, we propose to perform an automatic motion compensation using an ultrasound image-based visual servoing that synchronizes the probe and tissue motions during the strain map estimation process.

This chapter is divided in four sections. Section 3.1 introduces few works related where (s x , s y ) are the pixel width and height and (u cp , v cp ) are the pixel coordinates of the origin of the contact point frame, F cp , in the ultrasound 2D image.

In order to automatically compensate the moving tissue, we define the visual error was detailed in Chapter 2, and it is defined by the Equation (2.6). The fusion of the force control and the teleoperation task is achieved by using the redundancy control framework presented in Section 2.5. Therefore, the highest priority task is the force control defined as:

The secondary task is the teleoperation of the ultrasound probe. Then, to determine the second task, we must define a projector matrix P 1 that allows to project the second task onto the null space of the first task v 1 . The projector P 1 was defined in Equation (2.53), and the secondary task is expressed as

The secondary task v 2 does not disturb the primary task v 1 . Therefore, the fusion of the two robotic tasks leads to

where v p is now the velocity expressed at the frame F cp , and it includes the hierarchy of the two tasks. This velocity is applied to the robot using the Equations (2.7)-(2.9) in order to express them in the frame of the robot's end effector.

The control law presented in Equation (4.21) is always executed by the system. However, when the user switches to the impedance haptic mode in order to feel the force generated by the elastogram contained in the ROI, then we deactivate the secondary task by setting v 2 = 0 6×1 . This implies that the primary task is always running and the secondary task is paused so the elastogram can be estimated.

Experimental results

This section presents the results of the teleoperation and haptic feedback system previously described. First, we define the experimental setup as illustrated in Figure 4.14.

The haptic device used for the experiments is the Virtuose 6D (Haption S.A.). The robot and the ultrasound systems are the same used in the experiments presented in the previous chapters. We use the 3D ultrasound probe in 2D mode and the abdominal phantom ABDFAN US-1B. For more details about the equipments, please refer to Section 2.1.

The implementation of the system was coded in C++ and divided in five threads as shown in Figure 4.15. In the main thread is a graphic interface (coded with Qt libraries),

Résumé

Cette thèse se situe dans le contexte de la robotique médicale et porte sur l'élaboration d'un système robotisé permettant d'assister le processus d'élastographie ultrasonore. La solution proposée consiste à commander par retour d'effort et asservissement visuel un bras manipulateur actionnant une sonde échographique afin d'automatiser le mouvement de palpation nécessaire à la génération des images d'élasticité des tissus. La solution permet de réaliser l'imagerie élastographique de tissus sujets à des mouvements au moyen d'une tâche de compensation par asservissement visuel. Une approche innovante a également été proposée pour fournir à l'utilisateur un retour d'effort lui reflétant la sensation d'élasticité des tissus observés par l'intermédiaire d'un dispositif haptique. Les résultats expérimentaux des trois approches robotiques obtenus sur des fantômes constitués de tissus démontrent l'efficacité des méthodes proposées et ouvre des perspectives intéressantes pour l'élastographie ultrasonore assitée par robot.

Mots-clés -Robotique médicale, élastographie ultrasonore, asservissement visuel, haptique.

Abstract

This thesis concerns the development of a robotic control framework for quantitative ultrasound elastography. Ultrasound elastography is a technology that unveils elastic parameters of a tissue, which are commonly related with certain pathologies. This thesis proposes three novel robotic approaches to assist examiners with elastography. The first approach deals with the control of a robot actuating an ultrasound probe to perform palpation motion required for ultrasound elastography. The elasticity of the tissue is used to design a servo control law to keep a stiff tissue of interest in the field of view of the ultrasound probe. Additionally, the orientation of the probe is controlled by a human user to explore other tissue while elastography is performed. The second approach exploits deformable image registration of ultrasound images to estimate the tissue elasticity and to help in the automatic compensation by ultrasound visual servoing of a motion introduced into the tissue. The third approach offers a methodology to feel the elasticity of the tissue by moving a virtual probe in the ultrasound image with a haptic device while the robot is performing palpation motion. Experimental results of the three robotic approaches over phantoms with tissue-like offer an excellent perspective for robotic-assistance for ultrasound elastography.