# Data-Driven Interactive Simulations for Clinical and Robotic Applications

Jérémie Dequidt

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

## To cite this version:

Jérémie Dequidt. Data-Driven Interactive Simulations for Clinical and Robotic Applications. Modeling and Simulation. Université de Lille, 2023. tel-04144574

HAL Id: tel-04144574

https://inria.hal.science/tel-04144574

Submitted on 28 Jun 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Université de Lille

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

École Graduée Mathématiques, Sciences du numérique et de leurs interactions

Habilitation à Diriger des Recherches en Informatique

Simulation guidée par l'image pour
applications cliniques et robotiques

Data-Driven Interactive Simulations for Clinical
and Robotic Applications

Jérémie Dequidt

Composition du Jury

<table><tr><td>Rapporteur</td><td>Fernando Bello (PR)</td><td>Imperial College London</td></tr><tr><td>Rapporteur</td><td>Robert D. Howe (PR)</td><td>Harvard Paulson S.E.A.S</td></tr><tr><td>Rapporteuse</td><td>Maud Marchal (PR)</td><td>Université de Rennes</td></tr><tr><td>Examinatrice</td><td>Marie-Odile Berger (DR Inria)</td><td>Centre Inria Nancy - Grand Est</td></tr><tr><td>Examinateur</td><td>Pierre Renaud (PR)</td><td>INSA Strasbourg</td></tr><tr><td>Examinateur</td><td>Laurent Thines (PUPH)</td><td>CHU Besançon</td></tr><tr><td>Examinateur</td><td>Stéphane Cotin (DR Inria)</td><td>Centre Inria Nancy - Grand Est</td></tr><tr><td>President</td><td>Rochdi Merzouki (PR)</td><td>Université de Lille</td></tr><tr><td>Garant</td><td>Christian Duriez (DR Inria)</td><td>Inria Centre at the University of Lille</td></tr></table>

Vendredi, 10 février 2023
Jérémie Dequidt

Data-Driven Interactive Simulations for Clinical and Robotic Applications

Habilitation à Diriger des Recherches, Friday, February 10th 2023

Université de Lille

École Graduée Mathématiques, Sciences du numérique et de leurs interactions
Abstract

The research works briefly reviewed in this dissertation deal with questions related to the use of data in interactive simulations for applications in medicine and soft robotics. Geometrical modeling, (bio-)mechanical modeling as well as optimization techniques to reduce the computational footprint are presented. Using patient-specific or environment data with optimization techniques and simulation enables more realistic and predictive models and opens the way for new applications. Several examples are presented in the manuscript where our simulations were fed with data from medical images, videos streams or pre-operative data in order to provide more accurate medical simulations to provide better diagnoses or in the context of soft-robotics where visual servoing provides a more robust and stable control of soft-robots.
<figure><img src="image_4.png" /></figure>
Résumé

Les travaux de recherche examinés dans ce manuscrit traitent des questions liées à l'utilisation des données dans les simulations interactives pour des applications en médecine et en robotique déformable. La modélisation géométrique, la modélisation (bio-)mécanique ainsi que les techniques d'optimisation pour réduire l'empreinte de calcul sont présentées. L'utilisation de données spécifiques au patient ou à l'environnement avec des techniques d'optimisation et de simulation permet des modèles plus réalistes et prédictifs et ouvre la voie à de nouvelles applications. Plusieurs exemples sont présentés dans le manuscrit où nos simulations ont été alimentées avec des données provenant d'images médicales, de flux vidéo ou de données préopératoires afin de fournir des simulations médicales plus précises pour fournir de meilleurs diagnostics ou dans le contexte de la robotique déformable où l'asservissement visuel fournit un contrôle plus robuste et stable des robots souples.
<figure><img src="image_5.png" /></figure>
# Foreword

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

> Det spiller ingen rolle hvor sakte du går så lenge du ikke stopper.

> — Norwegian saying

First of all, I would like to thank my three reviewers: Maud, Fernando and Robert for accepted to review my work. I feel very honoured to have such prestigious scientists assessing my work and providing valuable feedbacks and inspiring comments about my work. I also would like to thank my defence committee (Laurent, Pierre, Marie-Odile, Stephane and Rochdi) for their nice words during the defence and for their interesting and relevant questions.

I also to thank Christian for his help during these past years. It has been a pleasure working with you since the Shaman/Shacra era and it has been a very passionate journey to help building DEFROST with you. I am very proud to be part of the team and to contribute to the many successes we had in medical simulations and soft robotics. Thank you also for being comprehensive about that the fact that an associate professor position implies many teaching and administrative duties.

Special thanks to Stephane: after my PhD, the perspective of having a career in academic was an unattainable destination and something I didn't consider for myself. You gave me the opportunity to join your team in Boston, you gave me full trust for developing my ideas and with you I unlocked my potential for doing research. We had many projects, many great achievements and many funny moments throughout the years. Even if things have changed since Mimesis / Defrost split, I am still very grateful for what you have done.

Many thanks for the numerous colleagues I have worked with since my PhD. I am taking advantage of the opportunity to thank the Alcove colleagues especially, Jowan and Jeremy which have now successful careers outside academics, the Shacra colleagues especially Nazim and Hugo (our Siggraph / post-Siggraph trip is still one of my best memories), Alex... and finally the Defrost ones Bruno, Damien, Alex K., Olivier, Eulalie, Thor, Quentin, Younes, Gang and the PhD students I have co-supervised. Thanks also for the people outside the Shacra/Defrost circle: Lambert, Marie-B, Axelle, Geraldine and of course thanks to our swiss-army knife / problem
solver wonder wom-Anne. Even if he has moved for sunnier places, I would like to thank David Simplot for being a role-model in many ways. Many thanks to the colleagues that are always ready for collaborative works: Erwan and P-F, Kenny and Miguel.

I would also like to thank colleagues from Polytech Lille because it has been an important place for my professional career and the office where I have spent most of my time lately (sorry Christian :/). Thanks to the teaching colleagues Alex, Florian, Nico, Sara, Bernard, Olivier, Othman, Rochdi, Laure... and François for sharing tips and tricks about managing an engineering department but also to the administrative colleagues Mumu, Melanie and Barbara for their efficient help and their responsiveness. It is nice working with all of you guys and having fun parties together. Special thanks to Manue (my partner in crime) for the fun and laugh you bring to this place and our fruitful collaboration. A final thank to Guy who learned me a lot, Nico for his trust (I hope to help you to do great things) and Cyril for his support.

Thank you also to my friends for their cheerful support and stress relief: Matthieu (my sidekick in sport), the Mexican mafia (Walter, Marcia, Vincent and Julien), the Houlette team (Oliv' and Caro, Steph and Marie, Jo and Eugen', Vincent and Sabine) and the polytech friends (Nour, Timmy, Pierre, Jean, Debo, Coline, Boris...) for the rare but precious moments. Thanks also to Clemence and Ju', Dom' and Louis, Pti' Pierre. Thanks also to the TNBM circle (Gylve, Ørjan, Vegard and Kristoffer), the DI6 circle (Douglas, Tony, Michael and Rose) and Lana for the inspiration.

And last but not least, I would like to thank my parents Gilles and Edith, my brother Geoffrey, my grand-mother Aniela and my family on the Putkownik and Dequidt side for their constant support and love. Thanks also to my step-family for their interest, support and accepting me as a full member of their clan: Bertrand, Françoise, Arnaud, Vincent and all the LSW. And finally a huge, full of love thank you to Marion for her love, her support and her inspiration that helps me reach this goal.
# Contents

<table><tr><td>Introduction</td><td>1</td></tr><tr><td>Curriculum Vitae</td><td>7</td></tr><tr><td>1 Interactive Physics-based Simulation</td><td>21</td></tr><tr><td>1.1 Motivations for Medical Procedures</td><td>22</td></tr><tr><td>1.2 Simulation Overview</td><td>24</td></tr><tr><td>1.3 1D Mechanical Model</td><td>28</td></tr><tr><td>1.3.1 Contributions</td><td>30</td></tr><tr><td>1.3.2 Coil Deployment in Collision-Free Environment</td><td>33</td></tr><tr><td>1.4 Asynchronous Simulation</td><td>36</td></tr><tr><td>1.5 Applications: Medical Simulators</td><td>38</td></tr><tr><td>1.5.1 Training System</td><td>39</td></tr><tr><td>1.5.2 Planning System</td><td>47</td></tr><tr><td>1.5.3 Discussion</td><td>52</td></tr><tr><td>2 Geometric Modeling for Interactive Simulation</td><td>53</td></tr><tr><td>2.1 Smooth surface using implicit surfaces</td><td>57</td></tr><tr><td>2.2 Models of vascular anatomy</td><td>63</td></tr><tr><td>2.2.1 Vessel Centerline Tracking</td><td>65</td></tr><tr><td>2.2.2 Implicit modeling: The LIM algorithm</td><td>70</td></tr><tr><td>2.2.3 Results</td><td>77</td></tr><tr><td>2.3 Models for the pelvic system</td><td>87</td></tr><tr><td>2.3.1 Overview of the approach</td><td>88</td></tr><tr><td>2.3.2 Results</td><td>96</td></tr><tr><td>2.3.3 Application: Pelvic System Mobility</td><td>100</td></tr><tr><td>3 Medical and Robotic Applications</td><td>113</td></tr><tr><td>3.1 Augmented Reality for Per-Operative Guidance</td><td>113</td></tr><tr><td>3.1.1 Mechanical Modeling</td><td>115</td></tr><tr><td>3.1.2 Non-rigid Registration</td><td>118</td></tr><tr><td>3.1.3 Implementation details</td><td>122</td></tr></table>
<table><tr><td>3.1.4 Experimental Results</td><td>123</td></tr><tr><td>3.2 Inverse FEM simulation for adaptive radiotherapy</td><td>123</td></tr><tr><td>3.2.1 Real-time Inverse FEM</td><td>125</td></tr><tr><td>3.2.2 Experimental validation</td><td>127</td></tr><tr><td>3.3 Motion Control for Interventional Radiology Procedures</td><td>129</td></tr><tr><td>3.3.1 Simulation with actuation and contacts</td><td>131</td></tr><tr><td>3.3.2 Control design</td><td>133</td></tr><tr><td>3.3.3 Experimental validation</td><td>136</td></tr><tr><td>3.4 Perspectives</td><td>137</td></tr><tr><td>Conclusion</td><td>139</td></tr><tr><td>Bibliography</td><td>143</td></tr></table>
# Introduction

> If the doors of perception were cleansed
> everything would appear to man as it is, infinite.
> For man has closed himself up till he sees all
> things thro' narrow chinks of his cavern.

> — William Blake

Digital technologies have greatly jolted modern medicine. For instance, the combina-
tion of electromagnetic radiation and digital geometry processing has originated the
*X-ray computed tomography* (or *X-ray CT*) that allows to create a three-dimensional
image of the targeted anatomical area. Another example is the use of computer
vision algorithms to fuse (*register*) data from different imaging modalities such as
*MRI* (*Magnetic Resonance Imaging*) that provides high-quality but slow-to-create
images and *Ultrasonography* (*US*) that delivers real-time but noisy images. This
fusion of information offers the opportunity to have high-quality and real-time visual
feedbacks of anatomical structures. These two examples are part of a larger set of
technologies such as *Radiography*, *Endoscopy*, *Positron Emission Tomography*... that
are usually referred by the generic term of **Medical Imaging**. There are many direct
benefits of Medical Imaging: it is non-invasive (no or small incision is required),
accurate, efficient, interactive and provides rich information. Location of diseased
areas is more accurately identified and tracked, evolution of pathologies is more
easily monitored allowing physicians to establish precise diagnoses. Initially aimed
at pre-operative steps of surgery, the availability of this technology has redefined
some surgical procedures where treatments can be delivered deep within the body
without surgery: for instance, bone biopsies can be performed in a minimally inva-
sive manner under image guidance, which avoid blood loss, complications, scars
and long recuperation [Jel+02]. They are now fully integrated in clinical routines
and have led to safer, less morbid and more cost-effective medical interventions.
The main drawback of such approach is that the complexity of the procedures has
significantly increased: anatomical structures manipulation is indirect through tele-
scopic devices or flexible ones and the clinicians observe the on-going procedure
through monitors. Several studies (for instance [Sil+09]) reports several disadvan-
tages related to the use of video-assisted surgery such as the loss of visual contact
with instruments that may cause disruptions or is more prone to complications, the usually fixed camera induces spatially transformed and suboptimal pictures and the loss of three-dimensional depth cues creates psychomotor problems of hand-eye coordination. Having inaccurate, suboptimal and non-consistent feedbacks (be they tactile, haptic or visual) adds a level of complexity to already stressful and complex procedures. To overcome the difficulties induced by minimally invasive surgery, two complementary approaches are considered. The first one is to make the surgeons more familiar with these procedures by providing a dedicated training that will first focus on the improvement of psychomotor skills. The second one is to provide computer or robotic assistance that will augment the surgeon capabilities to conduct the surgical procedures such as the daVinci® robot that allows to filter out hand tremors and to scale movements and forces.

Virtual Reality (VR) applications have been adequate training tools for improving hand-eye coordination in case of minimally invasive surgery because they offer a reproducible, controlled, supervised and ethical environment. The logical progression of these applications had been to train high-level tasks (such as suture) or even perform a complete surgical procedure by using simulation to mimic the behavior of the human anatomy. However, despite several governmental incentives [AC04] few clinical curricula require training on simulators. De Visser et. al. [De +11] report that while the benefits of VR simulators are indisputable to train novice-to-intermediate skills, they are yet not suitable for the training and maintenance of advanced skills. Insufficient level of physical realism, incomplete performance assessment and modest case complexity are suggested by De Visser et. al. [De +11] to be the main limitations to use simulators for surgical curricula. First, the insufficient level of physical realism relates to the fact that the modeling and simulation of human anatomy remains challenging. Several physical behaviors should be considered such as biomechanics, electro-physiology... and they are difficult to model. For instance, surveys in the field [Hum03] detail the complexities of the microstructure and the biomechanical behaviour of biological cells and tissues. These complexities can be captured, reproduced or sometimes predicted thanks to an important computation time. In the context of interactive VR applications, the trade-off between simulation accuracy versus simulation speed always leans toward the later producing fast-response applications but unable to faithfully reproduce the richness of the human anatomy. Second, limited case complexity can be partially related to the fact that the building of a virtual 3D model of an anatomical structure relies on many patient-specific parameters to be acquired. Some can be obtained through a time-consuming process which involves medical image segmentation (i.e. to delimit the different anatomical structures in an image), 3D reconstruction (i.e. to build a
3D model from segmented images) and 3D model tuning (i.e. trade-off between model complexity and desired accuracy) to be suitable for interactive applications. Some other parameters, for instance mechanical or physiological parameters, can be hardly acquired without invasive (and therefore ethically questionable) probes or sensors. As a consequence, these measurements are often conducted on animals or cadavers in order to have a rough estimation of the parameters. And some manual or semi-automatic post-processing can be needed to adjust these parameters as they are essential to obtain accurate simulations. Given the long time scale and the human interventions needed to acquire relevant data, VR simulators often propose a small set of virtual patients and scenarios. The length and the complexity of the acquisition pipeline do not only hinder the use of training simulator in the context of modern medicine but it also slows down the development of other applications. For instance, patient-specific interactive simulations could be used by surgeons to rehearse the procedure on a specific patient, to try different strategies, to face possible complications... in order to determine the safest procedure to be performed on the real patient.

Per-operative guidance could also be considered as a way to help the surgeon during an actual procedure by giving additional information such as the presence of weakened structures, the location of targeted area which are hardly visible... Augmenting the perception and the capabilities of the surgeons is the second approach used to reduce the complexity of minimally invasive procedures. In that case, patient-specific data plays a key role (and there are several other research problems to consider) for taking virtual reality applications to their next step: to provide, in real time, feedbacks about the geometry, mechanical or physiological status of the biological structures being manipulated. Moreover patient-specific data could also opens the way to new medical applications. For instance soft robotics, the scientific field that targets the design and the building robots made of soft materials, such as silicone, is expected to produce important outcomes in medicine. Surgical robots are made of metal alloys that interact with delicate biological tissues and are much stiffer than any anatomical structures (except for bones whose stiffness is of the same magnitude as metal alloys<fnref n="1" />). This stiffness gap results in robots damaging living tissues due to excessive contact pressures or stress that may lead to medical complications to the patient. Soft materials tighten this stiffness gap reducing potential damages when in contact with anatomical structures thanks to their ability to deform and to conform to the environment they are in contact with. This ability to deform raises however several scientific challenges; some are generic like how to accurately model

<footnote n="1">¹Stiffness of common materials can be found here: http://www.engineeringtoolbox.com/young-modulus-d_417.html</footnote>
and control a soft robot, some are specific to surgery where the control law of the
robot should take into account the possible contacts with the living tissues and thus
requires an accurate computational model of the human body.

These two complementary goals to help surgeons to perform surgery induce sev-
eral research problems for which we have proposed some original contributions.
This manuscript will be focused on medical applications since they were my main
application field since my PhD. But the contributions detailed here could be rele-
vantly used in other domains and several examples will be provided throughout this
document. Our aim is to tackle the problem of patient-specific data for interactive
simulations and the rationale of our approach is to mostly rely on images (be it
medical images or video streams) to estimate geometric and physical parameters.
We focus on images as the main source of patient-specific data for two reasons. First,
images are non-invasive and routinely used in medical procedures: we have already
mentioned that medical diagnosis is a big consumer of medical images (CT Scan,
IRM, UltraSound...), minimally invasive surgery (surgery through small incisions)
relies on tiny cameras or endoscope and finally interventional radiology procedures
are also supported by X-ray imaging. Second, images can be used as ground-truth to
validate our simulation and moreover can be also used in an optimization framework
to (dynamically) refine the parameters of the simulation. More specifically, the work
we have conducted over the last few years relates to the acquisition and optimization
(from a simulation standpoint) of patient-specific data, to the use of images to
validate experimentally simulation and to the use of images to infer or estimate
relevant parameters for the simulation. We want to emphasize that the scope of
our works lies at the intersection of Computer Vision, Geometric and Physics-based
Modeling and finally Interactive Simulation and that our main motivation is to
strive to enrich generic medical simulations to make them patient-specific. Bridging
the gap between simulation and images also implies several technical constraints
we have to deal with, such as the quality and availability of images or the limited
computational footprint of the physics-based model in order to achieve interactive
simulations. These constraints lead to not having a single solution relevant to every
kind of simulation or imaging techniques. However while the solutions are different
we want to illustrate that our methodology (i.e. the use of a fast physics-based
model as a backbone) allows to define relevant and efficient strategies for different
kind of medical procedures. Those works have thus been possible through tight
collaborations with surgeons and radiologists and have been declined in various
medical fields such as Interventional Radiology, Laparoscopic Liver Surgery, Pelvic
System Diagnosis and Radiotherapy. Finally modern medicine could benefit from
the use of soft robots. Soft Robotics is a recent cross-disciplinary field that binds
Applied Mathematics, Automatic Control, Mechanical Engineering, Robotics, Mecha-
tronics... In this field, the control of soft / deformable structures is a challenge and
we will demonstrate that the work that has been conducted on medical simulation
could be adapted to the control of soft robots.

The work presented in this manuscript covers the main results obtained following my
PhD in 2005. I first spend one year in the SimGroup (CIMIT/Harvard Medical School
/ Mass General Hospital, Cambridge, USA) in 2006–07 where I have been working
in geometrical modeling and contact modeling of vascular networks. This work
has been continued in 2007–08 when I joined a collaboration between Inria teams
Alcove and Magrit about the development of interventional radiology simulation.
Since 2008, I am Assistant Professor at University of Lille and I have been involved
in the research team Shaman / S.h.a.m.a.n / Shacra where the main focus was about
surgical simulations. There, I have worked on modeling, parameter-estimation of
digital models from medical images, simulation and augmented reality and the
natural trend of my research work has been to integrate more and more data
(medical images, video streams, sensors information...) in order to characterize
digital models and to enrich surgical simulations. Since 2015, I am a member of
Defrost research team whose main objective is to provide software tools to design,
simulate and control soft robots. In Defrost, I will extend the works on parameter
estimation to increase the robustness of the digital model of the soft-robot and to
enable a closed-loop control of the robot. The work presented here is a collective
effort conducted in about 10 years with colleagues in my successive research teams
and some fruitful collaborations. I had also the opportunity to advise some PhD
students: Ahmed Yureidini [Yur14], Nazim Haouchine [Hao15], Zhifan Jiang [Jia17],
Zhongkai Zhang [Zha19], Pierre Schegg [Sch22] and Flavie Przybylski who will
defend her thesis in 2024.

This manuscript is structured as follows. The first chapter [Interactive Physics-based Simulation] deals with interactive simulations: how to optimize the computational footprint and the possible use of digital simulations for medical applications. The work related to the extraction of dedicated, optimized geometric modeling from medical images is described in the second chapter [Geometric Modeling for Interactive Simulation]. The third chapter [Medical and Robotic Applications] details advances clinical and robotics applications that have been built using the works developed in the first two chapters.
<figure><img src="image_7.png" /></figure>
# Curriculum Vitae

## Jeremie Dequidt

<table><tr><td><strong>Nationality</strong></td><td>French</td><td><strong>Date of Birth</strong></td><td>December, 4th 1979</td></tr><tr><td><strong>Website</strong></td><td>www.dequidt.me</td><td><strong>Place of Birth</strong></td><td>Dechy (59), France</td></tr></table>

* Associate Professor (Maître de Conférences) at University of Lille, CNU Section 27 (Computer Science)

* Member the DEFROST team (joint team between CRISTAL and Inria)

* Head of the Embedded Systems Department at Polytech Lille

## Background

<table><tr><td>2008-now</td><td><strong>Associate Professor, University of Lille</strong><br><em>Research:</em> Geometrical and Biomechanical Modeling, Animation and Digital Simulation, Virtual and Augmented Reality, Soft Robotics<br><em>Teaching:</em> Programming, Computer Theory, Computer Graphics, Numerical Analysis, Digital Physics-based Simulation</td></tr><tr><td>2007-2008</td><td><strong>Post Doctorate, Inria</strong><br><em>Research:</em> Mechanical modeling of surgical tools for the simulation of interactive radiology interventions</td></tr><tr><td>2006-2007</td><td><strong>Research Fellow, CIMIT / Harvard Medical School, Cambridge (MA)</strong><br><em>Research:</em> Mechanical modeling of surgical tools for the simulation of interactive radiology interventions</td></tr><tr><td>2005-2006</td><td><strong>Adjunct Professor, Polytech Lille</strong><br><em>Research:</em> Computer Graphics and Simulation, Distributed Computation<br><em>Teaching:</em> Operating Systems, Databases, Programing</td></tr><tr><td>2002-2005</td><td><strong>PhD Student, University of Lille</strong><br><em>Research:</em> Computer Graphics and Simulation, Distributed Computation<br><em>Teaching:</em> Operating Systems, Networks, Programing</td></tr></table>
# Doctoral Co-Supervision

## Defended Thesis

* Ahmed Yureidini: Robust Blood Vessel Modeling for Interactive Simulations, defended 13/05/2014

* Nazim Haouchine: Image-guided simulation for augmented reality during hepatic surgery, defended 26/01/2015

* Zhifan Jiang: Évaluation des mobilités et modélisation géométrique du système pelvien féminin par analyse d'images médicales, defended 27/01/2017

* Zhongkai Zhang: Vision-based Calibration, Position Control and Force Sensing for Soft Robots, defended 10/01/2020

* Pierre Schegg: Autonomous Catheter and Guidewire Navigation for Robotic Percutaneous Coronary Interventions, defended 25/05/22

## On-going Thesis

* Flavie Przybylski (start 1/10/2021): Simulation, design and fabrication of wine-robots

## Undefended Thesis

* Francois Dervaux (start 1/10/2012), has leaved 18 months after the beginning of the thesis for a job in the industry

## Post-Doctorate and Master Supervision

### Post-Doctorate

* 2017-18, Thomas Morzadec, geometric modeling of soft-robots
## Master Students

* 2018, Meichun Lin (Bio Medical Engineering Master2), training simulator for brachytherapy procedures

* 2017, Félix Vanneste (ISEN Engineer), bridging between SOFA and ROS

* 2016, Piotr Sordyl (Erasmus Master), sensibility analysis of catheter simulations

* 2014, Matthias De Bie (Polytech Lille Engineer), bridging rendering and simulation softwares

* 2013, Pierre-Jean Petitprez (Polytech Lille Engineer), real-time X-ray rendering

* 2013, Mélanie Lelaure (Polytech Lille Engineer), user interface for interventional radiology simulations

* 2012, François Dervaux (Master2 IVI), asynchronous simulation in SOFA

* 2011, Nazim Haouchine (Master2 IVI), physics-based simulation and augmented reality

## Temporary Engineers

* 2018, Meichun Lin, training simulator for brachytherapy procedures

* 2016, Christo Gnonnou, sensibility analysis of catheter simulations

* 2015, Mario Sanz Lopez, device for tracking neuro-surgical devices

## Highlights

* Recipient of the Young Researcher Funding Programme (French Research Agency) in 2012, 175 keuros

* Recipient of the Dirk Bartz Prize for Visual Computing in Medicine in 2015

* Recipient of the Doctoral Supervision and Research Bonus (PEDR) for 2015–2018 and 2019–2022

* Runner-up Best Paper Award in ISMAR 2013

* TVCG Top Ranked AR/VR paper (presented at SIGGRAPH 2015)

* Best Paper Award in SIMPAR 2016

* Best Paper Award in Advanced Robotics 2018
# Collaborations and Responsibilities

## Collaborations

* University of Copenhagen, Kenny Erleben and Sune Darkner (on-going project about Digital Twins)

* Universidad Rey Juan Carlos, Miguel Otaduy (on-going project about Digital Twins)

* Vrije Universiteit Brussel, Bram Vanderborght (on-going project about Digital Twins)

* Magrit Team (Inria Nancy), Erwan Kerrien and Marie-Odile Berger (2 PhD co-supervising)

* Mimesis Team (Inria Nancy), Stéphane Cotin (2 PhD co-supervising), Igor Peterlik (publications on asynchronous simulation)

* LML, Mathias Brieu, Pauline Leconte, Jean-Francois Witz (1 PhD co-supervising)

* COL, Nick Reynaert (publication on deformable registration using FEM)

* King's College London, Hongbin Lu (publication on visual servoing for robotized catheter)

* CI2S (CRISTAL), Rochdi Merzouki (project on robotized brachytherapy)

* CHR Huriez, Jean-Francois Rouland (various publications on cataract surgery simulation)

## Implication in the scientific communities

### International Program Committee

* International Symposium on Visual Computing 2018, 2019, 2020

* Workshop on Virtual Reality Interaction and Physical Simulation (VRIPHYS) 2013, 2014, 2015 (+ Co-Organizer), 2018

* International Symposium on Biomedical Simulation 2010 and 2014 (+ General Chair)
## Projects

* Member of the CNRS project Visadapt (PEPS 2021-2023

* Member of the European project CoBra (Interreg 2 Seas) 2018 - 2022

* Member of the ADT (Inria Software Development Initiative) SIMILAR, 2018 - 2020

* Member of the European project COMOROS (FERDER grant), 2017 - 2020

* Member of the ANR project SORCERY, 2017-2020

* Co of the ADT (Inria Software Development Initiative) SOROS, 2016 - 2018

* Prime Investigator of the ANR project IDeaS (Image-Driven Simulation), 2012 - 2016, 175 keuros

* Member of the ADT (Inria Software Development Initiative) SORBET, 2014 - 2016

* Member of the ADT (Inria Software Development Initiative) SIMPLE, 2007 - 2009

## Technology Transfers

* 2018: Technology transfer of an interventional radiology simulator to Navworks

* 2015: Technology transfer of a neurosurgical simulator to InSimo

* 2013: Technology transfer of a cataract surgery (MSICS) simulator to InSimo

* 2011: Industrial contract with HelpMeSee foundation (NGO)

## Expertise

* Vice-President of the hiring committee for Assistant Professor positions (Polytech Lille) 2021, 2022

* Member of the evaluation committee for the french research agency (CE33: Interaction Robotique) 2021-2023

* Reviewer for french research agency (2017 generic call)

* Member of external hiring committees for Assistant Professor positions Calais (2010, 2011), Poitiers (2011)

* Elected Member of the internal hiring committee for Assistant Professor positions (since 2011)

* Reviewer for Soft Robotics, MICCAI, Media, Computer Graphics and Applications, IEEE Transaction on Haptics, IEEE/RSJ IROS, IEEE World Haptics, IEEE Transactions on Medical Imaging, Computer Graphics
## Teachings

* Head of the Computer and Electrical Engineering (IMA) Department at Polytech-Lille (since 2017-2022)

* In charge of the coordination of the Cross-Disciplinary modules at Polytech-Lille since 2015 (200 students)

* Coordinator of Ima SC (Communicating Systems) track at Polytech-Lille (2011-2017)

* Coordinator of Ima5i track (from 2008 to 2011... now replaced by IMA SC track)

* In charge of 3 modules in computer science at Polytech Lille

* Active participant to the teaching reforms at Polytech Lille

* Member of the Administration board at Polytech Lille

* Teaching charge of 300h/year

## Personal Publications

[Anx+09] René Anxionnat, Federico Rocca, Serge Bracard, Jeremie Dequidt, Erwan Kerien, Christian Duriez, Marie-Odile Berger, and Stephane Cotin. "Evaluation of a computer-based simulation for the endovascular treatment of intracranial aneurysms". In: *10th congress of the World Federation of Interventional and Therapeutic Neuroradiology*. 2009.

[Ben+13] Jan Bender, Jeremie Dequidt, Christian Duriez, and Gabriel Zachmann, eds. *VRIPHYS 13: 10th Workshop on Virtual Reality Interactions and Physical Simulations, Lille, France, 2013. Proceedings*. Eurographics Association, 2013.

[Bil+11] Alexandre Bilger, Jeremie Dequidt, Christian Duriez, and Stephane Cotin. "Biomechanical simulation of electrode migration for deep brain stimulation". In: *Medical Image Computing and Computer-Assisted Intervention*. 2011, pp. 339-346.

[Bos+14] Julien Bosman, Nazim Haouchine, Jeremie Dequidt, Igor Peterlik, Stephane Cotin, and Christian Duriez. "The Role of Ligaments: Patient-Specific or Scenario-Specific?" In: *6th International Symposium on Biomedical Simulation (ISBMS)*. 2014, pp. 228-232.
[Bou+10] Nadia Boubchir, Stephane Cotin, Olivier Comas, Frederick Roy, Christian Duriez, Jeremie Dequidt, Jeremie Allard, and Jean-Francois Rouland. "Computer-Based Simulation of Cataract Surgery: Toward a New Teaching Paradigm". In: *Annual Meeting of the American Academy of Ophthalmology*. 2010.

[Bou+09] Nadia Boubchir, Stephane Cotin, Christian Duriez, Jeremie Dequidt, Jeremie Allard, and Jean-Francois Rouland. "Computer-Based Simulation of IOL Injection: Toward a Full Featured Cataract Surgery Training System". In: *Annual Meeting of the American Academy of Ophthalmology*. 2009 (cit. on p. 43).

[Coe+17] Eulalie Coevoet, Thor Morales-Bieze, Frederick Largilliere, et al. "Software toolkit for modeling, simulation, and control of soft robots". In: *Advanced Robotics* 31.22 (2017). Best Paper Award, pp. 1208-1224.

[Coe+14] Eulalie Coevoet, Nick Reynaert, Eric Lartigau, Luis Schiappacasse, Jeremie Dequidt, and Christian Duriez. "Introducing interactive inverse FEM simulation and its application for adaptive radiotherapy". In: *Medical Image Computing and Computer Assisted Intervention*. 2014, pp. 81-88.

[Coe+15] Eulalie Coevoet, Nick Reynaert, Eric Lartigau, Luis Schiappacasse, Jeremie Dequidt, and Christian Duriez. "Registration by interactive inverse simulation: application for adaptive radiotherapy". In: *International Journal of Computer Assisted Radiology and Surgery* 10.8 (2015), pp. 1193-1200.

[Cot+05] Stephane Cotin, Paul Neumann, Xunlei Wu, Sylvère Fonteneau, Pierre-Jean Bensoussan, Damien Marchal, Jeremie Dequidt, Laurent Grisoni, and Sylvain Karpf. "Collaborative Development of an Open Framework for Medical Simulation". In: *Insight Journal*. 2005.

[Deq02] Jeremie Dequidt. "Détection de collisions entre objets physiques autonomes". MA thesis. University of Lille 1, 2002.

[Deq17] Jeremie Dequidt. "Modeling interactions with soft objects". In: *Workshop on Soft Morphological Design for Haptic Sensation, Interaction and Display*. 2017.

[Deq05] Jeremie Dequidt. "Objets autonomes en simulation physique temps-réel". PhD thesis. University of Lille 1, 2005.
[Deq+15] Jeremie Dequidt, Eulalie Coevoet, Laurent Thines, and Christian Duriez. "Vascular neurosurgery simulation with bimanual haptic feedback". In: *12th Workshop on Virtual Reality Interaction and Physical Simulation (VRIPHYS)*. Nov. 2015.

[Deq+13] Jeremie Dequidt, Hadrien Courtecuisse, Olivier Comas, Jeremie Al- lard, Christian Duriez, Stephane Cotin, Elodie Dumortier, Olivier Wavreille, and Jean-Francois Rouland. "Computer-based training sys- tem for cataract surgery". In: *Simulation: Transactions of the Society for Modeling and Simulation International* 16 (2013), pp. 1–15.

[Deq+09] Jeremie Dequidt, Christian Duriez, Stephane Cotin, and Erwan Ker- rien. "Towards interactive planning of coil embolization in brain aneurysms". In: *Medical Image Computing and Computer-Assisted In- tervention*. 2009, pp. 377–385.

[DGC04] Jeremie Dequidt, Laurent Grisoni, and Christophe Chaillou. Asyn- chronous Interactive Physical Simulation. Tech. rep. 5338. ALCOVE (INRIA / LIFL), 2004.

[DGC05a] Jeremie Dequidt, Laurent Grisoni, and Christophe Chaillou. "Collabo- rative interactive physical simulation". In: *GRAPHITE '05: Proceedings of the 3rd international conference on Computer graphics and interac- tive techniques in Australasia and South East Asia*. 2005, pp. 147–150 (cit. on p. 37).

[DGC05b] Jeremie Dequidt, Laurent Grisoni, and Christophe Chaillou. "Time- Critical Animation of Deformable Solids". In: *Computer Animation and Social Agents*. 2005.

[DGM04] Jeremie Dequidt, Laurent Grisoni, and Philippe Meseure. "Détection de collisions entre objets rigides convexes autonomes". In: *Revue de CFAO et d'informatique graphique* 18.2 (2004), pp. 183–195.

[DLC07] Jeremie Dequidt, Julien Lenoir, and Stephane Cotin. "Interactive contacts resolution using smooth surface representation". In: *Medical Image Computing and Computer-Assisted Intervention*. 2007, pp. 850– 857.

[DMG05] Jeremie Dequidt, Damien Marchal, and Laurent Grisoni. "Time-Critical Animation of Deformable Solids". In: *Computer Animation and Virtual Worlds* 16 (2005), pp. 177–187.
[Deq+08] Jeremie Dequidt, Maud Marchal, Christian Duriez, Erwan Kerien, and Stephane Cotin. "Interactive simulation of embolization coils: Modeling and experimental validation". In: *Medical Image Computing and Computer-Assisted Intervention*. 2008, pp. 695-702.

[Der+13] Francois Dervaux, Igor Peterlik, Jeremie Dequidt, Jeremie Cotin, and Duriez Christian. "Haptic Rendering of Interacting Dynamic Deformable Objects Simulated in Real-Time at Different Frequencies". In: *IEEE/RSJ International Conference on Intelligent Robots and Systems*. 2013, pp. 2010-2016 (cit. on p. 37).

[Dum+10] Elodie Dumortier, Stephane Cotin, Jeremie Dequidt, Christian Duriez, Jeremie Allard, and Jean-Francois Rouland. "Computer-Based Training System for Cataract Surgery". In: *International ICST Conference on Electronic Healthcare for the 21st Century*. 2010.

[Dum+11a] Elodie Dumortier, Stephane Cotin, Jeremie Dequidt, and Jean-Francois Rouland. "A Prototype of Simulation System for Cataract Surgery Training". In: *European Society of Cataract and Refractive Surgeons Winter Meeting*. 2011.

[Dum+11b] Elodie Dumortier, Stephane Cotin, Jeremie Dequidt, and Jean-Francois Rouland. "Assessment Metrics For A Prototype Of Simulation System For Cataract Surgery Training". In: *Annual Meeting of the Association for Research in Vision and Ophthalmology*. 2011.

[Dum+11c] Elodie Dumortier, Stephane Cotin, Jeremie Dequidt, and Jean-Francois Rouland. "Realisation d'un Prototype de Simulateur Fonctionnel pour la Chirurgie de la Cataracte". In: *117e CONGRES de la Société Française d'Ophtalmologie*. 2011.

[Dur+16] Christian Duriez, Eulalie Coevoet, Frederick Largilliere, T Morales-Bieze, Zhongkai Zhang, Mario Sanz-Lopez, B Carrez, Damien Marchal, Olivier Goury, and Jeremie Dequidt. "Framework for online simulation of soft robots with optimization-based inverse model". In: *Simulation, Modeling, and Programming for Autonomous Robots (SIMPAR), IEEE International Conference on*. Best Paper Award. 2016, pp. 111-118.

[Gou+16] Olivier Goury, Yann Nguyen, Renato Torres, Jeremie Dequidt, and Christian Duriez. "Numerical Simulation of Cochlear-Implant Surgery: Towards Patient-Specific Planning". In: *International Conference on Medical Image Computing and Computer-Assisted Intervention (MIC-CAI)*. 2016, pp. 500-507.
[GD05] Laurent Grisoni and Jeremie Dequidt. "Asynchronous Interactive Physical Simulation". In: *1st ICT Workshop on Virtual Reality*. 2005.

[GDC04] Laurent Grisoni, Jeremie Dequidt, and Christophe Chaillou. "Virtual Reality Applications and Frame-Rate Control". In: *Development and Application Systems*. 2004, pp. 461-465.

[Hao+15a] Nazim Haouchine, Alexandre Bilger, Jeremie Dequidt, and Stephane Cotin. "Fracture in Augmented Reality". In: *ACM SIGGRAPH 2015 Posters*. SIGGRAPH '15. Los Angeles, California: ACM, 2015, 22:1-22:1.

[Hao+14a] Nazim Haouchine, Stephane Cotin, Jeremie Dequidt, Erwan Kerrien, Berger, and Marie-Odile. "Réalité augmentée pour la chirurgie minimalement invasive du foie utilisant un modèle biomécanique guidé par l'image". In: *Actes de la conference RFIA 2014*. 2014.

[Hao+15b] Nazim Haouchine, Stephane Cotin, Igor Peterlik, Jeremie Dequidt, Mario Sanz Lopez, Erwan Kerrien, and Marie-Odile Berger. "Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery". In: *IEEE Transactions on Visualization and Computer Graphics* 21.5 (May 2015), pp. 584-597.

[Hao+13a] Nazim Haouchine, Jeremie Dequidt, Marie-Odile Berger, and Stephane Cotin. "Deformation-based Augmented Reality for Hepatic Surgery". In: *Medicine Meets Virtual Reality*. 2013, pp. 182-188.

[Hao+15c] Nazim Haouchine, Jeremie Dequidt, Marie-Odile Berger, and Stephane Cotin. "Monocular 3D Reconstruction and Augmentation of Elastic Surfaces with Self-Occlusion Handling". In: *IEEE Transactions on Visualization and Computer Graphics* 21.12 (Dec. 2015), pp. 1363-1376.

[Hao+14b] Nazim Haouchine, Jeremie Dequidt, Marie-Odile Berger, and Stephane Cotin. "Single View Augmentation of 3D Elastic Objects". In: *Mixed and Augmented Reality (ISMAR), 2014 IEEE International Symposium on*. 2014, pp. 229-236.

[Hao+12] Nazim Haouchine, Jeremie Dequidt, Erwan Kerrien, Marie-Odile Berger, and Stephane Cotin. "Physics-based Augmented Reality for 3D Deformable Object". In: *VRIPHYS - Virtual Reality Interaction and Physical Simulation*. 2012, pp. 31-38.

[Hao+13b] Nazim Haouchine, Jeremie Dequidt, Igor Peterlik, Erwan Kerrien, Marie-Odile Berger, and Stephane Cotin. "Image-guided Simulation of Heterogeneous Tissue Deformation For Augmented Reality during Hepatic Surgery". In: *Mixed and Augmented Reality (ISMAR), 2013*
IEEE International Symposium on. Runner-up Best Paper Award. 2013, pp. 199–208.

[Hao+14c] Nazim Haouchine, Jeremie Dequidt, Igor Peterlik, Erwan Kerrien, Marie-Odile Berger, and Stephane Cotin. "Towards an accurate tracking of liver tumors for augmented reality in robotic assisted surgery". In: Robotics and Automation (ICRA), 2014 IEEE International Conference on. 2014, pp. 4121–4126.

[JDC15] Zhifan Jiang, Jeremie Dequidt, and Stephane Cotin. "Détection de contours des organes pelviens dans des images médicales par modèles de B-spline". In: Congrès Français de Mécanique. AFM, Association Française de Mécanique, 2015.

[Jia+14] Zhifan Jiang, Witz Jean-Francois, Lecomte-Grosbras Pauline, Jeremie Dequidt, Stephane Cotin, Christian Duriez, and Mathias Brieu. "Segmentation, Geometric Modeling And Motion Tracking For Evaluation Of Pelvic System Mobilities". In: International Conference on Experimental Mechanics. 2014.

[Jia+17] Zhifan Jiang, Jean-Francois Witz, Pauline Lecomte-Grosbras, Jeremie Dequidt, Michel Cosson, Stephane Cotin, Christian Duriez, and Mathias Brieu. "Multi-organ Motion Tracking in Dynamic Magnetic Resonance Imaging for Evaluation of Pelvic System Mobility and Shear Strain". In: Strain 53.2 (2017) (cit. on p. 88).

[Jia+18] Zhifan Jiang, Jean-Francois Witz, Pauline Lecomte-Grosbras, Jeremie Dequidt, Michel Cosson, Stephane Cotin, Christian Duriez, and Mathias Brieu. "Virtual Image Correlation of Magnetic Resonance Images for 3D Geometric Modelling of Pelvic Organs". In: Strain 55.3 (2018), pp. 1230–5 (cit. on p. 88).

[Jia+15] Zhifan Jiang, Jean-Francois Witz, Pauline Lecomte-Grosbras, Jeremie Dequidt, Christian Duriez, Michel Cosson, Stephane Cotin, and Mathias Brieu. "B-spline Based Multi-organ Detection in Magnetic Resonance Imaging". In: Strain 51.3 (2015), pp. 235–247 (cit. on p. 88).

[KBD08] Erwan Kerrien, Marie-Odile Berger, and Jeremie Dequidt. "Refining the 3D surface of blood vessels from a reduced set of 2D DSA images". In: Augmented environments for Medical Imaging and Computer-aided Surgery. 2008, pp. 61–69.

[Ker+17] Erwan Kerrien, Ahmed Yureidini, Jeremie Dequidt, Christian Duriez, Rene Anxionnat, and Stephane Cotin. "Blood vessel modeling for
<table>
  <tr>
    <td>[Lar+15]</td>
    <td>Frederick Largilliere, Valerian Verona, Eulalie Coevoet, Mario Sanz-Lopez, Jeremie Dequidt, and Christian Duriez. "Real-time control of soft-robots using asynchronous finite element modeling". In: Robotics and Automation (ICRA), 2015 IEEE International Conference on. May 2015, pp. 2550-2555.</td>
  </tr>
  <tr>
    <td>[Lec+13]</td>
    <td>Pauline Lecomte-Grosbras, Mouhamadou Nassirou Diallo, Jean-Francois Witz, Damien Marchal, Jeremie Dequidt, Stephane Cotin, Michel Cosson, Christian Duriez, and Mathias Brieu. "Towards a better understanding of pelvic system disorders using numerical simulation". In: Medical Image Computing and Computer-Assisted Intervention. 2013, pp. 307-314.</td>
  </tr>
  <tr>
    <td>[Mar+06]</td>
    <td>Damien Marchal, Fabrice Aubert, Christophe Chaillou, Laurent Grisoni, and Jeremie Dequidt. "Towards Cuttable Deformable Objects". In: 2nd ICT Workshop on Virtual Reality. 2006.</td>
  </tr>
  <tr>
    <td>[Rab+07]</td>
    <td>James Rabinov, Stephane Cotin, Jeremie Allard, Jeremie Dequidt, Julien Lenoir, Vincent Luboz, Paul Neumann, Xunlei Wu, and Steve Dawson. "EVE: Computer Based Endovascular Training System for Neuroradiology". In: 45th Annual Meeting &amp; NER Foundation Symposium. 2007, pp. 147-150.</td>
  </tr>
  <tr>
    <td>[San+14]</td>
    <td>Mario Sanz-Lopez, Jeremie Dequidt, Erwan Kerrien, Christian Duriez, Marie-Odile Berger, and Stephane Cotin. "Testbed for Assessing the Accuracy of Interventional Radiology Simulations". In: 6th International Symposium on Biomedical Simulation (ISBMS). 2014, pp. 106-111.</td>
  </tr>
  <tr>
    <td>[Sch+22]</td>
    <td>Pierre Schegg, Jeremie Dequidt, Eulalie Coevoet, Edouard Leurent, Remi Sabatier, Philippe Preux, and Christian Duriez. "Automated Planning For Robotic Guidewire Navigation In The Coronary Arteries". In: IEEE 5th International Conference on Soft Robotics. 2022 (cit. on p. 140).</td>
  </tr>
  <tr>
    <td>[SDG12a]</td>
    <td>David Selosse, Jeremie Dequidt, and Laurent Grisoni. "A Sketch-Based Interface for Annotation of 3D Brain Vascular Reconstructions". In: Computer Graphics International. 2012.</td>
  </tr>
  <tr>
    <td>[SDG12b]</td>
    <td>David Selosse, Jeremie Dequidt, and Laurent Grisoni. A Sketch-Based Interface for Annotation of 3D Brain Vascular Reconstructions. Tech. rep. 00698832. SHACRA (INRIA LNE / LIFL), 2012.</td>
  </tr>
</table>
[Tal+15] Hugo Talbot, Nazim Haouchine, Igor Peterlik, Jeremie Dequidt, Christian Duriez, Herve Delingette, and Stéphane Cotin. "Surgery Training, Planning and Guidance Using the SOFA Framework". In: *EG 2015 - Dirk Bartz Prize*. Ed. by H.-C. Hege and T. Ropinski. prize. The Eurographics Association, 2015.

[Wei+12] Yiyi Wei, Stephane Cotin, Jeremie Dequidt, Christian Duriez, Jeremie Allard, and Erwan Kerrien. "A (Near) Real-Time Simulation Method of Aneurysm Coil Embolization". In: *Aneurysm*. InTech Publisher, 2012.

[Yur+11] Ahmed Yureidini, Jeremie Dequidt, Erwan Kerrien, Christian Duriez, and Stephane Cotin. "Computer-based simulation for the endovascular treatment of intracranial aneurysms". In: *Living Imaging (Imaging and its applications in medicine and life sciences)*. 2011.

[Yur+12] Ahmed Yureidini, Erwan Kerrien, Jeremie Dequidt, Christian Duriez, and Stephane Cotin. "Local Implicit Modeling of Blood Vessels for Interactive Simulation". In: *Medical Image Computing and Computer-Assisted Intervention*. 2012, pp. 553-560.

[Zha+19a] Zhongkai Zhang, Jeremie Dequidt, Junghwan Back, Hongbin Liu, and Christian Duriez. "Motion Control of Cable-Driven Continuum Catheter Robot through Contacts". In: *Robotics and Automation (ICRA), 2019 IEEE International Conference on*. 2019.

[Zha+19b] Zhongkai Zhang, Jeremie Dequidt, Junghwan Back, Hongbin Liu, and Christian Duriez. "Motion Control of Cable-Driven Continuum Catheter Robot through Contacts". In: *IEEE Robotics and Automation Letters* 4.2 (2019), pp. 1852-1859 (cit. on p. 135).

[ZDD18a] Zhongkai Zhang, Jeremie Dequidt, and Christian Duriez. "Vision-based Sensing of External Force for Soft Robots using Finite Element Model". In: *Robotics and Automation (ICRA), 2018 IEEE International Conference on*. 2018, pp. 1529-1536.

[ZDD18b] Zhongkai Zhang, Jeremie Dequidt, and Christian Duriez. "Vision-based Sensing of External Force for Soft Robots using Finite Element Model". In: *IEEE Robotics and Automation Letters* 3.3 (2018), pp. 1529-1536 (cit. on p. 129).

[Zha+16] Zhongkai Zhang, Jeremie Dequidt, Alexandre Kruszewski, Frederick Largilliere, and Christian Duriez. "Kinematic Modeling and Observer Based Control of Soft Robot using Real-Time Finite Element Method". In: *IEEE/RSJ International Conference on Intelligent Robots and Systems*. 2016, pp. 5509-5514 (cit. on p. 129).
[Zha+17] Zhongkai Zhang, Thor Enrique Morales Bieze, Jeremie Dequidt, Alexandre Kruszewski, and Christian Duriez. "Visual Servoing Control of Soft Robots Based on Finite Element Model". In: *IEEE/RSJ International Conference on Intelligent Robots and Systems*. 2017, pp. 2895-2901 (cit. on p. 129).

[Zha+19c] Zhongkai Zhang, Antoine Petit, Jeremie Dequidt, and Christian Duriez. "Calibration and External Force Sensing for Soft Robots using an RGB-D Camera". In: *Robotics and Automation (ICRA), 2019 IEEE International Conference on*. 2019.

[Zha+19d] Zhongkai Zhang, Antoine Petit, Jeremie Dequidt, and Christian Duriez. "Calibration and External Force Sensing for Soft Robots using an RGB-D Camera". In: *IEEE Robotics and Automation Letters* 4.3 (2019), pp. 2356-2363.

[Zha+22] Zhongkai Zhang, Antoine Petit, Jeremie Dequidt, and Christian Duriez. "SofaGym: An open platform for Reinforcement Learning based on Soft Robot simulations". In: *Soft Robotics Journal TBA.TBA* (20x22), TBA (cit. on p. 140).
Interactive Physics-based
Simulation

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

Dissimuler est feindre de ne pas avoir ce qu'on a. Simuler est feindre d'avoir ce qu'on n'a pas. L'un renvoie à une présence, l'autre à une absence.

— Jean Baudrillard

Simulacres et Simulation

Physics-based simulation or computer simulation<fnref n="1" /> is a standard tool for engineers and scientists. Thanks to the availability of computing technologies and their increasing performance, it has been widely used in many domains for the design and improvement of complex systems or to understand / predict natural phenomenon. For these kind of applications, the main target is accuracy because users want to be able to predict precisely what will happen in the real world. This accuracy is achieved by using complex physics models where the governing equations are highly non-linear and the parameters are numerous preventing to obtain an analytical solution of the problem. Numerical solutions are therefore sought after by discretizing the phenomenon in space and time leading to a high number of unknowns to solve. As a consequence, typical simulations in engineering exhibit high computation time or require clusters of computer to provide results in a reasonable time [Ste13].

Over the past decade I have been focusing my efforts in using computer simulation for medical applications. While scientific deadlocks related to general simulation still exist, modeling and simulating living tissues or organs raises specific challenges. First, the characterization of living materials is complex as few methods exist to perform harmless, ethical, non-invasive, in-vivo probing. Elastography [WMK06], an imaging technique that relies on the Doppler Effect to provide tissue movements or strain, may become the standard technology for biomechanical characterization but technological improvements are still needed to be effective. Second, properties of living materials exhibit lot of intra-variability (property is not constant in the whole material) and inter-variability (property is not constant for all beings). Third, experimental validation is also a challenging task and shares the same difficulties

<footnote n="1">¹https://en.wikipedia.org/wiki/Computer_simulation</footnote>
than characterization: the scarcity and the complexity to acquire data. To these
intrinsic challenges, I also wanted that these computer simulations target specific
medical applications: training simulators and computer-assisted software for surgery.
Therefore this raises the need for interactive simulation that allows virtual interaction
of a user on the simulation. Because this user can move, probe or manipulate the
virtual models, this leads to unexpected solicitations on the model that cannot be
pre-computed. Moreover the user expects the virtual simulation to immediately
react after each of his/her actions as it will happen in a real surgical procedure. As a
consequence, bounding computation time is the main has been another issue I have
tackled in my research works.

The intent of the previous paragraphs was to provide a *big picture* of the large scope
of scientific deadlocks that are related to computer simulation. My research activities
in the field has been devoted to a much narrower spectrum. First of all, I have
been mainly interested in mechanical / biomechanical modeling. Because most of
living materials are soft and because some surgical devices are also soft, I have been
working on the modeling on soft materials and the modeling of interactions between
soft materials and rigid ones. Second to reduce the computational footprint or to
provide the optimal trade-off between simulation accuracy and speed, I have worked
on asynchronous strategies that ensure mechanical consistency while being less
demanding with respect to computation power. Finally fast and accurate computer
simulations require some work on software / hardware integration. The remaining
of this chapter will present some scientific contributions on these different items
to illustrate my work. It is not intended to be exhaustive as the reader will find a
complete list of publications related to computer modeling and simulation at the
end of the manuscript.

## 1.1 Motivations for Medical Procedures

Until the 1960's, Radiology was only used as a way to diagnose vascular diseases
through angiography. In January 16, 1964, the first intervention (a dilatation of
a stenosis) through Radiology was performed by Charles Dotter [RKK03]. *Interventional Radiology* became popular after that first initiative as an efficient way to
treat vascular disease. It involves the insertion of thin flexible devices (*catheters* or *guidewire*) than are pushed in the vascular system until the reach of the targeted
area and many pathologies can now be treated through this technique.
<figure><img src="image_9.png" /><figcaption>**Fig. 1.1:** Photography of a radiology intervention: thin flexible devices are inserted in the femoral artery and are navigated inside the vascular network until reaching the diseased area. Then different treatments can be performed. Monitoring the motion of the device in the patient is performed via X-ray imaging. Photography appears courtesy of www.dicardiology.com</figcaption></figure>

Two decades later, surgical procedures have also faced a paradigm shift. The introduction of laparoscopic cholecystectomy (the surgical removal of the gallbladder) [BA11] encountered immediate acceptance by patients and surgeons. Compared to open cholecystectomy where a surgical incision of approximately 10 to 15 cm is typically made below the edge of the right ribcage, laparoscopic cholecystectomy only requires 4 small incisions (named operating ports) to let tubes of diameter from 5 to 10mm be inserted in the abdominal cavity. These tubes are then used to place instruments and a video camera in the patient belly. The surgeon monitors the procedure through a TV screen and acts by manipulating the instruments through the operating ports 1.2. This technique has retrospectively been proved to be safe and effective [BM94].

As a consequence and with the widespread of new technologies such as micro-instruments or imaging modalities, modern surgical procedures tend to be as less invasive as possible. Although this kind of surgery has tremendous advantages, offering quicker recovery, less pain, and greater safety, it presents new challenges for surgeons as they do not touch directly the patient but operate through specific (sometimes robotized) devices and have indirect visual feedback through video monitors. Because of this new complexity in medical procedures, medical institutions demand
<figure><img src="image_10.png" /><figcaption>**Fig. 1.2:** Illustrations of laparoscopic surgery procedure: 4 small incisions are required to insert surgical instruments and camera (left). The surgeons then monitor the procedure on a video screen while manipulating the long thin surgical instruments (right). Left image appears courtesy of Sunshine Coast Advanced Laparoscopic Surgery, right image appears courtesy of www.chronogram.com.</figcaption></figure>

new ways of training the interns before they face real patients. In that computer simulation is an ideal candidate as it provides a safe, reproducible, configurable, assessable test-bed and some simulators (classified as *Level-D simulators* by national aviation authorities) have already proven to be very effective for the training of pilots<fnref n="2" />.

## 1.2 Simulation Overview

When applied to surgical procedures, computer-based training systems rely on a combination of hardware and software. Hardware usually include a computational system with visual system (from video monitor to stereographic system) and an acquisition device meant to mimic surgical devices. A popular acquisition device is a robotic arm that will track the user motion and will be able to provide force feedback when touching objects in the simulation. The software part is obviously a software that allow to integrate all the hardware devices and runs the virtual simulation. An example of such system is provided in Figure 1.3 which shows a training system for neuro-surgical procedures.

<footnote n="2">²https://en.wikipedia.org/wiki/Full_flight_simulator</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 1.3:** Example of a neuro-surgical training system by NeuroTouch©. The system includes robotic arm equipped with stylus (Geomagic© Touch) to acquire user gesture and to render forces in case of contact and a stereo-vision display to mimic the use of binoculars.</figcaption></figure>

A computer software for bio-mechanical simulation can be summarized in a loop (called **Simulation Loop**) of 4 repeated steps 1.4: *Deformation Computation, Collision Detection, Contact Resolution* and *Rendering*. The *Deformation Computation* step aims at solving bio-mechanical equations in order to find out how the soft material has moved given specific forces applied to it. The *Collision Detection* step finds colliding (*touching*) or inter-penetrating (*object inside another*) objects involved in the simulation. Colliding objects can be virtual surgical tools manipulated by the user, therefore this step will acquire the current position of the tools. The *Contact Resolution* step will take care of these objects and will correct their position in order to provide a solution which is mechanically consistent and realistic. Forces will also be computed that can be sent to any force-feedback device. Finally the *Rendering* step is performed to provide a 3D / 2D / stereo visualization of the current state of the simulation. The *Deformation Computation* and the *Contact Resolution* steps rely on mechanical equations while the *Collision Detection* step relies on geometric algorithms. The *Rendering* step relies on Computer Graphics techniques and will not be addressed in this manuscript as I do not have proposed any original contributions in this field.
<figure><img src="image_12.png" /><figcaption>**Fig. 1.4:** A simplified overview of the simulation loop for bio-mechanical applications. These 4 steps are repeated to produce feedbacks (visual, haptic, audio...) to the user.</figcaption></figure>

This simulation loop complies with several technical constraints (stability, robustness...) but the most important constraint is induced by the user acting on the simulation. It implies that first, very few hypotheses on the motion / gesture of the user can be made and second, the simulation should react *fast enough* to look realistic and comfortable for the user. Thus, acceptable visual feedback rates are between 30 Hz and 60 Hz while haptic feedback are usually between 300 Hz and several kHz [FU06]. This means that the simulation software should output a new image every [16-33] milliseconds and a new force every [1-3] milliseconds which imposes high-performance computation on the simulation or to choose an adequate trade-off between accuracy and computation speed.

The core of any computer-training system for surgery is the computation of the deformation and the mechanics of deformable solids is quite a large topic [Dog13]. Computational models for deformable solids are a very active research field in the mechanical engineering and computer graphics communities. The concerns in computer graphics community focus on controlling the computational footprint and the ease of implementation and versatility [Nea+06]. Deformable solids are governed by the equation of the dynamics (Newton's second law of motion):
$$ M\dot{v} = f^{\text{ext}} - f^{\text{int}}(q, \dot{q}) + H^{\text{T}}\lambda \quad (1.1) $$

where $q \in \mathbb{R}^n$ is the vector of generalized degrees of freedom of the solid, $q$ represents the velocity, $M$ is a mass (intertia) matrix, $f^{\text{ext}}$ and $f^{\text{int}}$ are respectively the vectors of the external and internal forces. $H^{\text{T}}\lambda$ gathers the constraint forces contributions, this later term will be detailed later. Forces $f^{\text{int}}(q, \dot{q})$ are expressed as a function of the position of the solid and its velocity. Solving this equation analytically can be achieved on limited, special cases. When dealing with solids with arbitrary shapes and position-dependent forces, this equation requires a numerical approach to be solved. Thus the time-line is discretized as a set of time intervals for which the deformation will be computed. The temporal discretization is chosen according to the time constraints imposed by visual and/or haptic feedbacks. A fixed period (labeled $h$) is considered between time-steps whose value is usually between 1 milliseconds to 50 milliseconds. Integration of equation (1.1) over a time interval $[t, t + h]$ leads to:

$$ \begin{gathered} \int_{t}^{t+h} M\dot{v}dt = \int_{t}^{t+h} (f^{\text{ext}} - f^{\text{int}}(q, \dot{q}) + H^{\text{T}}\lambda)dt \quad (1.2) \\ M(v_{t+h} - v_t) = h(f_{t+h}^{\text{ext}} - f^{\text{int}}(q_{t+h}, \dot{q}_{t+h}) + H^{\text{T}}\lambda_{t+h}) \end{gathered} $$

Subscript $t$ and $t + h$ have been added to explicitly indicate when equation variables are considered. We use implicit (or backward) Euler method that states for any ordinary differential equation $\frac{dy}{dt} = f(t, y)$, the sequence of approximations are built using: $y_{k+1} = y_k + hf(t_{k+1}, y_{k+1})$. Then we can write that: $q_{t+h} = q_t + h\mathbf{v}_{t+h}$. Since internal forces are non-linear functions of the degrees of freedom and their velocities, we apply a Taylor series expansion and make the first order approximation (a single linearization by time-step):

$$ f^{\text{int}}(q_{t+h}, \dot{q}_{t+h}) \approx f^{\text{int}}(q_t, \dot{q}_t) + \frac{\partial f^{\text{int}}}{\partial q} dq + \frac{\partial f^{\text{int}}}{\partial \dot{q}} d\dot{q} \quad (1.3) $$

Combining equations 1.2 and 1.3 and posing $dq$ as $q_{t+h} - q_t$ (which also equals to $h\mathbf{v}_{t+h}$), we obtain the following equation:

$$ \left( \frac{M}{h^2} + \frac{\partial f^{\text{int}}}{\partial \dot{q}h} + \frac{\partial f^{\text{int}}}{\partial q} \right) dq = f^{\text{ext}} - f^{\text{int}}(q_t, \dot{q}_t) + \left( \frac{M}{h} - \frac{\partial f^{\text{int}}}{\partial \dot{q}} \right) \dot{q}_t + H^{\text{T}}\lambda \quad (1.4) $$
<figure><img src="image_13.png" /><figcaption>**Fig. 1.5:** Examples of manufactured surgical devices than can be modeled as one-dimensional solids: Surgical Thread (®AD Surgical), Endoscope (®International Light Technologies), Catheter (®Paralogic), Detachable coil (®MedTronic)</figcaption></figure>

which is a linear system of equations $A dq = b$ where terms $A$ and $b$ are estimated at each time-step and are specific to the mechanical model of the solid and the constraints we consider. Solving this linear system enables to compute the motion and the deformation of the solid over the time of simulation. We will detail in the following section how this equation could be applied to simulate a one-dimensional deformable object.

## 1.3 1D Mechanical Model

We choose to detail in this manuscript our contributions about one-dimensional deformable models. From an *application* standpoint, one-dimensional model can be used to simulate a large variety of objects in surgical training systems (see figure 1.5 for some devices): for instance surgical thread (for suturing tasks), minimally invasive or endovascular devices such as endoscopes, catheters but also anatomical structures such as blood vessels, ligaments or fiber muscles. And from *modelling* and *simulation* standpoints, it raises some of the challenges related to characterization and computational optimization.

In the clinical context of interventional radiology, detachable coil embolization is a minimally invasive procedure to treat aneurysm<fnref n="3" /> 1.6 where the vascular network is used as a medium to reach the disease vessel. The interventional radiologist starts by inserting a catheter (a long, thin and flexible tube) into the femoral artery. This catheter is then manipulated through the arterial system until the aneurysm location is reached. When complex branchings occur, a guidewire (a long, thin wire, stiffer than a catheter) can be inserted inside the catheter to help reaching the desired artery. Once in position, the physician places one or more small coils through the catheter

<footnote n="3">³An aneurysm often consists in an abnormal bulge that appears close to branching arteries</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 1.6:** Example of a coiling procedure in an aneurysm (blood vessel bulge): a thin wire (coil) is inserted in the aneurysm reducing the blood flow thus the blood pressure and allowing the blood to clot.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 1.7:** Examples of the diversity of manufactured detachable coils used in interventional radiology. Length, diameter and shapes (helical, 3D, 360 degrees) may differ and are used given the pathology encountered.</figcaption></figure>

into the aneurysm. The body responds by forming a blood clot around the coil, thus blocking off the aneurysm and considerably reducing the risk of rupture. Although coil embolization is less invasive than open surgery, and allows treatment of cerebral aneurysms that previously were considered inoperable, such procedures are very difficult to perform and require an important experience from the interventional radiologist. This is particularly true when treating brain aneurysms, where any mistake occurring during the coil deployment or when reaching the aneurysm could lead to a brain hemorrhage. Yet, even in the case of a successfully performed procedure, the choice of the coil plays a key role in the long term success of the procedure. Most of them have a core made of platinum, and are sometimes coated with another material or a biologically active agent and their principle characteristics are the shape, the length, the diameter (see figure 1.7). All types are made of soft platinum wire of less than a millimeter diameter and therefore are very soft. The softness of the platinum allows the coil to conform to the often irregular shape of the aneurysm, while the diameter, length and shape of the coil are chosen based on the shape and volume of the aneurysm, as well as the size of the neck of the aneurysm. In most cases, several coils are required to completely fill the aneurysm and maximize the chances to clot (see figure 1.8).
<figure><img src="image_16.png" /><figcaption>Fig. 1.8: Aneurysm coiling requires several coils to fill the aneurysm in order to reduce the blood flow and thus prevents the aneurysm from rupturing. Image courtesy of Stryker Neurovascular.</figcaption></figure>

### 1.3.1 Contributions

Given this clinical context, we proposed a model for simulating very thin and flexible devices such as coils. An original approach was also introduced for describing the complex rest shape of coils. Co-axial insertions of endovascular devices (coil inside a micro-catheter for instance) are achieved using a *composite model*, which is based on a geometric combination of the characteristic rest shapes of each model. This composite model allows to use a single model to simulate several instruments together and then saves the computational burden of having several interacting objects to simulate.

The proposed model uses a series of serially linked beam elements, similarly as proposed by Duriez et al. [Dur+06] for simulating catheters and guidewires. A beam element, as introduced by Przemieniecki [Prz68] allows to model mechanically a piece of curve that can handle bending, stretching and twisting. Both extremities of the beam are called *nodes* and are described by six degrees of freedom, three of which correspond to the spatial position, and three to the angular position of the node in a global reference frame (see figure 1.9a). An elementary stiffness matrix $K_e$ is a 12 × 12 symmetric matrix that relates spatial and angular positions of each end of a beam element to the forces and torques applied to it. The notable parameters of the matrix $K_e$ are the mechanical properties of the beam (E the Young's modulus and $\nu$ the Poisson's ratio) and the geometric properties of the beam (length, cross-section area and cross-section moments of inertia). The reader may refer to [Dur+06] or [Prz68] for more details about the construction of $K_e$ given these parameters.

In the reference frame, only deformations (bending, torsion, elongation) are measured. A transformation matrix $\Lambda$ is then defined to change the local frame of
reference to a global coordinate system. This leads to the following relationship between $\overline{\mathbf{K}}_e$ in local coordinates and $\mathbf{K}_e$ in a global frame:

$$ \mathbf{K}_e = \Lambda^T \overline{\mathbf{K}}_e \Lambda \quad (1.5) $$

As the beam deforms, only $\Lambda$ changes and needs to be recomputed, while $\overline{\mathbf{K}}_e$ remains constant, as long as the deformation of the beam in its local frame remains small. To model a complex structure such as a coil, which undergoes large displacements, we need to discretize it as a series of beam elements to ensure that each beam deformation will remain small (see figure 1.9). For the entire structure describing a catheter or guidewire, the global stiffness matrix $\mathbf{K}$ is recomputed, at each time step, by summing the contributions of each beam element, through its elementary stiffness matrix $\mathbf{K}_e$. Assuming lumped masses at the nodes, the mass matrix $\mathbf{M}$ is a diagonal matrix. A damping matrix $\mathbf{D}$ is also introduced as a linear combination of the stiffness and mass matrices $\mathbf{D} = \alpha \mathbf{M} + \beta \mathbf{K}$, known as Raleigh damping. This leads to express the internal forces from equation 1.2 as:

$$ \mathbf{f}^{\text{int}}(\mathbf{q}, \dot{\mathbf{q}}) = \mathbf{D}\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} \quad (1.6) $$

However, we introduce several modifications to take into account for the particular nature of interventional devices. Boundary conditions are specified by defining a particular translation or rotation for the first node of the model to represent user control of the device (the coil is manipulated by pushing and twisting a wire). Even if $\overline{\mathbf{K}}_e$ is a linear approximation of Hooke's law, our model can handle large displacements: given an initial configuration of a beam, a current deformed one and a reference frame of the first extremity, the beam theory defines the way the local frame of the second extremity will be updated. Since the first node of the model is constrained, the first beam equation will be used to update the local frame for the second node, thus allowing the second beam to be computed in a reference frame where no rigid transformation occurs. By repeating this process through the whole structure, we can compute $\Lambda$ for each beam element and therefore determine $\mathbf{K}_e$. This method is closer to the co-rotational approach [Fel00] than the incremental approach proposed in [Dur+06] and permits to model the geometric non-linearities that occur during the deformation of the coil.

Devices such as catheters, guide-wires and coils are not subject to elongation but mainly bending and twisting. They are also characterized by their rest shape, which plays a very important role in the delivery of the therapy. We model these complex shapes by changing the local frame at each node of the model, as illustrated in figure
<figure><img src="image_17.png" /><figcaption>Fig. 1.9: Modeling of a wire using linked beam elements. (a) Given an initial rest shape and forces applied on DOFs, internal forces are computed using a local reference frame to handle large deformations. Starting from the first beam, reference frames are updated to suppress rigid transformation from the expression of the beam deformation. (b) Modification of the initial shape by moving frames or positions allows the generation of various deformed shapes. (c) A composite model (co-axial combination of two wires) is obtained by geometrically mixing different rest shapes on one geometric model.</figcaption></figure>
1.9b. Obviously, an important feature is to dynamically modify these frames at run-time to describe more complex behaviors. For instance, to simulate a coil being deployed from a micro-catheter, we avoid computing the complex interactions the take place between the coil and the micro-catheter. Instead we propose to simulate this behavior by using a composite model. The Halpin-Tsai [HK76] equations provide a framework to update mechanical properties of a reinforced fiber. This approach was applied to the simulation of catheters and guide-wires by Lenoir et. al. [Len+05]. We propose an extension of this framework to take into account two or more local frames and combine them using weighted sum (see figure 1.9c):

$$q_{cmp} = \frac{q_1 E_1}{E_1 + E_2} + \frac{q_2 E_2}{E_1 + E_2}$$

where each $q_i$ is the quaternion representing the transformation from local frame global frame and $E_i$ represents the beam Young modulus of the considered rest shape. As a consequence, such weighted sum will give more influence on the final composite shape to stiff materials.

From an optimization standpoint, it is worth mentioning that since all beam elements are serially connected, the resulting stiffness matrix $K$ is a tridiagonal matrix with a band size of 12 (since each $K_e$ is a 12 × 12 matrix). Since the mass and damping matrices are also diagonal, we solve the linear system using the algorithm proposed by Kumar et al. [KP93]. The solution can be obtained in $O(n)$ operations instead of $O(n^3)$. This allows to drastically reduce the computational footprint of the system resolution.

### 1.3.2 Coil Deployment in Collision-Free Environment

In order to validate our one-dimensional model, we want to compare our simulation against real data. However, exact mechanical properties of coils are difficult to obtain since they are not shared by device manufacturers. Among the parameters required to simulate our coil model, some of them are usually provided, such as coil length, diameter of the cross section and possibly some hints about the rest shape (circle diameter for helical-shaped coil for instance). Other parameters such as mass, Young's modulus and Poisson's ratio can be found in mechanical engineering handbooks. Coils are made of platinum and other materials such as titanium, giving us an initial guess for what Young modulus and Poisson ratio values to use. Eventually, the rest shape is the decisive parameter we need to identify because this feature is very important for embolization coils.
<figure><img src="image_18.png" /><figcaption>**Fig. 1.10:** Experimental setup: (a) the coil is deployed in a contact-free environment. (b) Volumetric data is obtained by 3DXA and a marching cube is performed to get a mesh of the real coil. (c) Central line of the coil is extracted and a continuous formulation is built using splines.</figcaption></figure>

In order to estimate these missing parameters we build an experimental setup that consists of a box filled with water and a transversal fixed guide, defining the path for the catheter. The catheter is first introduced, followed by the helical coil as in the real procedure. The three-dimensional shape of the coil at different stages of the deployment is obtained from 3DXA (3D X-ray angiography) images using a marching cubes algorithm. The different steps of the experiment are illustrated in figure 1.10. In these experiments, the coil is only subject to a gravity force, and contacts with the box were avoided.

Once the real coil positions have been tracked and the shape segmented for three different stages of the coil deployment, an algorithm using a B-spline formulation was applied to the real coil data to obtain a continuous description of the shape. Given this dataset, we run an optimization (brute-force gradient descent) algorithm to estimate the rest shape of the coil as well as mechanical parameters. A validation of the coil simulations is reported in the next table for three different stages of the deployment. Then, simulations of the coil were performed using the same boundary conditions as measured during the actual coil deployment. In order to compare our simulation with the actual data, the B-spline curve was discretized using the same number of points as the number of nodes in the simulation. The error metric used to validate the coil simulation is the relative energy norm error [ZZ87]. The metric is defined as the ratio of the displacement between the simulated and the real coil and the real displacement. The error is computed for each node and for the three different steps of the deployment. The results are given in table 1.1. The relative energy norm error is small and exhibits that our simulation is close from real data. The absolute mean distance between simulated and real positions after deployment is equal to 1mm for each stage, explaining why the error increases if the coil is deployed on a small length. An illustration comparing the simulated coil and the real coil is proposed in the figure 1.11 showing the accuracy of the model.
<table><caption>**Tab. 1.1:** Error measurement. The average and the standard deviation (SD) of the relative energy norm error are given in percentage in the first column. The second column contains the mean displacement of the coil for the three different steps.</caption><thead><tr><th rowspan="2"></th><th>Relative Energy norm error</th><th>Mean total displacement</th></tr><tr><th>Average % (SD)</th><th>(mm)</th></tr></thead><tbody><tr><td>Stage 1</td><td>9.80% (4.05)</td><td>18.2</td></tr><tr><td>Stage 2</td><td>6.19 % (5.05)</td><td>24.99</td></tr><tr><td>Stage 3</td><td>4.17 % (1.42)</td><td>39.11</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>**Fig. 1.11:** Visual comparison of our coil simulation (in red) with the reconstructed coil model from our experiments (in yellow) at different stages of the deployment. Quantitative errors are available in table 1.1.</figcaption></figure>
1.4 Asynchronous Simulation

Given mechanical models for deformable solids, simulation of these models can be
included in a single virtual reality environment. In the context of medical simulations,
a common scenario considers a set of soft anatomical structures (soft tissues, organs,
blood vessels...) which are manipulated or cut by surgical devices. Compared to
these anatomical structures, the surgical devices appeared to be much stiffer and can
sometimes be modeled as rigid objects (for micro-surgical blades for instance). But
some devices like catheters, guide-wires or surgical threads, while being composed
of stiff materials, are thin enough to be bended and might be subject to high-speed
deformations or vibrations. Combining *soft* and *stiff* objects in a simulation raises a
computational issue:

* *stiff* objects are usually simulated with dedicated models (such as described in section 1.3) and are composed of several degrees of freedom (DOFs). They can be therefore simulated interactively with a small time-step to capture the deformations occurring at high-frequency.

* *soft* objects are built with generic FEM models with a large number of DOFs mainly because of the number of geometric primitives used to describe their arbitrary shapes. They therefore require a large time-step to be simulated interactively.

Having the two objects interacting in the same simulation, like the *stiff* needle puncturing the *soft* liver, requires the objects to share a single time-line and cannot be simulated in parallel. Such an approach can be acceptable in the case where the simulation is not required to be run in a real-time: in this case a small simulation time step is chosen to account for the high frequencies in the needle response. Thereby, it's assumed that it will run much more slowly than the real clock time, as in each small time-step, the costly dynamics of the bulky soft objects has to be recomputed.

In the context of interactive simulations, such a solution cannot be considered. The
strategy we have proposed is a multi-rate approach where dynamics of each object
in the scene is integrated in real-time with its own time-step. This approach is
especially suited for simulations with haptic rendering. In this case, the dynamics of
the stiff object (e.g. the flexible needle) which does not require heavy computations
can be simulated in the haptic loop running usually at 1000Hz, while the dynamics
of the soft highly-deformable body (e.g. the liver) is simulated on a much lower
(typically synchronized with the visual thread running about 30 or 60Hz), since it
requires more computationally expensive calculations. This approach was already studied for example in our previous works such as [DGC05a] where the objects were simulated at different frequencies, however, at the cost of simplification in modeling the interactions. The contributions we present in this sections are based on [PDC11] where a method allowing for multi-rate interactive simulation of objects simulated at two different frequencies has been proposed. However, the method was suitable only for a quasi-static scenario: it required that compliance matrices needed to simulate the interactions between the objects do not depend on time step.

Our framework relies on simulation loops running at different frequencies. A simple explanation of our strategy can summarized as follows: time-consuming operations will be performed in the low-rate simulation loop while fast computations will be performed in the high-rate loop:

* *low-rate loop:* unconstrained motions, proximity queries,

* *high-rate loop:* unconstrained motion (stiff objects), constraint solving (violation evaluation, compliance computation, resolution, corrective motion), haptic-force computation

The details and more specifically the equations are detailed in [Der+13] but we want to emphasize on the following features. The models of both stiff and soft objects are linked by the constraint based expression of their mechanical interactions. Among the interactions between the two deformable bodies, we need to detect the collision and provide an adequate response (without inter-penetration). The collision detection (or more generally the geometrical determination of the contact points) is often time consuming and is not compatible with the high rate loop. In our approach, we use an algorithm based on proximity queries between meshes that is computed in the low rate loop. It supposes that we have access to a geometrical position of the stiff object in this low rate loop. The algorithm places contact constraints where the local minimal distance between the meshes is small enough to have a potential collision between current time *t* and *t* + *h*. Other type of constraints, for instance bilateral constraints (attachment, sliding...) can also be set. After computing the queries, a matrix *H* (see equation 1.1) containing the directions of unilateral and bilateral constraints for each object is assembled. To solve the reaction forces related to constraints, we first to need evaluate the constraint violation when applying an unconstrained dynamic motion to both deformable objects (i.e. solving equation 1.1 with $\lambda$ set to 0). This unconstrained motion is also called *free motion* and allows to obtain the constraint violation $\delta^{free}$ (which is computed at low and high rates). Solving a MCP (mixed complementarity problem) allows to find the final constrained positions.
<figure><img src="image_20.png" /><figcaption>Fig. 1.12: Asynchronous simulation: details on the steps performed by the two loops that run at different frequencies. The central blocks represent shared buffers that are needed to compute mechanical interactions.</figcaption></figure>

To validate this method, we created a simulation scene where a cylinder is placed on a horizontal plane. The cylinder is modeled by tetrahedral corotational FEM and is simulated at 50 Hz. A wire is attached to the upper part of the cylinder and the other endpoint of the thread is driven by a haptic device as shown in figure 1.13. The wire is modeled with 10 serially-linked beam elements based on Tymoshenko formulation and its dynamics is computed at 1000 Hz. The collisions between the plane, cylinder and the thread are being detected and resolved during the simulation. During the interaction, the simulation remains stable for any admissible configuration: the position of the cylinder on the plane can be changed by pulling the wire, or the cylinder can be even pulled out of the plane, so it's hanging freely in the space been hung on the wire. At the same time, a detailed and realistic haptic force-feedback is delivered to the user.

## 1.5 Applications: Medical Simulators

The methodology employed through this chapter has been applied to other scientific problems related to mechanical modeling or simulation interactivity. In order to illustrate the potential and the relevance of the contributions, we have built several proof-of-concept medical applications that rely on interactive simulations. Some have been a necessary step to open new collaborations or to build more ambitious projects and other have been transfered to start-up companies to be commercialized. In this section, we will focus on two specific applications: first, a simulator to train
<figure><img src="image_21.png" /><figcaption>**Fig. 1.13:** Simulation of a deformable cylinder (simulated at 50 Hz) and a thread (simulated at 1 kHz) interacting with a solid plane.</figcaption></figure>

cataract surgery procedure and second, an application that help radiologists to plan coiling interventions.

### 1.5.1 Training System

The development of a training system for ophthalmology has been a good opportunity to showcase the various contributions in interactive simulation while strengthening a collaboration with the CHRU of Lille that has proven to be very fruitful.

When designing the simulator, our objectives were:

1. be able to train on the three most important steps of the procedure: capsulorhexis, phacoemulsification, and intra-ocular lens injection which has not, to our knowledge, been treated previously
2. develop advanced models based on biomechanics that can be easily parametrized using material properties and that can be computed quickly

3. have a working environment and devices as close as possible to real operating rooms.

## System overview

The environment consists in a microscope where stereo glasses have been fixed in order to render the virtual scene in 3D, a full-body mannequin where replaceable silicon eyes that support cataract surgery may be inserted. Monitors are also available in the room to provide additional informations (like vital constants) or feedbacks. This serves a pedagogic purpose because it increases the feeling of immersion while allowing to train and memorize the real procedure with the same devices (microscope, pedals...) and the same surgical tools design using rapid prototyping. The tools are tracked in real-time and their 3D position is used to feed our real-time simulation. Figure 1.14 illustrates the various components that are part of the training system and how they are connected.

<figure><img src="image_22.png" /><figcaption>Fig. 1.14: Functional view and picture of the training system. The system relies on devices used in real procedures (microscope, pedals...), replicas of surgical tools that are tracked in real-time and on a simulation engine that renders a realistic and interactive simulation of cataract surgery.</figcaption></figure>
**Interaction Devices** Real surgical tools could have been used to interact with the simulation, however they suffer from two important drawbacks. First they are expensive and second they are metal-based which makes the tracking more difficult and lower the accuracy (diffuse reflection area may appear on the images). Therefore, the use of replicas of surgical tools based on polyamide (nylon) was made. They were first modeled using standard CAD software and then were built using fast prototyping. Figure 1.15 shows the comparison between some of the real surgical tools and the ones used in our simulation. Reflective markers may disturb the clinician, especially for grasping tasks, but the markers are carefully placed to limit disturbance.

<figure><img src="image_23.png" /><figcaption>Fig. 1.15: Examples of tools used in cataract surgery and the ones we built for interacting with the simulation. From left to right: the phaco chopper, the phacoemulsification hand-piece and the capsulorhexis forceps.</figcaption></figure>

**Robust 3D Tracking** Our setup includes six infra-red cameras that targets the area around the eye. Two cameras may be sufficient to provide a 3D tracking but this is not a robust set-up. Indeed, intensive manipulations in that area may completely or partially hide instruments preventing the real-time tracking to work. With six cameras, the operative area is well covered and during any step of the procedure at least two cameras are able to track the instruments. The setup is shown on figure 1.16. The redundancy involved with the six cameras make the calibration more time-consuming but once it has been performed, the tracking algorithm is able to provide 3D position of markers in real time. We are also using a feature proposed by *Tracking Tools*<fnref n="4" /> software which allows to track rigid bodies. Each rigid body is defined as a set of reflective markers and the distance between markers do not change over time and the tracking is then more simple. Three markers may be sufficient to define a 3D position but again for robustness issues, we introduce redundancy by placing four markers on the instrument. Using different signatures

<footnote n="4">⁴http://www.naturalpoint.com/optitrack/products/tracking-tools/</footnote>
(e.g. the way markers are placed on the tools), we can manage to track several tools simultaneously (see figure 1.17) which is a prerequisite for several steps of the cataract surgery procedure.

For a better balance of the computation power needed to run our simulation system, two computers are used. The first one, which does not require large computational resources, is connected to the cameras and runs the tracking software while the second one runs the simulation engine. The two computers shared tracking informations by using the VRPN library<fnref n="5" /> which is a lightweight and low-latency client-server protocol dedicated to virtual reality applications. The simulation engine is therefore on a single computer (with high-performance graphics cards) and can provide fast response to user inputs.

<figure><img src="image_24.png" /><figcaption>**Fig. 1.16:** Overview of our setup: mannequin and microscope where 6 infrared cameras are attached. These cameras allow to track in real-time the motion of the surgical tools.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 1.17:** Screenshots of the tracking of surgical tools. With the reflective markers, the tracking system is able to compute the 3D position of the tool in real-time and using different marker signatures, we can manage to detect several tools simultaneously without ambiguity.</figcaption></figure>

<footnote n="5">⁵http://www.cs.unc.edu/Research/vrpn/</footnote>
Simulation considerations

**Bio-mechanical modeling** The deformable lens is modeled using a co-rotational FEM formulation with linear elasticity. It is represented by a tetrahedral mesh that can be tessellated to produce small elements (see Figure ??). We chose an implicit time integration scheme in order to achieve large time-steps and provide a stable simulation even when subjected to challenging interactions. We use the approach detailed in [Cou+10], where the resulting system of equations is solved iteratively using a Conjugate Gradient solver implemented on the GPU.

<figure><img src="image_26.png" /><figcaption>Fig. 1.18: Tessellation algorithm that can dynamically refine tetrahedra (two images on the left). Two-levels tetrahedral topology for lens phacoemulsification (two pictures on the right).</figcaption></figure>

One important aspect of this method is that the mechanical matrix of the full system is not explicitly computed, but matrix-vector products are instead parallelized directly over the mesh elements. This allows to easily change the topology of this mesh, as is necessary in this application to implement the interactions with the phacoemulsifier explained below.

The approach used to model the complex deformations of the implant during its
insertion in the lens capsule is based on the use of triangular shell elements and
a co-rotational formulation [CCD10; Bou+09]. The combination of both leads to
an accurate, yet computationally efficient, shell finite element method featuring
both membrane and bending energies. In addition, the polynomial shape functions
employed to compute internal forces in our FEM formulation are also used to
process contacts and interactions with the curved triangles. The benefits of our
approach (fast computation and possible interactions with other objects) allow us
to successfully simulate the complex step of the insertion and deployment of the
intra-ocular implant.

**Interactions** Besides the biomechanical modeling of the anatomical structures of the eye and the mechanical modeling of the surgical devices, several specific interactions happened during the procedure:
<figure><img src="image_27.png" /><figcaption>**Fig. 1.19:** An actual intra-ocular implant (courtesy of Alcon) and the triangular mesh used in our simulations. Notice the higher density of elements in areas where large deformations will take place.</figcaption></figure>

**Aspiration:** The phacoemulsifier proposes two modes of interactions to the surgeon. The first one, when activated, pumps liquids from the eye into the instrument, creating a force pushing the lens toward the tip of the phacoemulsifier.

**Fragmentation:** The second one, a probe at the tip vibrate at about 40,000 cycles per second, creating ultrasonic waves that break up the lens into small pieces that are sucked into the instrument by the aspiration and removed from the eye.

**Tearing:** *Continuous Curvilinear Capsulorhexis*, a technique used to create a circular opening in the lens capsule. Once a small initial incision has been made, this technique relies essentially on the application of shear and stretch forces to create the opening. When performed correctly, the continuous curvilinear capsulorhexis creates a smooth opening, reducing the risk of tear when forces are applied to the capsule during surgery. The main risk with the continuous curvilinear capsulorhexis is to either create an opening that is too large or to create an edge that could promote a tear proceeding outwards.

The aspiration aspect is handled by defining a spherical volume around the type and finding parts of the lens within it. When the aspiration is activated, an artificial force is added to pull the enclosed particles toward the tip. When fragmentation is also activated, we remove tetrahedral elements from the geometrical model of the lens. The removed elements should be as fine as possible to provide an accurate visual feedback to the surgeon. However, the mechanical mesh need to be relatively coarse to maintain good performances. To solve this issue, we use a finer tessellated mesh for rendering and interactions with the instrument. This is done by splitting each tetrahedron in the simulation mesh into eight smaller tetrahedra (Fig. 1.18). Therefore finer elements are removed due to aspiration. And if the eight finer elements of a coarse tetrahedron are removed, the coarse tetrahedron is also removed in order to mechanically take into account the loss of matter in the
lens. This allows the use of a relatively coarse mechanical mesh while using a more
detailed visual mesh for the rendering (Fig. 1.18).

<figure><img src="image_28.png" /><figcaption>**Fig. 1.20:** Three screenshots of the phacoemulsification simulation. Local fragmentation of the lens create a cross-like pattern on the lens.</figcaption></figure>

In our approach, we propose to model the deformation and tearing of thin soft
tissue (such as membranes, capsules, etc.) by using a transversely isotropic elastic
model. The model is based on a co-rotational finite element method formulation,
using triangular elements on which a specific fiber directions can be defined (see
Figure 1.21). During deformation of the tissue, this fiber orientation, and the
principal strain direction on the element are used to determine the tear propagation.
Tearing or fracture within an element only occurs when a threshold is reached, which
is computed using an eigenvalue decomposition of the strain tensor in each element.
If the largest eigenvalue is above a given threshold (in our experiments, the value is
tuned according to the feedbacks provided by the surgeon), the element is identified
as "breakable" and the eigenvector associated with this eigenvalue corresponds to
the principal strain direction (see figure 1.21). Since tearing tends to propagate
from an already fractured location, we also account for the history of tear location
and direction in the overall computation of the current tear direction. Topological
changes on the FEM mesh are performed at each time step. Each triangle in the
neighborhood of the current tear location which strain is above our threshold will
be re-meshed. This local re-meshing subdivides a triangle into new triangles. There
is no particular requirement on the local re-meshing besides avoiding very small
triangles which lead to ill-conditioned systems.

Comments

The results presented in this section illustrate that the work required needed to
step from the research contributions to a working simulation prototype. These
contributions rely on advanced biomechanical models and intensive use of modern
graphics hardware to provide fast computation times. Moreover, surgical devices,
that are tracked in real-time thanks to infra-red cameras, have been reproduced in
<figure><img src="image_29.png" /><figcaption>**Fig. 1.21:** From left to right: tearing propagation throughout the tissue. Fibers (not visible on the meshes since they are defined per triangles) are oriented with a 45 degree angle, and forces are applied outward, on the left and right sides of the square mesh. Our fracture criteria is represented as the shaded area. One can see on the left image how the anisotropy affects the deformation of the tissue, and on the rightmost images how the fracture direction depends on the previous fracture direction, as well as fiber direction. The square mesh is composed of 1,500 triangles. Anisotropic deformation, fracture and remeshing are performed in real-time.</figcaption></figure>

order to faithfully mimic the operating room environment and the actual procedures.
Preliminary tests have been performed by experienced ophthalmologists to qualita-
tively look at the realism of the simulation and the faithfulness of the interactions
with the surgical tools.

As these feedbacks were positive a double-blind experiment was supposed to be
conducted among ophthalmologist interns to assess the effectiveness of simulation
training compared to current companionship training methods. The common ground
used to measure surgeons performance during cataract surgery is based on [Sal+07]
and is an objective performance rating tool derived from the Objective Structured
Assessment of Technical Skill (OSATS) [Mar+97]. OSATS has also been validated
for other surgical specialties. The work of Saleh et. al. on OSACSS (Objective
Structured Assessment of Cataract Surgical Skill) exhibits construct validity and is
therefore a valuable tool to assess the surgical skill of trainees [Sal+07]. OSACSS is
based on a set of 20 measurements or observations that can be programmed into the
simulation engine. For instance, observations include *Eye Positioned Centrally Within
Microscope View, Iris Protection or Capsulorhexis: Formation and Circular Completion*
and can be easily translated into computer instructions. Unfortunately this study has
been preempted before the end as interns willing to conduct the study left before
the end of the study.
<figure><img src="image_30.png" /><figcaption>**Fig. 1.22:** Simulation of capsulorhexis, a technique used to create a circular opening in the lens capsule. This technique relies essentially on the application of shear and stretch forces to propagate a fracture throughout the membrane. Top (a): the mesh used to support the computation of the deformation and topological changes. Bottom (b): final result as seen in the training system.</figcaption></figure>

## 1.5.2 Planning System

In the case of interventional radiology (see section 1.3 for the clinical context), the overall morbidity and mortality associated with endovascular repair is lower than for surgery (9.5% at one year, compared to 12.2%) the technical challenges associated with coil embolization remain important. It was however shown recently that the risk of complications with coil embolization of unruptured aneurysms decreases dramatically with physician experience [Sin+02]. This clearly emphasizes the need for simulation systems that can be used for physician training as well as for pre-operative planning. This section describes how a simulation tool can be used to select the proper characteristics of the coil based on the patient anatomy, but to also deploy interactively the virtual coil (s) to plan or rehearse the procedure. The interactive aspect in this planning plays a key role as the interventional radiologist constantly controls the motion of the coil and as such influences the final coil placement within the aneurysm.
<figure><img src="image_31.png" /><figcaption>Fig. 1.23: Reconstructed vasculature using 3DXA imaging modality. The 3D volume was segmented and reconstructed using triangles.</figcaption></figure>

Using the approach presented in section 1.3, the geometry of aneurysms and nearby vessels can be reconstructed from a volumetric 3DXA data set (figure 1.23). Then, an actual coil has been deployed inside the aneurysm of the digital model and its shape, during its deployment, has been reconstructed. We then simulated the deployment of a virtual coil in the reconstructed model of the aneurysm. The coil was represented using 100 beam elements. Our dynamic coil model, in combination with the fast solver and contact model introduced above, allowed an interactive and computationally efficient simulation of the coil deployment. In this simulation, a friction coefficient $\mu = 0.6$ was used, and due to the rigid nature of the phantom model, we did not integrate the compliance for the vessel walls or aneurysm during the contact response. While this may seem as a simplified problem, it is actually more complex to handle numerically. Yet, the simulated coil embolization we have obtained exhibits a similar dynamic behavior and final configuration as what can be observed during an actual procedure (see pictures 1.24). Computation times for the whole simulation including coil deformation, collision detection were of 50 ms in average, although they varied depending on the number of contacts.

Figure 1.25 is a close-up on the final stage of the coil deployment and illustrates that our whole approach prevents penetration of the coil inside the vasculature wall as well as self-penetrations of the coil which is mandatory to have a realistic simulation. The mechanical model also allows large deformations of the model making the coil possible to conform to the aneurysm shape.

Using the simulation for coil deployment enables to have additional information that can be useful for clinicians. Since the contacts between the surgical devices and the
<figure><img src="image_32.png" /><figcaption>**Fig. 1.24:** Examples of our simulation results: (two images on the left) real coil embolization; (two images on the right) our simulated coil embolization with 3D coils.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Fig. 1.25:** Examples of our simulation results: close-up on the coil deployment and wire-frame display to illustrate the deformation of the coil.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Fig. 1.26:** Examples of our simulation results: pressure of the coil on the aneurysm surface can be displayed using a thermal scale. Using a larger coil (right) induces more important pressure on the aneurysm wall which could lead to rupture.</figcaption></figure>

anatomy are handled with mechanical constraints, pressure on the surface can be computed. Figure 1.26 illustrates the kind of additional data that can be displayed for the radiologists: pressure was color-coded using a thermal scale (blue=low to red=high) to enable a better comprehension of the locations of highest pressure on the aneurysm wall during coil deployment in real-time. Contact forces were also depicted as yellow arrows whose length is proportional to the intensity of the forces. Choosing a larger coil implies a significant increase of pressure forces on the aneurysm wall.

Having the simulation software can be used by experienced or novice radiologists to test different strategies in a safe environment. For example, given a specific patient with an aneurysm, users can try different coils (helical, 3D...) or insertion strategies to select the one that seems to be optimal given the aneurysm features. Figure 1.27 is an example of this planning strategy: 3 different helical coils (4mm, 7mm and 10mm diameter) have been inserted in the aneurysm and the simulation illustrates the role of the diameter in the deployment: the smallest one cannot fit the aneurysm cavity and is floating inside the bulge making it very dangerous for the patient (can move into the parenting artery and thus into the brain), whereas the largest one induces increased pressure on the wall that can lead to rupture (and cause hemorrhage).

The placement of the first coil is an essential step to guarantee the success of the procedure: it should cover most of the surface of the aneurysm and should be firmly anchored on the aneurysm wall as it will served as a support for the following coils [Sat+98; Tam+02]. A common metric that is routinely used in the CER (Coil Embolization Ratio) or VER (Volume Embolization Ratio) that has been to be positively correlated with the long term stability of the embolization [Clo+00; YAG+05]. The
<figure><img src="image_35.png" /><figcaption>**Fig. 1.27:** Examples of our simulation results: using helical coils with different diameters: (left) 4mm helical coil which is too small for the aneurysm and the coil is floating inside the aneurysm and this can be dangerous for the patient, (middle) 7mm helical coil which is the optimal size for this aneurysm, (right) a 10mm helical coil which is too big and leads to extra pressure on the aneurysm wall.</figcaption></figure>

CER requires volumetric data of the aneurysm in order to be computed which cannot be done in every clinical set-up. For a fast metric that can be done with any imaging modality, we propose to derive the CER to the CESR (*Coil Embolization Surface Ratio*) because it combines the several advantages: close to the current clinical usage of radiologists, easy computation that can be done in the working imaging view during the procedures. Computation is a simple ratio between the surface (in pixels) occupied by the first coil and the surface (in pixels) of the aneurysm.

<table><caption>**Tab. 1.2:** CSER computation for two patients. *REF* corresponds to the CSER of the coil inserted by the clinician, *OK* represents the CSER of the simulation of the same coil, *NOK* represents a simulated helical coil. The relative errors to the *REF* are provided in parentheses.</caption><thead><tr><th></th><th>REF</th><th>OK</th><th>NOK</th></tr></thead><tbody><tr><td>Patient 1</td><td>81% (0%)</td><td>68% (16%)</td><td>55% (32%)</td></tr><tr><td>Patient 2</td><td>24% (0%)</td><td>17% (30%)</td><td>11% (54%)</td></tr></tbody></table>

Figure 1.28 illustrates the computation of CSER for two different patients and the results are reported in the table 1.2 where 2 patients have been treated with 3D coils. The first column represents the CESR from real road-map image (*REF*), the second column represents the CESR from simulation with the right 3D coil (*OK*), the third column represents the CESR from simulation with a wrong helical coil (same diameter) (*NOK*). Given this scenario, the following comment can be made: the CESR is able to distinguish between a right (3D) coil and a wrong (helical) coil. This difference appears clearly in the case of the bilobated aneurysm (patient 2). For patient 1, the rounded shape of the aneurysm justifies a lower relative error for the helical coil. These first results are interesting but the quite large relative errors for the right coil require a refinement of CESR as a quantitative score.
<figure><img src="image_36.png" /><figcaption>Fig. 1.28: Computation of the Coil Embolization Surface Ratio metric: it divides the surface of the coil by the surface of the aneurysm.</figcaption></figure>

### 1.5.3 Discussion

The previous application examples illustrate that our simulations perform qualitatively well compared to actual surgical or interventional procedures. However a quantitative assessment of the performance is difficult to conduct: many software components need to be evaluated (bio-mechanical models, interactions, environment modeling...), boundary conditions are badly-known and ground-truth data depends on imaging modalities that exhibit several limitations in terms of resolution, dynamic range and noise. In order to address these particular challenges, the following chapter will detail some contributions around the use of medical images to extract relevant information for the simulation.
Geometric Modeling for
Interactive Simulation

“Arithmétique ! algèbre ! géométrie ! trinité
grandiose ! triangle lumineux ! Celui qui ne vous
a pas connues est un insensé !”

— Lautréamont
Les Chants de Maldoror

Obtaining digital representations from medical images is a two-step process that combines first the segmentation of the image and second the reconstruction of the digital model. Segmentation consists in splitting a raw image into one or several structures of interest and isolating non-relevant information (i.e. the background). From 2D, resp. 3D, arrays of pixels, resp. voxels, segmentation consists in classifying each pixel / voxel according to the structure it belongs to. While manual segmentation is possible, it is a time-consuming task and subject to operator variability. Automatic segmentation is therefore encouraged and mainly relies on signal processing techniques. Reconstruction consists in building a digital representation from the processed image according to geometric modeling techniques. While this process is very common in the Computer Vision community, several constraints inherent to medicine makes this problem much more challenging: some are related to imaging, some are related to the patient itself.

Whatever the imaging modality and despite impressive progress in imaging, segmentation faces three image-dependent problems. First, images are noisy and may alter pixels intensity making the classification more difficult. Second, images may exhibit intensity nonuniformity where the intensity of a single structure varies over the extent of the image. Third, images have relatively low definition (often around 256/512 pixels for each dimensions), which induces *partial volume averaging* meaning that an individual pixel contains a blend of several structures and thus the intensity of a pixel may not be related to the intensity of pixels from a given structure. For dynamic imaging such as dynamic MRI, the acquisition rate is not fast enough to capture very fast motion leading to blurry data. And at least, some structures are not visible under specific modalities (for instance ligaments in MRI).
There is a huge inter-patient variability for anatomic structures which means that for a given structure an average representation or template cannot be used as an *a-priori* to ease the segmentation step. Moreover even manual segmentations performed by experts present uncertainties: *inter-observer variability*, variability where several experts perform the segmentation, as well as *intra-observer variability*, variability when a single export perform the same segmentation multiple times, has been recorded.

Thus, from a formal standpoint, the segmentation of medical images can be seen as under-determined problem where there are multiple classifications from a single image. The main challenge when developing segmentation algorithms is to propose methods that try to overcome this image limitations while producing acceptable results for experts. This is a very active field in research and the medical imaging community have proposed many different approaches. First approaches used image processing technique such as thresholding [SV16], region growing [HK98], edge tracing [HL00], while more recent works have used statistical analysis [BHC93], clustering [Pha01], deformable models / active contours [KWT88], neural networks [Bar+94] or atlases [Fis+02].

Segmentation has been extensively studied in the medical imaging community but remains on open-problem. Reconstruction has however received less attention and this can be explained by the first applications of medical imaging which were mostly used as a tool for diagnoses. Therefore the smoothness, the compactness, details outside the area of interest were not primary criterion for the reconstruction. As a result first models were only voxel soups who exhibit jagged surfaces 2.3. This was enough to give some hints or overview about the targeted area but the overall shape lacked realism. Besides these visual considerations, the smoothness of shapes is an important feature to have for mechanical simulations since irregular or jagged surfaces could lead to artificial sticking of objects instead of sliding. For instance, the figure 2.1 illustrates the influence the importance of friction in the context of catheter navigation: since the main behaviour of this navigation relies on the catheter sliding along the vasculature wall and the use the support of the vessel to take a specific vessel in case of branchings, the change of friction coefficient or artificial friction introduced can lead to very different simulation. To illustrate the importance of static friction phenomenon, Alderliesten et al. [AKN07] have conducted a thorough quantitative experimental validation of a quasi-static model of guidewires: 3D comparison through 3DRA is performed with several deployment scenarios.
<figure><img src="image_37.png" /><figcaption>**Fig. 2.1:** 200msec recorded sequences of catheter deployments in a small cylindrical vessels. Catheters share the same mechanical parameters and only the friction coefficient between the vessel wall and the tool is modified. This illustrates the importance of adequate friction estimation and smooth geometry.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Fig. 2.2:** Recorded sequence of a coil deployed at low and constant speed where image acquisition is performed at 15 fps through 3D Rotational Angiography. The sequence shows the effect of friction and then sliding as the coil completely changes its configuration between frames 3 and 4 (67 msecs between each frame).</figcaption></figure>

Dynamic friction is also to be considered because it can significantly change the trajectory followed by the devices. Even if the friction between the arterial wall and the tools may be negligible in certain circumstances (for instance slow motion, simple path to reach the targeted area), its effects are significant in tortuous vasculature or aneurysms. The figure 2.2 illustrates that friction implies high frequency motion for the deployment of a coil (thin coated platinum thread) in an aneurysm. Therefore, accurate models for the tools and for the interactions between tools and the blood vessel are mandatory.

A second feature that is less important for the accuracy of the simulation but is essential for optimization is the compactness of the virtual 3D model, i.e the number of primitives used to describe the virtual model: not only a more compact model implies a smaller memory footprint but since collision detection step usually has to inspect all the primitives, this is a time-consuming procedure. Without any optimization the cost of the collision step, is quadratic with respect to the number of objects in the virtual environment and quadratic to the number of their primitives. Obviously, optimization strategies exist [Zac01] but even with these optimization, it is one of the most costly step of the simulation loop. Indeed, dealing with deformable objects invalidates optimization based on geometry pre-processing because the geometry may change at every time-step. Having a compact model may reduce the computational burden of the collision detection. Dedicated ways of representing shapes such as distance fields or implicit surfaces may also be investigated.
<figure><img src="image_39.png" /><figcaption>Fig. 2.3: Reconstructed liver and its vasculature: segmentation and reconstruction have been performed without smoothing. The limited resolution of the imaging device leads to virtual 3D models that exhibit jagged surfaces which are not very organic.</figcaption></figure>

The contributions detailed in this chapter will illustrate two examples of geometric modeling: starting from medical images, specific and *simulation-aware* strategies will be used to obtain 3D geometric models that are optimized for mechanical simulation applications. But before digging into these two contributions, we will elaborate on implicit surfaces and illustrate that it is a representation that is very convenient for mechanical simulation.

## 2.1 Smooth surface using implicit surfaces

An implicit surface is based on a function $f$ and a scalar constant $c$, called *isovalue*. The function $f$ associates to any point $P(x, y, z)$ of the 3D space, a scalar value $f(P)$:

$$ f: \mathbb{R}^3 \rightarrow \mathbb{R} \\ \mathbf{P} \mapsto f(\mathbf{P}) \qquad (2.1) $$

The resulting space obtained by applying $f$ to $\mathbb{R}^3$ gives a scalar field called *implicit* or *potential* field. The equation $f(P) = c$ describes an isosurface, separating the values greater to $c$ of the implicit field from the values lower to $c$. Usually a $g$ function
is introduced which is $g(P) = f(P) - c$. Such surface is usually discretized and parametrized to be displayed.

An interesting property is that the gradient of $f$ defined as:

$$ \nabla f = \frac{\partial f}{\partial \boldsymbol{P}} \qquad (2.2) $$

gives the normal to the isosurface that passes through $f(p)$ and $-\nabla g$ points to the direction that reaches the faster the higher values of the field. Such property is useful for contacts resolution with implicit surface.

The choice of function $f$ is a complex step there is no convenient way to extract the function from an arbitrary shape. In the case of tubular structures, such as blood vessels, Bloomenthal [Blo95] introduced convolution surfaces which are skeleton-based implicit surfaces. The surfaces are described with few geometrical primitives (called skeleton $s$) and a convolution filter $h$ which is a function from $\mathbb{R}^3 \to \mathbb{R}^+$ and has a finite support or fast decay to 0. The resulting surface is defined as $\forall P \in \mathbb{R}^3, f(P) = h(P) \otimes s = \text{iso}$, where $\text{iso}$ is an isosurface value and $f(P)$ the potential at $P$. The convolution surface then produced is smooth and very appropriate to design organic 3D objects. The figure 2.4 represents a simple vessel-like network built with cylindrical primitives and a gaussian convolution filter. The surface is also easily modifiable to produce anatomical anomalies such as stenoses or aneurysms.

Given the previous analytic formulation of the surface, the collision detection is about to find where a point cross the surface. With $g(P) = f(P) - \text{iso}$, this is equivalent to root-finding problem on the interval $[P_t, P_{t+1}]$. A Newton-Raphson algorithm is used to grant quadratic convergence. Moreover to insure convergence criteria (i.e. monotonicity of $g$ and a single root on the interval), we subdivide $[P_t, P_{t+1}]$ according to the number of primitives and their radii. When the crossing point $P$ is found, we can easily compute an linear approximation of the surface as $-\nabla g(P) \times x = \|P\|$ (tangent plane) to use it for the contact algorithm. Figure 2.5 illustrates several steps of a simulated procedure where a catheter and a guide-wire are inserted into a virtual simple vascular model.

We have applied our method to a complex simulation: the navigation of a catheter and guide-wire inside a reconstructed vascular network to perform a virtual angiography and the deployment of a coil inside an aneurysm to perform an embolization. The catheter, guide-wire and coil models consist of a series of non-linear deformable
<figure><img src="image_40.png" /><figcaption>**Fig. 2.4:** Skeleton-based convolution surface: using a skeleton defined by cylindrical-shaped primitives (top, for clarity, only the center line of the cylinders are shown) and a Gaussian convolution filter. The potential field is displayed with a heat-map (left-bottom).</figcaption></figure>

beam elements. Each device is composed from 100 to 200 beam elements while the vascular network is constituted of more than 4,000 vessels and undergoes periodic deformations due to both cardiac and respiratory motions. When entering the cerebrovascular system, where the diameter of the vessels is very small, the catheter or guide-wire are constantly colliding and sliding along the vessels wall (see Figure 1). Similarly, when deployed within an aneurysm, the coil becomes highly constrained, and proper contact modeling becomes of prime importance to guarantee a correct behavior during the simulation (see Figure 2.6). We have performed a series of simulations on a standard desktop computer and obtained real-time computation rates (100 Hz). These timings include the computation and inversion of the system stiffness matrix K at each time step, as well as collision detection and collision response. Since the contacts are solved in the contact space, the size of the system is the number *c* of contact (defining *n* as the number of DOFs, *c* ≤ *n* and usually *c* << *n*). It is also important to mention that, in order to enforce the convergence and stability of the contact algorithm, we use a subdivision strategy where each time step is subdivided into a variable number of sub-steps. The initial time step is subdivided if not all contacts have been solved after *N* iterations of the main loop of the algorithm. This subdivision strategy allows us to solve complex contact configurations and to handle concave cases with a succession of convex cases.
<figure><img src="image_41.png" /><figcaption>Fig. 2.5: Example of catheter and guide-wire navigation inside a vascular phantom using convolution surfaces as the underlying geometric model.</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>**Fig. 2.6:** Example of catheter and guide-wire navigation inside a vascular network of a specific patient. The vascular model consists in more than 4000 blood vessels and the simulation can be run at interactive rates (100 Hz). Complex configurations (tortuous and small vessels) are handled correctly.</figcaption></figure>
<figure><img src="image_43.png" /><figcaption>**Fig. 2.7:** Catheter deployment in a closed capsule. (Left-column) the capsule is modeled as a polygonal mesh, (Right-column) the capsule is modeled as a convolution surface. The snapshots taken at the same simulation time-steps illustrate that artificial friction can be induced with polygonal meshes because the tip sticks into the perfectly smooth vessel instead of sliding along it (which is the case of the implicit surface).</figcaption></figure>

To illustrate how artificial friction can be induced with (relatively) polygonal models, we perform a simple experiment: we deploy a catheter-like tool inside a closed capsule with an ellipsoidal shape. This shape is modeled with 400 polygons (mainly quads but with several triangles) or with a skeleton-based convolution surface with a single primitive (**a segment**). The results are illustrated on 2.7 where the left-column displays the deployment inside the polygonal mesh while the right-column displays the deployment inside the implicit surface. Deployment inside the implicit surface is perfectly smooth without any hitch. On the contrary, the deployment inside the polygonal surface exhibits a moment where the tip of the catheter is stuck in the surface (on an edge between two polygons), then the deployment produces a spring catheter shape where a lot of internal energy is accumulated that can lead to fast and uncontrollable release at any time. Again such artificial friction is unnoticeable when objects in simulation have temporary contacts but when the simulation will contain sliding of objects on another, this becomes a major concern. Obviously this limitation of polygonal surfaces can be addressed by using smaller elements and thus producing a *smoother* surface but it will significantly improve the time needed for collusion detection.

This previous examples illustrate the potential of implicit surfaces for simulation and more specifically for applications where contacts with friction and sliding will be preponderant such as interventional radiology procedures are. However while it is straight-forward to derive implicit surfaces for basic shapes or skeleton-based
shapes, building an implicit surface from an arbitrary shape is complicated. Since we aim at patient-specific simulation (for planning applications for instance), we should address how to build accurate implicit surfaces for medical images.

## 2.2 Models of vascular anatomy

Our objective is the simulation of neuroradiology interventions. Procedural implicit surfaces appear best suited for precise and physically grounded contact computation. But the complex intracranial vasculature is likely to require a complex implicit function, putting interactive simulation rates at risk. Besides, vessel centerlines carry precious topological information to accelerate contact localization. But current implicit vessel representations cannot be used for simulation. Cylinder-like shapes are unable to capture details at the vessel surface, and cross-sectional descriptions suffer from a lack of smoothness in the longitudinal direction.

The Computer Graphics community has devoted a vast literature to the problem of implicit surface reconstruction from data points [Ber+13]. Such methods bear no difficulty to interpolate missing data points [Oht+05] but handling noise and outlier data points is still considered an issue [Mul+10]. Besides the advantages mentioned in the previous section, implicit surfaces tend to blend as they come close to each other. As a result, implicit modeling of a complete vasculature is very difficult since, for example, close or kissing vessel branches merge 2.8. Controlling the blend is still an open problem when more than two objects are merging or in the case of complex topologies. Moreover using generic-purpose segmentation algorithms may lead to other artifacts such as holes or discontinuing vessels (see Figure 2.8 bottom).

Therefore, we propose in this section a new blood vessel surface model, namely Local Implicit Modeling (LIM), that uses a specific implicit reconstruction from data points. To get these data points, a new blood vessel tracking algorithm is proposed that

* detects the vessel centerline;

* provides a local radius estimate;

* extracts a dense set of points at the blood vessel surface.

This algorithm is based on a RANdom SAmple Consensus (RANSAC) based robust fitting of successive cylinders along the vessel. Our method was validated against the Multiple Hypothesis Tracking (MHT) algorithm on 10 patients three-Dimensional
<figure><img src="image_44.png" /><figcaption>**Fig. 2.8:** Examples of reconstruction issues from segmented medical images. Since the contouring of the vascular structure do not encode any topological information many artifacts may occur such as blending between aneurysms and vessels (top) or holes and discontinuing vessels (bottom).</figcaption></figure>
(3D) Rotational Angiography (RA) data. Over 744 blood vessel of various sizes were considered for each patient. Given these data, LIM locally reconstructs the blood vessel surface thanks to information provided by its centerline and a dense sampling of the vessel surface (point-set). In other words, LIM locally performs the point-set reconstruction by fitting a skeleton-based implicit function - i.e. Blobby Model (BM) - through energy minimization, alternating with an original selection and subdivision scheme of BM primitives (iterative process). The idea of LIM is to enrich the blood vessel centerline by attaching a BM to a point on the centerline that locally reconstructs the vessel surface. The resulting blood vessel surface is supplied as a tree of local implicit surfaces generated by BMs.

## 2.2.1 Vessel Centerline Tracking

### Overview of the algorithm

We model the local vessel shape as a finite cylinder (C, r, d, h) where C is the center, r the radius, d the axis direction, and h the height. Our algorithm alternates between a prediction and an estimation step. State t is predicted from state t − 1 by:

$$C_t^* = C_{t-1} + s \frac{h_{t-1}}{2} \mathbf{d}_{t-1} \qquad (2.3)$$

and

$$(r_t^*, \mathbf{d}_t^*, h_t^*) = (r_{t-1}, \mathbf{d}_{t-1}, h_{t-1}) \qquad (2.4)$$

where s is a scaling factor for the height-dependent step size. We set it to 0.5 in order for successive cylinders to superimpose.

The estimation step is more complex. First, it extracts a set of candidate points at the vessel surface in the vicinity of the cylinder location predicted by ($C_t^*, r_t^*, \mathbf{d}_t^*, h_t^*$). Then the best cylinder ($C_t, r_t, \mathbf{d}_t, h_t$) is found that fits these points. A subset of candidate points are thereby found to lie at the blood vessel surface. We propose to add this subset of selected points to the vessel shape description. Hereafter, a focus is made on the different steps of the estimation process.
<figure><img src="image_45.png" /><figcaption>**Fig. 2.9:** Outline of RBT vessel tracking algorithm. (top left) Cut plane through the original 3DRA data; (top right) Prediction: Starting from $C_{t-1}$ (green dot) and $\vec{d}_{t-1}$ (white arrow), the new center $C_t^*$ (red dot) is found. Candidate points (white dots) are extracted at the vessel surface by casting rays in the volume, from $C_t^*$. Only the points around the displayed cut plane are shown. (bottom left) $N_d$ directions are tested for the cylinder axes (yellow arrows). One fitting cylinder is found per axis, using RANSAC. (bottom right) Best cylinder found ($C_t, r_t, \vec{d}_t$), together with its consensus set (green dots). The red dots were used to compute the cylinder height $h_t$: they are the set of points $P_t$, attached to the cylinder.</figcaption></figure>
## Candidate point extraction

In angiographic images, such as 3DRA data, the vessels are bright tubular-shaped objects, on a dark background. The points of minimal directional gradient along rays cast from the cylinder center $C_t^*$ are collected as candidate points at the local vessel surface (one point per ray). This extraction scheme is very similar to taking the points of maximal gradient along the columns of a bounded spherical projection image [WC07a] except that we used $N_r$ rays evenly sampled on the Gaussian sphere to avoid point clustering that occurs at the poles with spherical coordinates.

The prediction step, followed by the candidate point extraction is illustrated in figure 2.9, top right.

## Cylinder fit with RANSAC filtering

Most of the extracted points lie at the local vessel surface, even though this location is tainted with noise, but some points are outliers, due e.g. to neighboring structures or bifurcations. RANdom SAmpling Consensus (RANSAC) is a robust estimation procedure that has been successfully applied to fit cylinders to a point set [CG01]. It relies on a closed-form equation to retrieve the cylinder parameters from a minimal set of points. Various solutions exist depending on the parametrization of the cylinder [CW06]. In [CG01], the cylinder axis direction was first estimated and the center and radius of the cylinder were estimated using RANSAC. This parametrization has the advantages of simplicity (it resorts to fitting a circle to three points in the plane orthogonal to the axis) and low complexity (the minimal input set has only 3 points).

In our case, the axis direction can hardly be estimated from the point normals as in [CG01]. Indeed, the normals would be computed from image gradients, which are sensitive to noise, and would be untrustworthy around curved portions of the vessels. Therefore, a strategy similar to [Fri+10] was adopted: $N_d$ predefined axis directions are tested, centered on $\vec{d}_t^*$; the best cylinder (center and radius) for each direction is estimated by RANSAC (see figure 2.9, bottom left). The $N_d$ directions are reviewed by increasing angle with $\vec{d}_t^*$: as soon as a valid cylinder is found (see below), the search stops and the cylinder is kept. This pruning strategy speeds up the process.

In order to cope with a large scale of vessel sizes, the distance threshold was set proportional to the estimated radius: a point is an inlier if its distance to the cylinder is below $0.1r_t^*$. The set of inlier points is called the consensus set.
## Validity of a cylinder

RANSAC iteratively fits a cylinder to a set of 3 randomly sampled points. This cylinder is considered valid as soon as its consensus set encompasses more than $p_{int}$% of the input point set. The density constraint proposed by [CG01] was not applied in order to consider cylinders despite holes in the point set due to e.g. vessel bifurcations. However, RANSAC may exhibit cylinders of extreme radii (either too small or too large) as degenerate solutions. Therefore, a cylinder was considered as valid only if its radius was in $[r_t^*/2, 3r_t^*/2]$ [WC07a].

## Refinement and height estimation

The minimization, using Powell's algorithm, of the distance from the returned valid cylinder to its inlier point set is searched to tune its parameters. In particular, it enables the axis direction to be refined from its initial discrete value. Then, the position of the cylinder center is adjusted along the axis, to coincide with the median position of the inlier points. Thereby the new parameters ($C_t, r_t, \vec{d}_t$) are set. A new set of candidate points is extracted using these values and the points in consensus with the cylinder are kept: the height $h_t$ is set so that 75% of the inlier points are encompassed around $C_t$. These inlier points are returned by the estimation procedure. We will denote this set as $\mathcal{P}_t$ (see figure 2.9, bottom right).

RBT produces a representation of the vasculature as a tree of center points $C_t$, each equipped with a tuple {$r_t, \vec{d}_t, h_t, \mathcal{P}_t$}. We call this information a *node*: each node gives a local picture of the blood vessel around its center.

## Tracking the whole vessel tree

The previous paragraphs described one inner loop of our RANSAC-Based Tracking algorithm (RBT). The correction step of RBT, where a cylinder is fitted to candidate points, proved to be very robust to poor initialization of the cylinder parameters. As an example, the trackings for brain arteries, were initialized with the up vector direction $\vec{d}_0 = (0, 0, 1)$, the average radius reported in the literature ($r_0 = 2mm$) and $C_0$ manually placed in the proximal portion of the internal carotid artery. The tracking started with the correction step.

Moreover, RBT can track a single branch at a time, but it does not stop at bifurcations. Instead, it recursively takes one bifurcating branch and tracks it until its distal end. According with the algorithm properties, this branch depends on its alignment to the
<figure><img src="image_46.png" /><figcaption>**Fig. 2.10:** Adding a new branch: (left) Initialize with a user-input point (green point), project it on the parent centerline (violet point). Their connection initializes the vessel direction (arrow). Tune the initial (green) point, by moving it along the vessel direction, so that it is out of the parent vessel, but still close to the proximal end of the child vessel (blue point). (right) Tracks the child branch from this point and connects it to the parent branch (stale blue polyline)</figcaption></figure>

prior tracked centerline, its contrast and clarity, but may also be random, depending on the RANSAC sampling. To add a new branch, the user was asked to manually place a starting point $C_0^*$ (see figure 2.10). This initial location was automatically tuned by projecting $C_0^*$ onto the parent vessel centerline as $C_0'$: the initial direction was $\vec{d}_0 = (C_0^* - C_0')/|C_0^* - C_0'|$ and the initial radius $r_0$ was set equal to the radius of the parent vessel. The initial point $C_0$ was set so that $C_0 = C_0' + 1.25r_0\vec{d}_0$.

## Capturing aneurysms

Aneurysms do not comply well with the cylinder shape prior. A manual procedure was chosen to capture aneurysms, following [WNV00]: center points are placed within the aneurysm (1 or 2 points are enough in most cases), and rays are cast from each point, keeping the first location where the directional gradient is above a user-defined threshold to build the corresponding $P_t$. The radius is estimated by fitting a sphere to each point set $P_t$, in the least-squares sense, and center points are linked to complete the vessel centerline. The subtree corresponding to the aneurysm is attached to the RBT tree.
## 2.2.2 Implicit modeling: The LIM algorithm

RBT models the vascular tree as a tree of local point sets, around center points associated with a vessel radius estimate. This depiction is too rough and approximate for an accurate and realistic tool simulation to be conducted. The central idea of local implicit modeling (LIM) is to further enrich each node on this tree with a continuous and smooth representation of the local vessel surface as an implicit surface. Figure 2.11 displays the main steps in the reconstruction process for one node in the vicinity of a bifurcation. The same process is applied to each node in the vascular tree.

### Implicit formulation

Implicit surfaces generated by skeletons are particularly well-suited for modeling smooth free-form shapes in a very compact and versatile way [BTC95]. They present good locality properties: skeleton objects far apart have no, or at least negligible, influence on each other [DOG04]. In particular, point set skeletons produce procedural implicit functions which are computationally very efficient.

An implicit iso-surface generated by a point-set skeleton [Bli82] is expressed as the zero-level set $S$ of a function $f$, a sum of implicit spheres:

$$f(X; p) = \sum_{j=1}^{N_b} \alpha_j \phi \left( \frac{|X - C_j|}{\rho_j} \right) - T \quad (2.5)$$

where $X \in \mathbb{R}^3$, $T \in \mathbb{R}^+$ is the iso-surface threshold, $\{\alpha_j\}$ are positive weights (height), and $\{C_j\}$ is the point-set skeleton. Each implicit sphere $\#j$ is defined by a symmetric spherical function, centered on $C_j$, of width $\rho_j$ and following a radial profile. This profile (or kernel) is a function $\phi: \mathbb{R}^+ \to \mathbb{R}^+$, rapidly and monotonically decreasing to 0 and such that $\phi(0) = 1$, and $\phi$ has a bell-shaped graph [She99]. $C^1$ continuity around the center point imposes that $\phi'(0) = 0$. Though our implementation supports several kernels, all our experiments were realized using the Cauchy kernel: $\phi(x) = (1 + x^2/5)^{-2}$.

Muraki [Mur91] was the first to use a Gaussian-based model in the context of object reconstruction. Such models were called differently depending on the kernel used [She99]. Following the seminal work of [Mur91], we shall use the terms *blob* for an implicit sphere, and *Blobby Model (BM)* as a generic name for the implicit models $f$.
<figure><img src="image_47.png" /><figcaption>**Fig. 2.11:** Outline of the LIM algorithm. (top left) Input data: RBT for two branching arteries on 3DRA data results in points on the centerlines (blue and black points) equipped with points locally extracted at the vessel surface (orange dots); (top right) The red node is chosen. (middle left) Closest neighbor nodes are selected w.r.t a given topological distance (of one, here). (middle right) Local point sets are concatenated. (bottom left) The local implicit surface (blobby model) is initialized with two close spheres. (bottom right) LIM successively fits and splits blobs to fit the input points. 30 subdivisions were necessary to produce the resulting implicit surface.</figcaption></figure>
In our particular simulation context, in order to help predict collisions, and have the function give a valid interpenetration direction, the algebraic value $f(X; p)$ at point X should relate monotonically to the geometric distance of X to the surface. We set $\alpha_j = \rho_j$, which obviously establishes the sought relation in the case of a single blob. Furthermore, redundancy in the implicit formulation of $f$ is dismissed since height $\alpha_j$ and width $\rho_j$ are dual in their role. As a consequence, the implicit function $f$ is parametrized by $p^T = \{C_j^T, \rho_j\}_{j=1...N_b}$, the vector of dimension $4N_b$, gathering centers and widths of the $N_b$ implicit spheres composing the implicit function.

## Energy Formulation

Fitting a surface to $N_p$ points $\{P_i\}$ can be written as an energy minimization problem [Mur91; BTC95; TBC95]. We propose to combine 3 energy terms: $\mathcal{E} = \mathcal{E}_d + \alpha\mathcal{E}_c + \beta\mathcal{E}_a$ where $(\alpha, \beta) \in \mathbb{R}^{+2}$, and:

* $\mathcal{E}_d = 1/N_p \sum_i f(P_i; p)^2$ translates the algebraic relation between data points and the zero-level set. It gives a raw expression of the approximation problem.

* $\mathcal{E}_c = 1/(N_b(N_b-1)) \sum_{j \neq k} \left( \frac{s\sqrt{\rho_j\rho_k}}{|C_j-C_k|} \right)^{12} - 2 \left( \frac{s\sqrt{\rho_j\rho_k}}{|C_j-C_k|} \right)^6$ is Lennard-Jones energy. Each term is minimal (with value -1) for $|C_j - C_k| = s\sqrt{\rho_j\rho_k}$. $\mathcal{E}_c$ is repulsive for blobs closer than this distance, and attractive for blobs further away. It imposes some cohesion between neighboring blobs to avoid leakage where data points are missing, while preventing blobs from accumulating within the model.

* $\mathcal{E}_a = 1/N_p \sum_i \kappa(P_i)^2$ $\kappa(P)$ is the mean curvature. It can be computed in a closed form at any point in space from the implicit formulation [Gol05]

$$ \kappa(P) = \frac{\nabla f^t H_f \nabla f - |\nabla f|^2 \text{trace}(H_f)}{2|\nabla f|^3} $$

where $\nabla f$ is the implicit function gradient and $H_f$ its Hessian matrix, both computed at point $P$. This energy smoothes the surface according to the minimal area criterion. In particular, the wavy effect that could stem from modeling a tubular shape with implicit spheres, is reduced.

Behind the rather classical form given above for the energy terms, it is important to notice that the whole energy is known under a closed-form expression. As a
consequence, closed-form expressions were derived for its gradients with respect to the BM parameters {$ρ_j$} and {$C_j$}.

## Selection-Subdivision

Finding the appropriate number of blobs is very hard. Starting from one or few blobs, a subdivision scheme enables to progressively reach a compromise between a good model accuracy and a low complexity. The blob subdivision procedure proposed in the seminal work of Muraki was exhaustive and time consuming [Mur91]. A blob selection mechanism was added in [TBC95], measuring the contribution of each blob to $E_d$ in a user-defined window, and choosing the main contributor. User input is not an option in our context where thousands of BMs should be handled. Moreover, we experimentally noted that this technique was prone to favor small blobs, thus focusing on details, before dealing with areas roughly approximated by one large blob. This behavior is caused by the selection mechanism using the algebraic distance to the implicit surface. Instead, our criterion relies upon the geometric distance approximation proposed by [Tau91] (see figure 2.12, middle left). First, we find the point $P_{i*}$ farthest to the surface:

$$i^* = \arg \max_{1 \le i \le N_p} \frac{|f(P_i; p)|}{|\nabla f(P_i; p)|} \qquad (2.6)$$

Next, we select the blob $\#j^*$ whose isosurface is the closest to $P_{i*}$:

$$j^* = \arg \min_{1 \le j \le N_b} \frac{|f_j(P_{i*})|}{|\nabla f_j(P_{i*})|} \qquad (2.7)$$

with

$$f_j(P_{i*}) = \rho_j \phi \left( \frac{|P_{i*} - C_j|^2}{\rho_j^2} \right) \qquad (2.8)$$

We thereby select the blob $\#j^*$ whose optimization will most likely reduce the error at the worst fitted point $P_{i*}$.

The subdivision step then replaces this blob with two new ones (figure 2.12, middle right). Their width $ρ'_{j*}$ is chosen such that two blobs, centered on $C_{j*}$, of width $ρ'_{j*}$ would have the same isosurface as one blob centered on $C_{j*}$, with width $ρ_{j*}$. The first new blob is centered on $C_{j*}$, while the second is translated by $ρ_{j*}/10$ towards $P_{i*}$. This small-shift procedure stabilizes the cohesion energy $E_c$ which otherwise would become large when two blobs occupy the same location.
<figure><img src="image_48.png" /><figcaption>Fig. 2.12: The elementary reconstruction step in 2D. (a) Points (black) are to be reconstructed with two blobs (red contour). (b) First, all blobs are fine tuned over the centers and then, over the widths. (c) Second, find the point farthest to the contour (cyan) and the blob closest to it (purple). (d) Third, replace this blob with two blobs. The first (black contour) is placed at the same position (black) and the second (cyan contour) is translated in the direction of the farthest point (cyan). Widths are calculated so that both blobs centered at the black point generate the same isosurface as the fissioned blob (purple). (e) Fourth, the translated blob (cyan) is tuned over its center then width. (f) Finally, the BM is updated by integrating the recently tuned blob.</figcaption></figure>
## Optimization

Such an incremental subdivision procedure may lead to a dramatic increase in the number of blobs, and hence the size of the optimization problem. The locality of the kernel $\phi$ allows us to focus the optimization onto the newly created pair of blobs (figure 2.12, bottom left). More exactly, only the new blob that is slightly misplaced is optimized, the other blobs remaining constant. The energy is minimized using Polak-Ribiere (PR) gradient descent algorithm, taking advantage of the closed-form expressions of both the energy and its gradients. A single minimization loop consists in one PR minimization over the center (3 variables), followed by one on the width (1 variable). The optimization stops as soon as a maximum number of loops or a threshold accuracy $t_g$ (as measured by Taubin's distance, see eq 2.6) is reached. In practice, a maximum of 5 loops proved sufficient. Threshold $t_g$ will be discussed in Section 2.2.3.

## Cleaning

Redundant blobs are either ejected far away from the surface during the optimization over the centers, or their width is reduced to almost zero during the minimization over the widths. A simple clean-up procedure was applied after each minimization loop: blobs relatively far from the node center (above 20 mm), and blobs whose width was below a certain size (0.02 mm i.e. one tenth of the voxel size) were removed from the BM. In combination with the subdivision process, this clean-up enables the algorithm to simulate large blob displacements, which is hardly possible with a local minimization, thereby adding robustness to poor initialization.

## Initialization and global modeling

The vascular tree segmentation obtained with RBT is taken as input to our algorithm as a tree of nodes. In order to ensure smooth transitions between adjacent node BMs (figure 2.11, top right), the local data point sets $P_t$ are concatenated so as to create overlaps between neighbors. The input point set considered to be fitted by LIM is therefore $\{P_i\}_{t = \bigcup_{k \in N_t} P_k}$ with $N_t$ a neighborhood of the node of interest, defined by all nodes within a maximum topological distance $\mathcal{T} = 2$ in the tree in our implementation.

First, the initial BM is constituted of two blobs placed on either side of the node center $C_t$ along the vessel direction $\vec{d}_t$ so that an elongated shape is formed. The
<figure><img src="image_49.png" /><figcaption>**Fig. 2.13:** Implicit modeling of an aneurysm. The points {$P_i$} are in red. (From left to right) Initialization with a single blob; after the first minimization; after 25 subdivisions; final result (100 subdivisions)</figcaption></figure>

width $\rho_t$ is scaled to match the shape of a sphere with radius $r_t$ (see figure 1.14, bottom middle. This match depends on $\phi$).

Then, an implicit surface which fits the point set {$P_i$}$_t$ is produced through selection, subdivision and optimization of the blobs as described above. The process is repeated until the maximum distance of the blobby model to the point set falls below threshold $t_g$ or a maximum number $N_s$ of blobs is reached. Each local BM is treated independently of others, thus enabling parallel computation. Each initial node is enriched with a local BM $B_t$, resulting in a global model of the vasculature as a tree of implicit functions.

Figure 2.13 illustrates the essential steps of this fitting algorithm and its capacity to model complex shapes, such as aneurysms, even from rough initialization.
### 2.2.3 Results

To show the strength of the proposed framework, we decomposed our assessment in two parts. First, we compared RBT against state-of-the-art Multiple Hypothesis Testing (MHT) algorithm [Fri+10] on clinical data. Next, the LIM algorithm was tested on synthetic and clinical data. A more complete description of the experimental protocol is available on our journal publication in *Medical Image Analysis* [Ker+17].

We considered 10 patients scheduled for the endovascular angioplasty of an intracranial aneurysm at the Department of Interventional Neuroradiology, University Hospital of Nancy, France. Each patient data set consisted of a 3D Rotational Angiography (3DRA) acquired on a vascular C-arm (Innova 4100, GE Healthcare) during the intra-arterial injection of the internal carotid artery. 3DRA volumes presented as a 512³ isotropic voxel cubes, with a voxel size between 0.18-0.22 mm.

#### Evaluation of the RBT algorithm

**Assessing the tracking ability** A tracking was considered as successful if the resulting centerline remained inside the vessel of interest. All 744 vessels were visually inspected to compute this Success Rate (SR) in percents.

The MHT tracking failures were then visually, manually cut before the problem occurred (see figure 2.14 for an example). No such cut was ever required for RBT. The lengths of the centerlines were computed and summed up to give the Total tracked Length (TL) in millimeters.

Each vessel was thereafter classified as: short, medium, or long, if the centerline was respectively correct along less than 33%, 33-66%, more than 66% of the visually estimated length of the targeted vessel. The success rate was recomputed, considering as successes only medium and long vessels. We denote this score SRL in percent.

Then the Common path Length (CL, in millimeters) could be computed as the total tracked length on those vessels classified as medium or long by both algorithms.

We assessed the distance between two centerlines with a Modified Average Symmetrical Distance (MASD):

$$MASD = \frac{1}{N' + M'} \left\{ \sum_{i}^{N} d_{w}(P_{i}, Q) + \sum_{j}^{M} d_{w}(Q_{j}, P) \right\} \quad (2.9)$$
<table><caption>Tab. 2.1: Evaluation of RBT compared to MHT. SR and SRL are given in percents, TL and CL are given in mm. The mean value was computed over the 10 patients.</caption><thead><tr><th rowspan="2"></th><th colspan="4">MHT</th><th colspan="4">RBT</th></tr><tr><th>SR</th><th>SRL</th><th>TL</th><th>CL</th><th>SR</th><th>SRL</th><th>TL</th><th>CL</th></tr></thead><tbody><tr><th>mean</th><td>64.9</td><td>68.8</td><td>778.6</td><td>680.9</td><td>93.5</td><td>88.7</td><td>1446.1</td><td>727.6</td></tr><tr><th>min</th><td>48</td><td>57</td><td>555.4</td><td>389.0</td><td>81</td><td>77</td><td>833.9</td><td>403.6</td></tr><tr><th>max</th><td>74</td><td>77</td><td>1086.9</td><td>1074.1</td><td>100</td><td>99</td><td>2373.6</td><td>1164.7</td></tr></tbody></table>

where $P = \{P_i\}_{i \in [1,N]}$ is the centerline produced by MHT, and $Q = \{Q_j\}_{j \in [1,M]}$ is produced by RBT. $d_w$ is the classical Euclidean distance between a point and a curve, but thresholded by $w$: it returns 0 whenever the Euclidean exceeds $w$. Consequently, $N'$ and $M'$ are the number of points respectively in $P$ and $Q$ where the distance was not thresholded. The MASD measure automatically selects the portions common to both centerlines and computes the mean distance between them (see last row of figure 2.14).

Two values were used for $w$: $w = 3$ and $w = \bar{r}$, the average radius of the tracked vessel. Similarly low results would prove the capacity of the measure to consistently select the correct common portions of the vessel centerlines.

Table 2.1 compares MHT to RBT according to the above measures, averaged over the 10 patients. Mean values are reported along with their minimum and maximum values.

RBT displays a SR of 93.5% compared to only 64.9% for MHT. This success rate is confirmed by the SRL, with two precisions. First, a slight increase was observed for MHT (68.8%) since the results were edited. It reintroduced vessels which, for example, were well tracked but for MHT failing to stop (see e.g. figure 2.14, top row, middle). Second, a small decrease in SRL was observed for RBT, confirming its tendency to sometimes stop prematurely. The 20-25% difference in the success rate scores was due to MHT often failing to track a vessel that came close to another dense structure, so-called Kissing Vessel (KV) issue, and to leak into that second structure. Figure 2.15 illustrates the difficulty of MHT to sometimes handle KV, while RBT demonstrates a high robustness to that problem and can therefore track even very tortuous vessels.

When both algorithms succeeded, on average, RBT went further (CL = 725.5 mm) than MHT (CL = 680.9 mm). Nonetheless, MHT definitively delivered a lesser number of successful trackings which impacted negatively the average tracked length per patient (TL) down to 778.6 mm against 1446.1 mm for RBT.
<figure><img src="image_50.png" /><figcaption>**Fig. 2.14:** Top: (left) Maximum Intensity Projection of one patient data. Centerline de-lineation produced by MHT (middle) and RBT (right). Bottom: (left) When the MHT centerline (red) does not remain inside the vessel wall, it is cut – for TL computation – before the problem arises (green). (middle) MHT centerline (red) and the points (cyan) selected for MASD computation ($w = \bar{r}$). (right) The resulting points for MHT (red) and RBT (blue, all the original points) used for MASD computation.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**Fig. 2.15:** (left) A tortuous branch of the posterior choroidal artery (0.6 – 0.7 mm of diameter) in Patient 3. (middle) MHT is misled by neighboring structures in the proximal portion of the artery whereas (right) RBT tracks the tortuous vessel and successfully handles the Kissing Vessel (KV) issue.</figcaption></figure>
The average MASD between both centerlines, with both 3 voxels and the mean radius of the targeted vessel, was 0.18 mm (down to 2 figures of precision) and was below one voxel for all patients (min=0.16 mm, max=0.21 mm with a voxel size between 0.18 and 0.22 mm in our patient data set). This showed that RBT produced similar results as MHT.

## Evaluation of our LIM algorithm

After applying RBT to the carotid artery of one patient, one node of the centerline was considered. The energy weights were tuned to provide energy terms of the same order of magnitude. The isolevel value *T* was chosen to be approximately half the voxel size. Other values were chosen to produce visually good results. In particular, for Lennard-Jones energy, *s* was set to 2: Two blobs with the same width *ρ* are thereby separated by 2*ρ* at the minimum for this energy. Values are summarized in Table ?? and where used for all synthetic and patient tests.

**Synthetic experiment** Four synthetic shapes were used in a first experiment. They consist in skeleton-based convolution surfaces:

* a capsule (C) based on a 4 cm segment skeleton

* the same capsule opened at both extremities (oC)

* an arc (A) based on a 4 cm radius half-circle

* a bifurcation (B) composed of two portions of arcs.

These shapes are displayed in blue on figure 2.16. There were tuned so that their width (distance to their skeleton) was 1 mm. Each shape was sampled on 300 points using an equi-distributed point-based sampling technique [witkin94]. These points served as input to the LIM algorithm (point set called $P_t$ in section 2.2.2).

**Protocol** Each point in $P_t$ was perturbed by a white Gaussian noise along the normal to the reference surface. The standard deviation of the noise varied from 0.025 to 0.25 mm (step 0.025, 10 levels of noise). LIM was applied 30 times for each noise level, with $t_g$ set to 3 times the standard deviation. $N_s$ was set to a large value (100) in order to evaluate the capacity of LIM to produce compact models (i.e. blobby models with a reduced number of blobs). LIM initial BM contained 2 (C and oC) or 5 (A and B) blobs. The accuracy of LIM was evaluated by robust statistics on the Euclidean distance between $P_t$ (noise free version) and the resulting blobby
<figure><img src="image_52.png" /><figcaption>**Fig. 2.16:** Reference surface (blue) and the resulting BM (purple) for maximal values in Table 2.2. The 300 noisy input points are in red.</figcaption></figure>

model: The output BM was sampled on 1000 points (again using [witkin94]) to compute the distance of each point in $P_t$ to the BM and the 90th percentile was taken: this score is denoted $d$. The number of blobs in the output BM was also recorded.

LIM relies on the distance approximation to an implicit surface proposed by [Tau91]. We took here the opportunity to evaluate the adequacy of this approximation and computed it together with the Euclidean distance. The 90th percentile is also reported, and denoted $d_a$.

**Result** Table 2.2 reports the minimum, median and maximum of both $d$ and $d_a$ for the 300 tests, as well as the number of blobs in the final BM.

Excellent distance figures were recorded, with a low median value close to the minimum and a maximum value below the maximum standard deviation of the applied noise. A higher median value was observed for the open capsule oC which exposes a more challenging fitting problem due to its open extremities. Figure 2.16 displays LIM output for the recorded maximum $d$ value. We can see that the points at maximum distance are located in the vicinity of the open extremities. This issue is tackled in our case by letting successive BMs overlap.

The number of blobs in the final BMs is typically low and demonstrate the capacity of LIM to produce compact models. Overly large BMs were seldom obtained when the noise process randomly applied large displacements to some points. LIM is indeed not robust but the filtering task is devoted to RBT.

The approximate distance $d_a$ gave very similar results to $d$ except for large values which it tended to underestimate. The ordering of points according to their distance to the implicit surface was nevertheless always the same with both distances, which makes the use of the approximate distance appropriate in LIM.
<table><caption>**Tab. 2.2:** Evaluation of LIM accuracy and compacity on 4 synthetic shapes (capsule C, open capsule oC, arc A and bifurcation B) perturbed by a Gaussian noise. *d* is the 90th percentile of the Euclidean distance of noise free points in *P*<sub>t</sub> to the output BM. *d*<sub>a</sub> presents the same statistics using the approximate distance proposed by [Tau91]. *N*<sub>b</sub> is the number of blobs in the output BM. Minimum (min), median (med) and maximum (max) values are reported over 300 tests.</caption><thead><tr><th>Shape</th><th></th><th>C</th><th>Co</th><th>A</th><th>B</th></tr></thead><tbody><tr><td rowspan="3">d</td><td>min</td><td>0.015</td><td>0.039</td><td>0.016</td><td>0.048</td></tr><tr><td>med</td><td>0.042</td><td>0.119</td><td>0.063</td><td>0.078</td></tr><tr><td>max</td><td>0.189</td><td>0.233</td><td>0.174</td><td>0.334</td></tr><tr><td rowspan="3">d<sub>a</sub></td><td>min</td><td>0.015</td><td>0.038</td><td>0.016</td><td>0.046</td></tr><tr><td>med</td><td>0.041</td><td>0.103</td><td>0.064</td><td>0.079</td></tr><tr><td>max</td><td>0.164</td><td>0.175</td><td>0.146</td><td>0.211</td></tr><tr><td rowspan="3">N<sub>b</sub></td><td>min</td><td>5</td><td>5</td><td>8</td><td>12</td></tr><tr><td>med</td><td>7</td><td>9</td><td>15</td><td>15</td></tr><tr><td>max</td><td>89</td><td>83</td><td>102</td><td>15</td></tr></tbody></table>

**Reconstruction on patient data** Input data was provided by RBT algorithm which was first run on the carotid artery (the stem vessel of the vascular tree) and subsequently, all vessels were tracked (see Section 2.2.1). A total of 379 instances of RBT were needed (between 26-56 instances per patient). Moreover, a total of 11 aneurysms were segmented with the algorithm described in Section 2.2.1. A total of 87564 BMs were to be fitted by LIM on these patient data. Parameters used (see Table ??) offer LIM good convergence properties for nice visual results. We first assessed this visual quality.

The accuracy however depends on two parameters that required more care to be fixed. $N_s$ is the maximum number of subdivisions: the larger $N_s$, the more blobs, allowing for greater modeling versatility; $t_g$ is the maximum error allowed for LIM: this error is used to determine $i^*$ (see eq. 2.6). We tested 3 values for both $N_s$ ({30,50,100}) and $t_g$ ({0.2,0.3,0.5}, in mm). The resulting 9 configurations were evaluated based on precision and compacity of the output model.

**Visual assessment** Local BMs generated by LIM do not assemble into a representation suitable to the visualization of the global vessel tree. A rather good approximation can however be given by first meshing each local BM (we used marching cubes), then cropping each local surface into its Voronoi cell (see Section ?? for a precise definition). Finally all these local surfaces are visualized together for a given BM tree. An example is given in figure 2.17.
<figure><img src="image_53.png" /><figcaption>**Fig. 2.17:** Visual assessment of LIM: (top-left) isosurface of 3DRA data; (top-right) LIM result, showing much more vessels; (bottom-left) no unwanted blending at kissing vessels (encircled); (bottom-right) difficulty to sometimes model branching of a tiny vessel on a large one (encircled), small discontinuities are visible between adjacent BMs (arrow).</figcaption></figure>

LIM was able to reconstruct both tiny and large vessels, and capture well even tortuous vessels. Even though each BM is processed independently from its neighbors, our overlapping neighboring point sets produced smooth transitions between BMs. A visual comparison with a classical isosurface extraction revealed that more vessels are visible with our proposal (see figure 2.17, second to left). Kissing vessels were well handled and no unwanted blending occurred (see figure2.17, second to right). Bifurcations were correctly handled. Some convergence issues were recorded in some places where a very small vessel branches onto a large vessel (e.g. the anterior choroidal artery on the internal carotid artery, see figure 2.17, right).

**Precision** In the absence of ground truth, we resorted to the statistics on the distance approximation called *d*<sub>a</sub> in Section 2.2.3. Our synthetic experiment indicated *d*<sub>a</sub> was a good estimate of the score *d* based on the Euclidean distance.

Blobby Model outputs (87654) were separated into 4 classes based on their *d*<sub>a</sub> score:
* $C_1 \equiv d_a < 0.5 v_s$, *C_2 \equiv 0.5 < v_s \le d_a < v_s$, *C_4 \equiv v_s \le d_a < 2 v_s* and *2 v_s \le d_a*.
The voxel size *v_s* varies between 0.18 and 0.22 mm per patient and is related to the
<table><caption>**Tab. 2.3:** Distribution in % of the BMs according to the error of fit $d_a$ in 4 classes (left column): $C_1$ corresponds to BMs with $d_a < 0.5 v_s$, $C_2$ to $0.5 v_s \le d_a < v_s$, $C_3$ to $v_s \le d_a < 2 v_s$ and $C_4$ to $2 v_s \le d_{bm}$. $v_s$ is the voxel size. 9 algorithm configurations were investigated depending on parameters $t_g$ (targeted precision of fit) and $N_s$ (maximum model complexity). Computations implied 87564 BMs.</caption><thead><tr><th>t<sub>g</sub></th><th colspan="3">0.5mm</th><th colspan="3">0.3mm</th><th colspan="3">0.2mm</th></tr><tr><th>N<sub>s</sub></th><th>100</th><th>50</th><th>30</th><th>100</th><th>50</th><th>30</th><th>100</th><th>50</th><th>30</th></tr></thead><tbody><tr><th>C<sub>1</sub></th><td>51.37</td><td>51.37</td><td>51.37</td><td>88.65</td><td>88.63</td><td>88.54</td><td>97.15</td><td>96.96</td><td>96.58</td></tr><tr><th>C<sub>2</sub></th><td>47.67</td><td>47.64</td><td>47.61</td><td>11.00</td><td>10.96</td><td>10.97</td><td>2.60</td><td>2.73</td><td>3.01</td></tr><tr><th>C<sub>3</sub></th><td>0.84</td><td>0.85</td><td>0.86</td><td>0.25</td><td>0.28</td><td>0.34</td><td>0.15</td><td>0.20</td><td>0.26</td></tr><tr><th>C<sub>4</sub></th><td>0.12</td><td>0.14</td><td>0.16</td><td>0.10</td><td>0.12</td><td>0.15</td><td>0.09</td><td>0.11</td><td>0.14</td></tr></tbody></table>

point extraction noise in RBT. Table 2.3 reports the distribution, in percents, of $d_a$ among the 4 classes.

99% of BMs present sub-voxel $d_a$ ($C_1 \cup C_2$) for all configurations. Moreover, $N_s$ does not appear to have any influence on the precision of fit, which is governed by $t_g$. Besides, a precision of $t_g = 0.2$ mm led to less smooth models where LIM had a tendency to try and adjust to noise in the input point set. $t_g = 0.3$ mm was the best compromise between smoothness and precision.

**Compacity** A measure of compacity is given by the number of blobs in a BM. A first measure of compacity was the total number of blobs in the BMs required to model a patient's vasculature. We found out that this number did not depend on $N_s$ and we hereby report only figures for $N_s = 100$. The average number of blobs per patient was: 33930 ($t_g = 0.5$ mm, average of 3.9 blobs per BM), 56395 ($t_g = 0.3$ mm, average of 6.4 blobs per BM), and 85733 ($t_g = 0.2$ mm, average of 9.5 blobs per patient).

The number of blobs is the number of primitives required to encode the implicit surface. In the case of a triangular mesh, a similar measure is given by the number of triangles. Therefore, we computed for each patient, the ratio between the number of triangles in a mesh representation of the surface, and the total number of blobs in the BMs used to implicitly model the whole vasculature: It indicates the ability of LIM to compress the surface representation. Results are reported in table 2.4. For a targeted precision $t_g = 0.3$ mm ($N_s = 100$), the expected compression ratio is over 100.
<table><caption>**Tab. 2.4:** Ratio between the number of triangles #Δ in a triangulated surface and the number of blobs $N_b$ in the local implicit model for 3 values of $t_g$ ($N_s$ = 100). Minimum (min), median (med) and maximum (max) values are reported for 10 patients.</caption><thead><tr><th rowspan="2">t<sub>g</sub> (mm)</th><th colspan="3">#Δ / N<sub>b</sub></th></tr><tr><th>min</th><th>med</th><th>max</th></tr></thead><tbody><tr><td>0.5</td><td>159</td><td>234</td><td>308</td></tr><tr><td>0.3</td><td>73</td><td>104</td><td>138</td></tr><tr><td>0.2</td><td>44</td><td>58</td><td>77</td></tr></tbody></table>

<figure><img src="image_54.png" /><figcaption>**Fig. 2.18:** Data used for validation: capsule (left) and spiral (right) surfaces. Both synthetic shapes presented respectively dimensions (width×height×depth) in mm: 11.9 × 2.5 × 2.5 and 8.1 × 6.12 × 1.05. A primitive was placed at each red dot composing the skeleton (white).</figcaption></figure>
## Use for interactive simulation

One expected benefit of using implicit surfaces is the reduced collision management time. To assess the relevance of LIM, we consider two synthetic models. A capsule (figure 2.18 left) was synthesized from a 10 mm line segment skeleton sampled on 10 points, with a radius of 1.25 mm. It provides the most simple vessel shape. A spiral (figure 2.18 middle) was generated from 100 points on a curve line with polar equation $(\rho, \theta)$ s.t. $\rho = \frac{\lambda\pi}{\theta^2}$ with $\lambda \in [15.75, 63]$ and $\theta \in [\pi, 4\pi]$. The shape radius varied from 0.15 mm to 0.5 mm. A smooth surface was obtained in both cases as a convolution surface using Cauchy kernel.

Collision with an implicit surface is easily detected through the sign of the implicit function: in our case, a point crossed the surface if $f$, the local BM function, became negative (see Equation 2.5). A node must be selected in the tree output by RBT and LIM to determine the local BM. The node, whose center is the closest to the colliding point, is not appropriate. First, between two time steps, the point could occasionally jump from one artery to another one, e.g. in the case of kissing vessels, which is mechanically impossible. Therefore, after motion, a point is associated to either the same node as before, or one of its topological neighbors in the tree (see figure 2.19).

Second, selecting the closest node resorts to finding the cell the point belongs to in the Voronoi diagram of the node centers. This strategy fails where a small vessel branches onto a larger one: The BM modeling the proximal portion of the small vessel might be selected while the tool point is still inside the large vessel. Using a power diagram enables to adapt the size of the cell to the estimated vessel radius $r$ at the node. This diagram is generated with the power pseudo-distance :

$$d^p(X, C) = |X - C|^2 - r^2 \qquad (2.10)$$

where $X$ is the tool point, $C$ the node center and $r$ the estimated vessel radius at the node. After motion, the point is associated to its closest candidate node according to $d^p$.

At initialization, all the nodes in the tree are candidates for the proximal point of the tool. The $d^p$-closest node is selected for this point. Then, each point is successively associated to its $d^p$-closest node, in the neighborhood of the node selected for the previous point.

The average time taken to detect and solve one collision was measured both for implicit LIM and mesh models. Results are reported in Table 2.5.
<figure><img src="image_55.png" /><figcaption>**Fig. 2.19:** The surface of the local BM selected to solve contact at the tool tip point is displayed in wire-frame for 4 simulation steps. The overall vessel surface is shown in transparent red as a mesh.</figcaption></figure>

Collision management is significantly faster with LIM models than with triangular meshes (about two orders of magnitude). In fact, collision management time does not depend on the total number of primitives for LIM: the capsule have 12k blobs and the spiral 28k blobs but they display similar computation times. This time only depends on the number of blobs in the local selected BM, at most a few tens whatever the complexity of the global shape. BM selection is also done in constant time since we restrict the search to topological neighbors of the current BM.

## 2.3 Models for the pelvic system

The strategy conducted in the previous section shows the benefits of having geometric models dedicated for simulation. Skeleton-based implicit surfaces were convenient in the context of vascular geometry thanks to some assumptions that eases the process of defining the implicit function given a vessel network. However given an arbitrary geometry, the search for an implicit function may be computationally intensive. Other
<table><caption>**Tab. 2.5:** Computation timings for collision detection on various deployment scenarios, using state of the art implementations. For a comparable amount of primitives, LIM outperforms triangular meshes for collision detection and makes it easier to achieve interactive simulation.</caption><thead><tr><th></th><th># of primitives</th><th>Time (msec)</th><th># of contacts</th></tr></thead><tbody><tr><td>Capsule (Mesh)</td><td>12k triangles</td><td>4.89</td><td>51</td></tr><tr><td>Capsule (LIM)</td><td>12k blobs</td><td>0.0232</td><td>50</td></tr><tr><td>Spiral (Mesh)</td><td>2.74k triangles</td><td>1.09</td><td>79</td></tr><tr><td>Spiral (Mesh)</td><td>11.3k triangles</td><td>1.84</td><td>80</td></tr><tr><td>Spiral (Mesh)</td><td>45k triangles</td><td>5.58</td><td>80</td></tr><tr><td>Spiral (LIM)</td><td>28k blobs</td><td>0.0219</td><td>80</td></tr><tr><td>Phantom (Mesh)</td><td>100k triangles</td><td>16.23</td><td>80</td></tr><tr><td>Phantom (LIM)</td><td>42k blobs</td><td>0.133</td><td>82</td></tr></tbody></table>

limitations needs also to be addressed: the possibility to track several organs / soft tissues simultaneously and to track them over time like in dynamic medical imaging (dynamic MRI for instance). Therefore this section uses another representation (B-splines) to represent anatomical structures that will undergo mechanical simulation. This section takes place in the context of female pelvic medicine, where image-based diagnoses of pelvic floor disorders like prolapse or endometriosis are performed. These pathologies rely on mechanical indicators, such as mobilities of organs and shear displacements between organs but image data do not provide directly precise and global values of these indicators. Since the methodology used for this problem is the same as previous, a short overview of the approach and the results will be presented. The reader may refer to [Jia17; Jia+15; Jia+17; Jia+18] for details.

## 2.3.1 Overview of the approach

The chosen approach can be formulated as an optimization procedure. Four key parts are involved in such optimization:

* input data (3D static and 2D dynamic MR images of the patient),

* a mathematical model with variables to be optimized (B-spline Model),

* a cost function, which links the model to the input data (e.g. integration on the domain of the model), and

* an optimizer, which finds the optimal values of the parameters to minimize the cost function.
The optimization process can be described as iterative updates of its DOFs. As we deform the geometric model, we optimize the DOFs that describe the geometrical shape. The geometric model is considered as a mapping from DOFs to the geometry. We introduce two mappings: (1) an affine transformation for the coarse registration and (2) a cubic B-spline as the shape descriptor for the deformable registration. The control points are the researched DOFs. Thus, similar with the typical image registration, our model-to-image registration method is also based on a multi-scale framework.

**Input Data** The input experimental image data consist of T2-weighted static and dynamic MRI of the female pelvic system. The use of these images of volunteers for our research was approved by the National Ethics Committee (Comité d'Ethique de la Recherche en Obstétrique et Gynécologie CEROG 2012-GYN-06-01-R1). The static MR images are obtained in the three planes (sagittal, axial and coronal), which provide information on the 3D anatomical organization. The dynamic MR images (or Cine MRI) refer to a temporal sequence of 2D images, which are widely used for diagnosing pelvic sagging pathologies. For our study, these 2D dynamic images were obtained in the same midline sagittal plane of the patient during pushing, for a given frequency (about 1.2fps). The choice of this middle plane keeps a standard procedure in the clinical routine because it is the most representative of the organ shapes and motions. In order to increase the image contrast of pelvic organs, a gel had been injected into the vaginal cavity and the rectum while patient drank liquid to highlight the bladder. In this way, pelvic organs were easily observable during MRI acquisition. We were interested in 2D segmentation of the three organs (bladder, vagina and rectum), for which both 2D sagittal static and dynamic MR images were used (spatial resolution from 0.47 to 1.17 mm/pixel).

**B-Spline Model** For modeling the geometries of the three organs, we used 2D B-spline of 3-degree that is a classical geometric modeling tool. In our case, each organ is modeled by a closed B-spline $M: u \in [0, 1] \rightarrow M(u) \in \mathbb{R}^2$. Each position on the curve can be calculated by the value of u:

$$M(u, [\mathbf{P}_0, \mathbf{P}_1, \dots, \mathbf{P}_l]) = \sum_{i=0}^{l} \mathcal{N}_{i,p}(u)\mathbf{P}_i \qquad (2.11)$$

where $[\mathbf{P}_0, \mathbf{P}_1, \dots, \mathbf{P}_l]$ are the $l+1$ control points. However, the first and last control point of each curve are attached to form a closed curve: $M(0) = M(1) = \mathbf{P}_0 = \mathbf{P}_l$. Hence, our model is analytical and of $C^2$ continuity except at the end point. Under
<figure><img src="image_56.png" /><figcaption>**Fig. 2.20:** Sagittal dynamic magnetic resonance (MR) image (spatial resolution = 1.17 mm/pixel, 256*256 pixels. (top): First and last images of a dynamic MRI sequence of a healthy volunteer. (bottom-left): dynamic MRI of a prolapse where the uterus moved down. (bottom-right): static MRI of endometriosis</figcaption></figure>
this configuration, the *l* control points are supposed to be the DOFs in the registration algorithm. Another advantage of the B-spline is its locality property, which means each control point is highly independent. In addition, for the affine transformation step, we have the relevant mapping $M_A$, six DOFs of which are evaluated in this coarse step: $[T_{11}, T_{12}, T_{21}, T_{22}, T_x, T_y]$. Suppose $M^0$ is the mapping function corresponding to the initial B-spline curve, and C the center, then the mapping function of the affine transformation becomes:

$$ \mathcal{M}_A(u, [T_{11}, T_{12}, T_{21}, T_{22}, T_x, T_y]) = \begin{bmatrix} T_{11} & T_{12} \\ T_{21} & T_{22} \end{bmatrix} (\mathcal{M}^0(u) - \mathbf{C}) + \begin{bmatrix} T_x \\ T_y \end{bmatrix} + \mathbf{C} \quad (2.12) $$

**Cost Function** In order to define an appropriate cost function, we introduce a virtual image to form the first energy term. A so-called virtual image is generated from the model for calculation of a cost function, which finds the best correlation with the real image. This virtual image is created in the narrow band $\Omega$ along the contour with a certain width. The value of grey levels in this narrow band image is chosen to be similar to the ones in the real image (grey level MR image). This idea was first introduced and proved by Semin et al. [SFA10] to detect fiber-like objects in digital images. This approach was developed then by Réthoré et al. [RF14]. The intensity profile of the virtual image is defined by a levelset function $\mathcal{L}: v \in \mathbb{R} \to \mathbb{R}$ where $v$ is the distance from a position C to the contour (Equation eq 2.13).

$$ \mathcal{L}(v) = \begin{cases} A_{\max} & \text{if } v < 0 \\ A_{\min} + \frac{A_{\max} - A_{\min}}{2} * \left(1 + \cos\left(\frac{\pi v}{\tau}\right)\right) & \text{if } 0 \le v \le \tau \\ A_{\min} & \text{if } v > \tau \end{cases} \quad (2.13) $$

where $\tau$ is the estimated width of narrow band $\Omega$ and $A_{\max}$ and $A_{\min}$ the grey levels inside and outside organ. The value of $A_{\max}$ and $A_{\min}$ are chosen by clicking on the foreground and background in the real image.

$$ \mathcal{F}(D) = \begin{cases} A_{\max} & \text{if } D < -\lambda \\ \frac{A_{\max}}{4\lambda^2}(D - \lambda)^2 & \text{if } -\lambda \le D \le \lambda \\ 0 & \text{if } D > \lambda \end{cases} \quad (2.14) $$
The virtual image $I_Y$ is then computed on the discrete points cloud over the narrow band of the organ border. The discrete points form actually a bended grid, a position $X(u, v)$ of which is defined in the curvilinear coordinate system based on tangent and normal direction of the contour. The virtual image value $I_v(X)$ is computed in the analytical way (Equation 2.15) so that the grid spacing size can even be much smaller than the pixel spacing in real image. Unlike the virtual image, the real one has a pixel-wise square grid; we should then project $X$ and evaluate real image value $I_r(X)$ using a bi-cubic interpolation. Also, for these reasons, the virtual image is considered virtual and provides sub-pixel information compared with the real one. Furthermore, grey level values are merely evaluated in its definition domain $\Omega$ instead of the whole domain of the real image.

$$I_v(\mathbf{X}) = I_v(u, v) = \mathcal{L}(v) \qquad (2.15)$$

The positions $X$ can be evaluated by the analytical formulation of the parametric geometric model (Equation 2.16). Suppose a parametric curve: $\mathcal{M}: u \in \mathbb{R} \to \mathbb{R}^2$.

$$\mathbf{X}(u, v) = \mathcal{M}(u) + v\vec{N}(u) \qquad (2.16)$$

Thus, the computation of the virtual image is more natural and direct with a parametric model $\mathcal{M}$ but can be applied to any type of model. The cost function related to this virtual image for each organ can be computed with:

$$E_{\text{image}} ([\mathbf{P}_0, \mathbf{P}_1, \cdots, \mathbf{P}_l]) = \int_{\mathbf{X} \in \Omega(\text{Band})} [I_r(\mathbf{X}) - I_v(\mathbf{X})]^2 d\Omega \qquad (2.17)$$

Thanks to the virtual image, our method is exempt from image pre-processing and based merely on reliable original data (no need of pre-processing the images). The method requires less computation, and in a mathematical point of view, the virtual image can be considered as a smoothing of gradient field that avoids local optimum in this non-convex problem.

To avoid the inter-penetration of the analytical functions defining the contours of organs during deformable registration, we introduce a collision term to penalize the movement of DOFs (control points). The second energy term is computed by using signed distance maps $\mathcal{D}: \mathbf{X} \in \mathbb{R}^2 \to \mathcal{D}(\mathbf{X}) \in \mathbb{R}$. Each contour retains its own distance map describing how far a position in the image domain is from the contour: $\mathcal{D}_{bl}$ for the bladder, $\mathcal{D}_{va}$ for the vagina and $\mathcal{D}_{re}$ for the rectum. The
collision energy term is defined by a quadratic function with regard to the distance $D, \mathcal{F}: D \in \mathbb{R} \to \mathbb{R}^+$:

$$ 
\begin{aligned}
E_{Bl} &= \int_{\mathbf{X} \in \Omega_{Bl}(\text{Band})} \mathcal{F}(\mathcal{D}_{va}(\mathbf{X})) d\Omega_{Bl} \\
E_{Va} &= \int_{\mathbf{X} \in \Omega_{Va}(\text{Band})} [\mathcal{F}(\mathcal{D}_{bl}(\mathbf{X})) + \mathcal{F}(\mathcal{D}_{re}(\mathbf{X}))] d\Omega_{Va} \\
E_{Re} &= \int_{\mathbf{X} \in \Omega_{Re}(\text{Band})} \mathcal{F}(\mathcal{D}_{va}(\mathbf{X})) d\Omega_{Re}
\end{aligned} 
\quad (2.18) $$

where $\lambda$ defines the width of the possible intersection zone associated with the contour, so that it penalizes the possible intersections when one contour happen to approach the narrow band of another closely. As shown in Figure 2.20, the vagina is at the middle of the three organs, so the bladder and the rectum are never in contact. We use the following equations to compute the collision term for each contour:

$$ 
\begin{aligned}
E_{Bl} &= \int_{\mathbf{X} \in \Omega_{B/(Band)}} \mathcal{F}(\mathcal{D}_{va}(\mathbf{X})) d\Omega_{Bl} \\
E_{Va} &= \int_{\mathbf{X} \in \Omega_{Va}(Band)} [\mathcal{F}(\mathcal{D}_{bl}(\mathbf{X})) + \mathcal{F}(\mathcal{D}_{re}(\mathbf{X}))] d\Omega_{Va} \\
E_{Re} &= \int_{\mathbf{X} \in \Omega_{Re}(Band)} \mathcal{F}(\mathcal{D}_{va}(\mathbf{X})) d\Omega_{Re}
\end{aligned} 
\quad (2.19) $$

As the virtual image, $\mathcal{F}$ is used as a level-set function for computing the energy associated distance map. The map for each organ is redefined once at each iteration. Thus, for example, the value of $\mathcal{F}(\mathcal{D}_{va}(\mathbf{X}))$ for a given position $\mathbf{X}$ is obtained directly by projecting $\mathbf{X}$ onto the map of the vagina. We do not recompute the distance $\mathcal{D}_{va}(\mathbf{X})$ for each discrete point at the same iteration. Finally, these two energy terms are summed to form a global energy term as the cost function.

$$ E_{global} = E_{image} + \alpha E_{collision} \quad (2.20) $$

$E_{image}$ is to be associated with each organ. $E_{collision}$ is to be replaced by $E_{bl}, E_{va}$ and $E_{re}$, respectively. The coefficient $\alpha$ is chosen empirically during tests; the same value is used for the images of the same spatial resolution.

**Optimization** In our registration procedure, optimization is the step connecting the cost function (energy) and the model. It updates iteratively the DOFs of the model. For such a non-linear optimization problem, the algorithm is based on gradient. We opt for a gradient descent algorithm in this paper, which is suitable for optimizing the position and the number of the control points, and for controlling the convergence with regards to the adaptive refinement step. However, the Newton-Raphson method can also be computed similarly, The searched direction of the DOFs depends on the
derivatives of the cost function. Let Ti be one of the DOFs in Equation 2.12, the derivative is computed as follows:

$$ \frac{\partial E_{\text{image}}}{\partial T_i} = \int_{\mathbf{X} \in \Omega(\text{Band})} -2 [I_r(\mathbf{X}) - I_v(\mathbf{X})] \cdot \underbrace{\frac{\partial I_v}{\partial \mathbf{X}}}_{A} \cdot \underbrace{\frac{\partial \mathbf{X}}{\partial T_i}}_{B} d\Omega \quad (2.21) $$

As mentioned in the B-spline Model section, our model is optimized via a multi-scale registration procedure. At the initialization step, the user defines a point inside of each organ as the center of B-spline curve. We have tested a new initialization, while asking the user to put the point approximatively at the zone where the portion of organ is larger. We have observed a better increase of the initial contours at the first steps. Concerning the number of the control points, thanks to the adaptive refinement that we consider as one of our major contributions, we think it may not have a real importance because the points are to be added and deleted automatically after. And we give 10 by default. Then, a circle-like B-spline with 10 control points is created for each organ, which is a basic and simple one-click initialization (see figure 2.21).

For the coarse registration, the parameters of affine transformation (Equation (3)) are optimised (search direction computed by Equation (13)). In Figure 7B, one can observe that after the affine transformation, the three circle-like shapes are transformed to ellipse-like ones, which give an approximate initialisation of the detection procedure. Following the first step, a finer registration is applied using deformations based on B-spline model. Thus, the positions of control points are DOFs to be updated to fit the organ contours by means of minimising a cost function:

$$ \frac{\partial E_{\text{image}}}{\partial \mathbf{P}_i} = \int_{\mathbf{X} \in \Omega(\text{Band})} -2 [I_r(\mathbf{X}) - I_v(\mathbf{X})] \cdot \underbrace{\frac{\partial I_v}{\partial \mathbf{X}}}_{A} \cdot \underbrace{\frac{\partial \mathbf{X}}{\partial \mathbf{P}_i}}_{B} d\Omega \quad (2.22) $$

The partial derivative A is the gradient of virtual image that can be computed by the derivative of the level-set function. The partial derivative B is the mapping function, which means B-spline interpolation function in this case. Thanks to the power of analytical formulation, the required derivatives are simple to compute. In addition, the chain rule of derivation shows the independence of geometric model and the image energy computed by the difference between the virtual image and the real one. Hence, one can use the virtual image approach with different geometric representations (or mapping function). The collision term $E_{\text{collision}}$ can be computed in a similar way. Because of the nature of gradient descent, after certain iterations,
<figure><img src="image_57.png" /><figcaption>Fig. 2.21: Contour detection with the cost functions. (A) Three small circles initialised by user with 10 control points for each. (B) Positions after an affine transformation. (C) Final contours with optimal number of control points. (D–F) Illustration of the cost function metrics $E_{\text{image}}$ along the organ contours: $[I_r(u, v) - I_V(u, v)]^2$, with u evolves in the tangent direction of the contour and v in the direction perpendicular to the contour.</figcaption></figure>
it shows a phenomenon of local minimum (the energy does not evolve). We use this phenomenon as the criterion for adding control points. A test of local minimum is performed automatically for every 10 iterations. Let $E_{max}$ be the maximum value of the cost function during this 10 iterations, $E_{min}$ the minimum value and $E_0$ the initial value $\frac{\|E_{max}-E_{min}\|}{E_0} < \varepsilon$. If the criterion is satisfied, the optimization is supposed to reach a local minimum (low variation of the cost function); then, the curve is refined where the cost function needs to be reduced. The refinement step is treated automatically by inserting and removing control points based on the evaluation of cost function along the contour. Once the curve is refined, the optimization passes the local stability and continues to find the best correlation. This test and the adaptive refinement are performed repeatedly until the refinement can no more affect the cost function. By integrating insertion and removal of control points, our method can refine adaptively the geometric model with respect to images. The major advantage is that our B-splines fit properly the forms even in case of local aberrant curvature changes and avoid useless control points for description of the curve. This feature makes our method capable of a balance between model complexity and simplicity of optimization. The insert of control points does not change the B-spline curve by using its geometric algorithms. Every time when the gradient descent algorithm reaches its local minimum detected by a threshold of variation of the cost function, we refine the curve by subdividing the knot interval where there appears the maximum value of cost function. On the other hand, when two control points approach each other within a distance threshold, one is deleted for eliminating inter-cross of curve. In Figure 2.21F, we illustrate the distribution of cost function along the contour at the three stages (initialization, after affine transformation and the final position). The values of the cost function decrease significantly during the optimization. Moreover, the distance from a position on the contour to the organ border can be indicated by evaluating the metric map along contour. By zooming inside the contour, higher values occur when it is far from the real image border, which can be the indication for refinement.

## 2.3.2 Results

We validate our method of segmentation on a database of 19 patient images. The types of these patients include the witness, the prolapse and the endometriosis (a pathology of the pelvic system related to the hypo-mobility). Each patient data set consisted of a T2-weighted MRI of the pelvic system. We used for this 2D segmentation test the image in the mid-line sagittal plane of each patient, which is representative and useful for further study of the pelvic mobility. By the end of each
optimization process, a medical doctor manually corrected the shape of organs by modifying the positions of control points so that the final curves and the contours of organs, to the medical doctor's opinion, were well fitted. These final curves were used as the reference in our validation. In Figure 2.22, we illustrated the detection and the validation for each type. Then, numerically, we compared the curves before manual corrections with the final ones by evaluating error with several metrics: Dice coefficient, Hausdorff distance and average curve distance (ACD). Figure 2.23 shows the computation of these metrics by using the right contour (green) and an arbitrary contour (red, and manually created to show a larger error). In this paper, suppose that F is the set of edge points on the semi-automatically segmented region and G is one of the manually segmented region. The ACD values represent the average (L1-norm) of distances in millimeter between the edge points on the manually segmented organ and the automatically detected organ. The minimum distance for the ith point in F to the set G is $d_i^{FG}$, thus the ACD values can be computed as in Equation 2.23. The number $N_F$ and $N_G$ were chosen to make sure that approximately an edge point was used by pixel.

$$ 
\begin{aligned}
ACD^{FG} &= \frac{1}{N_F} \sum_{i=1}^{N_F} d_i^{FG} \\
ACD^{GF} &= \frac{1}{N_G} \sum_{i=1}^{N_G} d_i^{GF} \\
ACD &= (ACD^{FG} + ACD^{GF}) / 2
\end{aligned} 
\qquad (2.23) $$

In addition, for the Dice coefficient values, a higher Dice value assesses a better segmentation quality. On the 19 patients data, we got a similarity of over 90% for bladder and 80-90% for vagina and rectum especially because the ends of these two organs are difficult to define. In Figure 2.24, we illustrated some exceptions of the detection due to the unwanted effect presented in the images. Generally, the bladder is well highlighted under acquisition, while as the injected gel can exit the cavity with the patient is pushing, blurs and incertitude will occur in the image for the vagina and the rectum (Figure 2.24B). Because of the complexity of the shape, crease will occur in the 2D plane of the rectum (Figure 2.24A) and intrinsic material generates spots in the rectum in the image (Figure 2.24C). All these unwanted effects would make the algorithm limited in certain cases. In these cases, the algorithm detected contours. However, the contours are not the anatomical contours of the organs. The quantitative results of our tests are shown in Table ??.