This work is the result of three intense years within which I have evolved and grown as a scientist and as a person. The first year was extremely difficult in terms of motivation and also uncertain. But thanks to the support of the persons surrounding me, the realisation of this work has been possible. I would first like to thank my advisor Stéphane Cotin without who I would still be lost somewhere in the Amazonian forest in Peru. Thank you for your humanity and understanding my concerns. Thank you for your time and dedication, and for sharing and spreading your passion for research! I would like to express my gratitude to Elias Cueto, to José David Martín-Guerrero, to Yohan Payan and to Souvik Chakraborty for accepting being members of the jury of my PhD defense. I would like to particularly thank José David Martín-Guerrero's team and María José Ruperez Moreno for the inspiring scientific conversations and for being a reference starting point of my PhD.

I would also like to thank Raphael Sznitman and the ARTORG center in Bern, for the support in publishing my first PhD paper. Thank you Pablo Marquez-Neila for your precious help in the technical aspects of deep learning. Thank you Tatiana Fountoukidou and Jan Hermann for the patience when running the experiments with the OCT and the robot arm.

I feel grateful for the amazing time spent in the Mimesis team: you guys are awesome! Thank you for the inspiring conversations and for the sacred breaks playing Quoinche! I would like to particularly thank my deskmate Jean-Nicolas Brunet for all the precious advice, the patience when I did not understand, for being a real friend and listening when I needed to be heard. Thank you Antoine Petit and Sergei Nikolaev for not letting me alone at Clovis during the paper deadlines. Thank you Alban Odot for turning my research scripts into a usable optimized plugin and also for your bad jokes ;) Thank you Hugo Talbot and Rémi Bessard for the rooftop garden! And big thank you to all the team members for the smiles, the croissants, the expertise! I also would like to thank very important persons in my life whose support was essential for the realization of this work. Gracias Aaron por haber estado a mi lado ii a lo largo de este viaje, de principio a fín. Gracias Amoñi por estar cerca mío, ahora y siempre. Merci Gérard de m'avoir rappelé que je suis humaine et pas une machine. Sans toi, j'aurais brûlé toutes mes cartouches à mi-chemin. Thank you Lou for showing me the importance of listening to my body and to my intuition. Thank you Nibiru and Nébula for supporting my mental health, in particular during the lock-down. Thank you to all of my friends that loved me despite my bad humor and unavailability. In particular thank you Adèle and Gaby for being there. Gracias ama, aita y Mario, por creer y confiar en mí, y por haberme dado el tiempo y el espacio necesarios.

I would also like to thank my self for all the effort, the confidence and perseverance dedicated to this work.

iii

Notations and abbreviations

General notations

• < x, y >: dot product of x and y • : strain tensor

• c : right Cauchy-Green strain tensor

• g : Green-Lagrange strain tensor

• n: surface normal

• t n : stress vector linked to normal n

• dS: infinitesimal surface

• dV : infinitesimal volume

• V e : volume of element e

• σ: stress tensor

• σ c : Cauchy stress tensor

• σ p : first Piola-Kirchhoff stress tensor

• σ s : second Piola-Kirchhoff stress tensor

• C: constitutive matrix (stress-strain matrix)

• W ( ): strain energy function

• λ and µ: Lamé coefficients

• ρ: density of the material

• f s : surface forces

• f v : volume forces

• E: Young's modulus

• ν: Poisson's ratio

• N i : shape function at node i

• K: global stiffness matrix

• K: tangent stiffness matrix in the Newton-Raphson solver vii

• r: internal elastic force vector

• N e : number of elements in the finite-element mesh Chapter 3

• f : neural network function

• x = (x 1 , ..., x d ): input of f
• y: output of f

• θ: parameters of the network

• h j : artificial neuron j

• w j = (w j,1 , ..., w j,d ): weights linked to neuron j

• b j : bias of neuron j

• φ: activation function

• L: number of hidden layers

• k: layer index

• z (k) (x): output of a neuron in layer k before activation

• a (k) (x): output of a neuron in layer k after activation

• W (k) : weights matrix at layer k (number of rows equal to the number of neurons in layer k and number of columns equal to number of neurons in layer k -1

• b (k) : weights vector at layer k • h: network function

• n x × n y × n z : resolution of the grid

• C: input constraints (expressed in the grid nodes) viii

• U s : input surface displacement (expressed at grid nodes)

• U v : output volumetric displacement field (expressed at grid nodes)

• N : number of samples in training data set

• M : number of samples in testing data set

• n: number of degrees of freedom of the FE mesh (e.g. number of nodes in sparse grid)

• c: number of channels in the first layer of the U-Net

• k: number of steps in the U-Net

• Γ N : surface area on which forces are applied (its location varies such that the boundary of the computation domain is completely covered)

• Λ: number of samples per Γ N

• e: mean norm error of a sample (mean euclidian error in Chapter 8)

• ē: average error over the testing data set

• σ(e): standard deviation of the error over the testing data set Chapter 5

• p: sought parameter

• µ 0 and σ 0 : initial mean value and standard deviation of p

• n: number of nodes

• k: number of stiffness parameters

• m: number of elements in the stochastic state vector Chapter 6

• f p : external force due to intraocular pressure

• f n : needle pushing on the surface of the sclera

• f g : gravitational force

• S: surface area of the eye • CNN: convolutional neural network

• CT: computerized tomography

• FC: fully-connected

• FE: finite-element

• FEM: finite-element method

• FSS: feature space size

• FUS: freehand ultrasound system

• GPU: graphics processing unit

• H8: 8-node hexahedral elements

• HPOD: hyperreduced proper orthogonal decomposition

• IBM: immersed-boundary method

• ICP: iterative closest point

• IOP: intraocular pressure

• MIS: minimally-invasive surgery

• ML: machine learning

• MLP: multi-layer perceptron

• MRI: magnetic resonance imaging

General Introduction

There are many applications in engineering where the deformation of nonlinear structures needs to be simulated in real time, or would benefit from being computed interactively. Some important examples can be found in the field of medicine, in order to develop training systems for learning surgical skills [Ayache et al., 2006] or in the field of surgical navigation, where augmented reality combined with interactive simulations can bring significant improvements to clinical practice [Haouchine et al., 2013a]. Medical robotics, involving flexible robots or interactions with soft tissues, is another important area where real-time simulation of flexible structures is essential, in order to have a better control of the robot [Duriez, 2013].

In the context of minimally-invasive surgery (MIS) or laparoscopic surgery for instance, extensive surgical training is mandatory to support the success of the intervention [START_REF] Rosser | Minimally invasive surgical training solutions for the twenty-first century[END_REF]. During this type of surgery, usually three to four small incisions are done in the abdomen, through which surgical tools and a camera are inserted in the abdominal cavity. The abdomen is then inflated with gas, thus creating a working space for the surgeons to operate guided by the images from the camera, displayed on a screen. MIS offers several benefits to the patient such as reducing the bleeding and the risk of infection as well as shortening recovery time [START_REF] Li | Laparoscopic vs open hepatectomy for hepatolithiasis: An updated systematic review and meta-analysis[END_REF]. However, the reduced field of view of the endoscopic camera and the lack of tactile information, make it difficult to succeed in such interventions [START_REF] Plantefeve | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF]. Surgeons need a lot of dexterity so an extensive training is required. Instead of practicing on animals or on cadavers, virtual laparoscopic training allows surgeons to practice on a completely virtual environment. Through virtual reality and simulation, it is possible to reproduce with high fidelity what the surgeon experiments during the surgery. Apart from the deformation of the organs, the interactions between surgical tools and the tissues are also modeled using haptic devices to mimic tactile sensation. To obtain a stable haptic feedback, simulations must run at 500 frames per second (FPS) [START_REF] Courtecuisse | Haptic Renderingof Hyperelastic Models with Friction[END_REF]. Interventional radiology is another example of intervention requiring specific surgeon training due to its difficulty [START_REF] Messina | Integrated fellowship in vascular surgery and intervention radiology: a new paradigm in vascular training[END_REF]. Indeed, the anatomy is visualized through Xray images, and the surgical gesture is performed with flexible instruments (such as catheters) through the vascular tree towards the target. In such training systems, the surgical tools are deformable objects and the vessels are generally considered to be rigid. Such models need to be interactive, therefore the simulations must be computed in real-time (e.g. 30 FPS for visual interaction [Duriez, 2013]).

Computer-assisted interventions (CAI) also require the computation of nonlinear deformations in real-time. For example the use of surgical robots is revolutionizing the field of MIS as they can overcome several human limitations [Troccaz, 2012].

The control loop of the robot is often based on intraoperative images and it must take into account the interactions of the robot with soft tissues, as well as its own deformation in the case of soft robotics. In needle insertion procedures for instance, both the needle and the tissue usually deform as the insertion proceeds and authors in [START_REF] Adagolodjo | Inverse real-time Finite Element simulation for robotic control of flexible needle insertion in deformable tissues[END_REF] proposed to control the robot through real-time numerical simulations.

On another hand, in the context of hepatic surgery, the objective is to accurately locate the internal structures such as tumors and blood vessels (that need to be preserved for the post-operative regeneration of the liver tissue). While the initial position of these structures is known from the preoperative images, their actual position during surgery is often hidden or uncertain (even for well-trained and experienced surgeons). To help the surgeon overcome these difficulties, augmented reality (AR) techniques are used to enrich visual information through a fusion of intraoperative images and a preoperative 3D model of the patient's anatomy. However, correctly aligning preoperative data to intraoperative images remains an unsolved problem, especially when large deformations are involved and only sparse input data is available. This is typically the case when trying to provide an augmented view of an organ during surgery. In most cases, just about 30% of the surface of the organ is visible due to the limited field of view of the laparoscopic camera or size of the incision [START_REF] Plantefeve | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF].

In all the cited examples, the considered structures are deformable and can be very soft. Surgical manipulations cause the organs to deform and interact with the surrounding anatomy. Several works have demonstrated the benefits of physics-based models, particularly patient-specific biomechanical models, for accurate registration between different preoperative 3D anatomical model and intraoperative data [START_REF] Clements | Robust surface registration using salient anatomical features for image-guided liver surgery: Algorithm and validation[END_REF], Haouchine et al., 2013b, Suwelack et al., 2014, Alvarez et al., 2018]. While there are different numerical strategies for solving these models, we only consider the finite-element (FE) method throughout this manuscript, for its accuracy and ability to simulate a large range of materials on potentially complex domains. It gives a numerical approximation of a partial derivatives equation discretizing the considered object using nodes connected by elements. The choice of the constitutive model, and its parameterization, are obviously key to an accurate registration. It is usually acknowledged that a registration of internal structures below 5 mm is needed for best clinical impact, such as targeting relatively small tumors [START_REF] Ruiter | Model-based registration of X-ray mammograms and MR images of the female breast[END_REF]. Achieving the demanded accuracy within real-time constraints raises several challenges that are far from being solved.

When considering augmented reality or, more generally, real-time tracking of the organ deformation to provide up-to-date guidance, the computational efficiency of the FE method is essential. Indeed, if the simulations are not up-to-date, they might simply be ignored by the clinicians. In the case of augmented hepatic surgery for instance, intraoperative images are acquired at about 20 Hz leading to update times of less than 50 ms. During this small amount of time, acquisition and processing of the images as well as model update need to take place. As a result, FE computation times should require less than 30 ms. If only small deformations take place, achieving such computation times is feasible [START_REF] Meier | Real-time deformable models for surgery simulation: a survey[END_REF]. However, if large nonlinear deformations happen, computation times become incompatible with such constraints. A solution might be the use of the co-rotational FE method, where geometrical nonlinearities can be handled in real-time [Haouchine et al., 2013a, Petit et al., 2018]. Nevertheless, when more complex biomechanical models need to be used these optimizations no longer hold. Alternative solutions have been proposed with different trade-offs regarding the ratio between computation time and model accuracy [START_REF] Peterlík | Modeling and real-time simulation of a vascularized liver tissue[END_REF], Suwelack et al., 2014, Modrzejewski et al., 2018, Niroomandi et al., 2008, Johnsen et al., 2015, Allard et al., 2007].

[ [START_REF] Marchessau | [END_REF] proposed the Multiplicative Jacobian Energy Decomposition (MJED) that allows for fast and realistic liver deformations including hyperelasticity, porosity and viscosity. Also, [START_REF] Miller | Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF] introduced Total Lagrangian explicit dynamics (TLED) which can achieve real-time performances when coupled with explicit time integration and GPU-based solvers [START_REF] Joldes | Real-time nonlinear finite element computations on GPU-Application to neurosurgical simulation[END_REF].

In recent years, machine learning (ML) started to revolutionize several fields (vision, language processing, image recognition, genomics). With sufficient ground truth data, machine learning algorithms can map the input of a function to its output without any mathematical formulation of the problem. The high inferring speed of these methods can be useful for many applications where the prediction speed is of critical importance. Due to this characteristic and the fact that they are driven directly by data, these methods seem promising for the learning of the entire mechanical behavior of the anatomy without relying on prior models. Some first attempts that exploit learning methods to estimate the deformation of biological tissues have recently been made. By implicitly encoding soft tissue mechanical behavior in the trained ML models, they proved successful to predict the entire 3D organ deformation starting either by applied surface forces [START_REF] Morooka | Real-time nonlinear FEM with neural network for simulating soft organ model deformation[END_REF], Tonutti et al., 2017, Rechowicz et al., 2013] or by surface displacements [START_REF] Lorente | A framework for modelling the biomechanical behaviour of the human liver during breathing in real time using machine learning[END_REF]. The accuracy of a ML model highly depends on the network architecture and on the quality and on the amount of data used to train it. In an ideal scenario, such a model would be trained with an infinite amount of real patient-specific noise-free data, which is in practice not possible. Within this framework, FE simulations can be exploited to generate synthetic data that is highly representative of the reality, to be used as training samples.

Among the various ML techniques, the use of neural networks (NN) has considerably increased. This is due to the fact that they are the building blocks of deep learning, a class of methods which is able to learn data representations and has demonstrated strong abilities at extracting high-level representations of complex processes. For example, NN are used by [START_REF] Tonutti | A machine learning approach for real-time modelling of tissue deformation in imageguided neurosurgery[END_REF] and [START_REF] Rechowicz | Development and validation methodology of the Nuss procedure surgical planner[END_REF] to predict the displacement of brain tumors and of the rib cage surface respectively, starting from the acting forces. However, both these works do not predict whole volume deformation but only surface displacements. NN based methods have been also used to predict liver deformation in augmented surgery. [START_REF] Morooka | Real-time nonlinear FEM with neural network for simulating soft organ model deformation[END_REF] trained a NN to predict liver deformations for a given input force. They use their model together with principal-component analysis to compress the size of the output deformation modes, and thus reduce the training time. Although the model proved able to learn the deformation modes, it was only applied to simulated data and not to real cases. From all these works it emerges that the main advantage of using neural networks to predict anatomical deformations is that the prediction speed is in the order of few milliseconds and is not affected by the complexity of the model used to generate the training dataset.

To ensure the aforementioned requirements in terms of model parameterization and computational efficiency, we propose to combine patient-specific FE simulations with deep learning techniques. Within this thesis work, our main contribution is the use of deep neural networks for real-time numerical simulations of nonlinear deformations. This manuscript is divided in two main parts. The first part is dedicated to generic real-time simulations with deep neural networks. The first two chapters establish the theoretical foundations of the problems we are looking at. We first present the mechanical formulation of our problem and the finite-element method used for its numerical resolution. In a next chapter, the fundamentals of artificial neural networks are introduced and we give a brief explanation on how such networks are trained. In the third chapter, we present the main contribution of our work: the U-Mesh framework. It consists of a data-driven deep neural network that learns the desired biomechanical model to predict complex nonlinear deformations in real-time on simple geometries.

The second part of this thesis deals with adapting the finite-element method and the deep learning algorithms to patient-specific modeling. We look at the importance of correct model parameterization and we describe a stochastic assimilation method to identify the patient-specific parameters. Chapter 4 describes how to perform patientspecific modeling using Kalman filtering. Chapter 6 concerns medical robotics. A neural network for robotic force estimation is trained with patient-specific finiteelement simulations. The objective is to estimate forces during robotized intravitreal injections using a neural network for image classification. Chapters 7 and 8 extend the U-Mesh framework to patient-specific geometries. We demonstrate its potentiality in computer-assisted interventions, in particular augmented liver surgery and US-guided breast biopsy. In the last chapter of this manuscript, we propose several strategies in order to adapt the network's prediction to patient-specific properties.

Part I

Real-time numerical simulations 8

Chapter 2 Biomechanical Problem Formulation

Elasticity equations for soft tissue simulation

A biomechanical model is a mean to represent the physical behavior of an organ whose mechanical properties and boundary conditions are difficult to determine. In order to construct such a model, we explain the physical phenomena considered in this work using the theory of continuum mechanics (assumption that the studied object fills entirely the space it occupies). The goal of a simulation is to predict the behavior of a deformable body from the knowledge of the external forces f ext applied to the body. There are two types of forces undergone by an object: the external forces (such as gravity and contacts) and the internal forces that are given by the stress tensor. The internal forces maintain the coherence of the body: a deformable body cannot be infinitely stretched or compressed. Applying forces on a deformable body leads to the creation of constraints (stress) and deformations (strain) inside of it. Such quantities will be defined in next sections. From Newton's second law of motion we have that:

f ext + f int = ma, (2.1)
where f int are the internal forces, m the body mass and a the acceleration. In the static case we can rewrite this equation as:

f ext = -f int .
(2.2)

A measure of deformation: the strain tensors

We observe a continuum solid occupying a volume Ω whose boundary is Γ. We consider a Lagrangian description of the movement meaning that the deformed state of a particle p ∈ Ω at time t is given by p

= φ(p, t), p = φ(p, 0) (2.3)
where p is the initial position of particle p and φ is the transformation undergone by the solid with φ : Ω × R + → Ω.

A naive approach to quantify the deformation between two infinitesimally close points p1 and p2 is to measure their respective distances in the initial and deformed configurations. In the initial configuration, such distance reads as dp = p2 -p1 (2.4) and in the deformed configuration:

dp = p 2 -p 1 = φ( p1 + dp, t) -φ( p1 , t) (2.5)
where p 1 and p 2 are two infinitesimally close points in the deformed configuration and dp the distance between them. For each particle in the rest configuration, there exists a displacement vector u, u : Ω × R + → Ω, describing its deformed position such that u(p, t) = φ(p, t)p.

(2.6)

The displacement field u characterizes the movement of the solid. Nevertheless note that for a rigid deformation, u is different from zero even if the object is not deformed. Hence, a better definition of the deformation must be introduced.

Let us introduce the deformation gradient in R 3×3 :

F = ∇u + I, (2.7)
where I ∈ R 3×3 is the identity matrix and ∇u is the gradient of u. F describes the variation of the distance between the particles in all the directions. Equation (2.5) can be rewritten as dp = Fdp.

(2.8)

The square of distance dp leads to (dp

) 2 = dp • dp = (Fdp) • (Fdp) = (Fdp) T (Fdp) = dp T F T Fdp = dp T c dp,
where c = F T F is known as the right Cauchy-Green strain tensor. However, when the difference in angle between the rest and deformed configuration is null, c is not null. Indeed, if dp = dp, c = I. Therefore, we compute the difference between the squared distances (dp

) 2 -(dp) 2 = dp T c dp -(dp • dp) = dp T ( c -I)dp = 2dp T g dp,
where g = 1 2 ( c -I) is known as the Green-Lagrange strain tensor satisfying g = 0 when dp = dp. It can be rewritten as:

g = 1 2 (F T F -I) = 1 2 (∇u + I) T (∇u + I) -I = 1 2 ∇u T ∇u + 1 2 ∇u T + ∇u .
The Green-Lagrange tensor consists of a linear component linear = 1 2 ∇u T + ∇u and a nonlinear one nonlinear = 1 2 ∇u T ∇u . By linearization of g , for small deformations nonlinear can be neglected and g can be approximated by:

:= linear = 1 2 ∇u T + ∇u .
(2.9)

The infinitesimal strain tensor considerably simplifies the equations. However, as soon as the deformations get large, this linearization is no longer valid and the complete g must be considered.

A measure of internal forces: the stress tensor

We now need a way to account for the internal forces of an object that constrain its deformation (avoiding infinite stretch or compression for instance). Let t n be the stress vector at point p whose normal to the surface is n:

t n = f dS , (2.10)
where f is the average force applied on the infinitesimal surface dS. The stress vector t n depends on the normal n, hence it does not completely describe the impact of the force in the whole volume. If we consider an infinitesimal volume dV , for example a tetrahedron (see Figure 2.1), the force at the wider face whose normal is n is expressed as: where dS is the area of the considered face. On the other faces, the force is expressed as:

f = t n dS, (2.11) Biomechanical Problem Formulation 12
f i = -t i dS i i ∈ {1, 2, 3}, (2.12)
where t i is the stress vector and dS i is the area of the i-th face of the tetrahedron. The minus sign comes from the normals that are equal to the negative of the coordinate axis. Let f v dV be the volume force inside the tetrahedron. From the conservation of linear momentum we have that:

f + f 1 + f 2 + f 3 + f v dV = ρdV dv dt , (2.13)
where ρ is the mass density, v is the velocity of the tetrahedron. Using equations (2.11), (2.12) and expliciting dV , Equation(2.13) can be rewritten as:

t n dS -t 1 n 1 dS -t 2 n 2 dS -t 3 n 3 dS + f v hcdS = ρhcdS dv dt , (2.14)
where h is the height of the tetrahedron, c is a constant independent of h and (n 1 , n 2 , n 3 ) are the coordinates of the normal n. If we devide by dS and make h tend to 0 we get that:

t n -t 1 n 1 -t 2 n 2 -t 3 n 3 = 0. (2.15)
In other words, if the stress vectors acting on three faces of normal parallel to the coordinate axes is known, we can compute the stress vector for any normal direction n:

t n = t 1 t 2 t 3     n 1 n 2 n 3     =     σ 11 σ 12 σ 13 σ 21 σ 22 σ 23 σ 31 σ 32 σ 33         n 1 n 2 n 3     = σ c • n (2.16) Biomechanical Problem Formulation 13
where σ c is called the Cauchy stress tensor.

There exist other definitions for the stress tensor such as the first Piola-Kirchhoff tensor σ p : (2.17) where J = det(F) is the matrix allowing to pass from the rest to the deformed configuration. However, σ p is not symmetric and is based on the deformed configuration.

σ p = J • σ c • F -T ,
On the contrary, the second Piola-Kirchhoff tensor defines the stress based on the rest configuration:

σ s = JF -1 • σ c • F -T .
(2.18)

Constitutive law: linking stress and strain

There is a relation between the strain tensor and the stress tensor σ that depends on the mechanical properties of the material. Such function translates the link between the forces applied on the considered body and its deformation. Even if it can be very complex, for small displacements we can assume its linearization using Hooke's constitutive law. where C is the constitutive matrix (see [Slawinski, 2007] for detailed equations). If the material studied is isotropic, such equation can be written as: σ = 2µ + λtr( )I.

(2.20)

where λ and µ are the so-called Lamé coefficients, which in the present case are related to the material properties of the soft body (introduced later in Equation(2.41)).

Note that Hooke's law is linear in but also in ∇u because we used the linearized strain tensor. If we do not work with small displacements, we need to take account for nonlinearities of the displacement and of the material, and the relation between strain and stress is no longer linear in ∇u (cf Figure 2.2). The material should therefore be modeled using an hyperelastic constitutive law such as the Saint Venant-Kirchhoff, Ogden, Neo-Hookean or Mooney-Rivlin models.

Hyperelasticity: the Saint Venant-Kirchhoff model If we rewrite the Equation (2.20) with g instead of the linearized , we obtain the Saint Venant-Kirchhoff model where σ is not linear in ∇u anymore. It is the simplest and most efficient extension of a linear elastic material to the nonlinear regime. In the following chapters, we will most of the time use the Saint Venant-Kirchhoff model.

It goes without saying that the constitutive law can be made more complex in order to model a higher degree of nonlinearity but the computational needs also increase. For a given mesh resolution and configuration (constant material properties and boundary conditions) the computation time needed for solving linear elasticity equations can be much smaller than that needed for hyperelasticity equations. Therefore, it is crucial to adapt the choice of the constitutive model to the considered deformations (to avoid using unnecessarily complex and computationally expensive models).

Forces and weak formulation

In the previous sections, we have defined the strain tensor linking the deformation and the displacement. Then, we have described the stress tensor translating the internal forces undergone by the soft body. Finally we have characterized the relation function between the strain and the stress tensors. To balance the internal forces with the external forces, we can express the force equilibrium over the computational domain Ω such that:

div(σ) + f ext = ρü, (2.21)
where ρ is the density of the material and f ext the external forces.

(x). n(x) = f s , ∀x ∈ Γ N .
(2.24)

The boundary conditions will be further detailed in Section 2.3.2.

In order to distinguish the forces that are applied at the boundaries of the body such as pressure and contacts (denoted f s ) from the volume forces such as gravity (denoted f v ), we introduce the weak formulation of Equation (2.22). The weak form consists of multiplying the strong form by a test function v (verifying v = 0 on Γ D ) and integrating over the domain Ω. It is given by Equation (2.28) using Green's formula:

-

Ω div(σ).v dΩ = Ω f ext .v dΩ (2.25) ⇔ Ω σ : (v) dΩ - Γ σ. n.v dΓ = Ω f ext .v dΩ (2.26) ⇔ Ω σ : (v) dΩ - Γ N f s .v dΓ - Γ D σ. n.v dΓ = Ω f v .v dΩ (2.27) ⇔ Ω σ : (v) dΩ Material stiffness - Γ N f s .v dΓ Boundary conditions = Ω f v .v dΩ Volume forces
(2.28)

We remark that the spaces of the solution u and the test function v are different:

u = u D on Γ D , whereas v = 0 on Γ D .
In the following, for the sake of simplicity, we will set u D = 0.

These equations must be applied to all the points that form the body but this would take an infinite amount of time. Hence, we need a method that approximates the solution to this formulation: the finite-element method that will be explained in the next section.

The Finite-Element Method

The equations of the previous section must be integrated all over the domain Ω, which leads to an infinite number of equations. The finite-element method (FEM) is a numerical method that gives a numerical approximation of a partial derivatives equation discretizing the object using nodes connected by elements. Positions inside the elements are interpolated from the element nodes. Continuum mechanics produces a stiffness matrix for each element that relates nodal displacements to nodal forces. FEM has interesting properties for our application domain such as accuracy and robustness. However, its computation time and sensitivity to the mesh resolution raise several challenges in the context of real-time simulation of organs. In this section, we will present the FEM for the simulation of deformable objects.

Domain discretization

The first step is the domain discretization. It consists of generating a mesh composed of simple geometrical shapes (for example triangles or quadrilaterals in 2-D or tetrahedra or hexahedra in 3-D) that composes a non-gap and no intersection subdivision of the computational domain. We then obtain a finite number of points that are linked by elements. The main idea is to solve the equations at these elements and then approximate the solution over the entire domain Ω by a linear combination of element-wise continuous polynomial functions called the shape functions. The order of the polynomial gives the order of the element (usually linear elements or quadratic elements).

Tetrahedral finite-element meshes are usually preferred over hexahedral meshes as the latter are extremely difficult to generate automatically, and its generation takes several orders of magnitude longer than tetrahedral mesh generators [START_REF] Shepherd | Hexahedral mesh generation constraints[END_REF]. However, in this manuscript we will mainly use hexahedral finite element meshes. The main motivation for this choice is that we want to use the FE simulations to train Convolutional Neural Networks that work on regular grids (see Chapter 4 for details). Moreover, hexahedral meshes offer several benefits over tetrahedral finite element meshes. They present smaller error and less elements are usually needed to reach a given accuracy, and tetrahedral elements can exhibit a too stiff behavior [START_REF] Shepherd | Hexahedral mesh generation constraints[END_REF]. Moreover, trilinear hexahedral elements (Q 1 elements) provide similar accuracy and convergence characteristics as quadratic tetrahedral elements (P 2 elements) [START_REF] Cifuentes | A performance study of tetrahedral and hexahedral elements in 3-D finite element structural analysis[END_REF], Benzley et al., 1995]. For all these reasons, we will use trilinear hexahedral elements in our simulations.

To compute accurately the deformation of the considered object, the finite-element mesh needs to be sufficiently fine (to avoid discretization errors). The more we refine, the better the solution but at the cost of higher computational resources. Therefore, a compromise must be found in order to obtain the desired accuracy within an acceptable amount of time.

Shape functions

In this section, we will explain how to build the finite-element space in which we will search the numerical solution. For the sake of readability, we will only present the case of a linear constitutive law associated to the small displacement regime, since handling the linearized version is similar. We will now express the shape functions (also known as basis functions) that allow us to interpolate the value of the displacement field through each element based on the nodal values. Throughout this manuscript we only consider linear or trilinear shape functions e.g. first degree polynomials. Note that we could use higher degree polynomials but at the cost of higher computation times. For triangular or tetrahedral elements, the P 1 nodal shape function N i at node i reads as:

N i (x, y, z) = j+k+l≤1 α jkl x j y k z l , (2.29)
where α jkl are real coefficients defined such that N i = 1 at node i and N i = 0 at other nodes. For quadrilateral or hexahedral elements, the Q 1 shape function reads as: (2.30) where, again, α jkl are real coefficients defined such that N i = 1 at node i and N i = 0 at other nodes.

N i (x, y, z) = j,k,l≤1 α jkl x j y k z l ,
At any point p, the discretized displacement u h can be expressed as: (2.31) where u hi is the displacement at node i and N is the number of nodes in the mesh (excluding those on the Dirichlet boundary).There is one shape function N i for each node i in the mesh and such function must have a local support and be piece-wise continuous.

u h (p) = N i=1 N i (p).u hi ,
For readability reasons, we rewrite N i=1 N i as k N k in the following equations.

Biomechanical Problem Formulation 18

Build the system matrices

In order to write the system in a matrix form, we first need to write Equation (2.28) in its discrete form. The problem we aim at solving reads as:

∃?u h ∈ V h,0 such that Ω σ(u h ) : (v h ) dΩ - Γ N f s h .v h dΓ = Ω f v h .v h dΩ, ∀v h ∈ V h,0 .
(2.32) u h is the discretized version of the displacement u and is given by:

u h (p) = k N k (p).u hk .
(2.33)

f s h and f v h are the discretized forms of the surface and volume forces expressed with the shape functions. They read as follows:

f s h (p) = k N k (p).f s hk , (2.34) f v h (p) = k N k (p).f v hk .
(2.35) Space V h,0 is the discretization of the Hilbert space H 1 0 (Ω). If in addition, we apply the Galerkin theorem suggesting to replace the test function by the shape function N n (p), by linearity we get that:

k u hk Ω σ(N k (p)) : (N n (p)) dΩ - k f s hk Γ N N k (p).N n (p) dΓ = k f v hk Ω N k (p).N n (p) dΩ .
Let us define the stiffness matrix K:

K =     Ω σ(N 1 (p)) : (N 1 (p)) dΩ ... Ω σ(N 1 (p)) : (N N (p)) dΩ . . . . . . . . . Ω σ(N N (p)) : (N 1 (p)) dΩ ... Ω σ(N N (p)) : (N N (p)) dΩ     ,
(2.36) and the displacement vector U h :

U h =     u h1 . . . u hN     .
(2.37)

Let us also define the volume and surface force vectors as:

F s h =     f s h1 . . . f s hN     and F v h =     f v h1 . . . f v hN     .
(2.38)

Biomechanical Problem Formulation 19
Then, the full system reads as:

KU h = F s h + F v h .
(2.39)

To derive the stiffness matrix we made the assumption of small displacements. However, most of the times, the stiffness matrix cannot be expressed explicitly due to the material nonlinearity (relation linking the stress and the strain tensors) or to the geometrical nonlinearity (relation between the strain tensor and the displacement). Hence, it is some times impossible to separate the material stiffness from the displacement u. In such case, the weak formulation will be solved iteratively thanks to the Newton-Raphson algorithm (presented in next section) and K will be replaced by the associated tangent stiffness matrix.

Note that in the case of small displacements, we suppose that the matrix K does not change over time (meaning that we can compute its inverse once for all at the beginning of the simulation thus accelerating the computations considerably). For more complex materials, K will need to be recomputed at every simulation step. Intuitively, the matrix K explains the influence of a point on the other points of the mesh. It will be used to know the internal forces generated over the nodes by their displacement.

Solving the system of equations

Equation (2.39) represents the static system that aims at studying the equilibrium states of the deformation. The solution is the displacement for which the internal forces are equal to the external forces applied to the object. In this case, any linear system solver could be used for solving the equation. In the nonlinear case, we can obtain the solution of the system by iterating in the Newton-Raphson method for example. Basically we build a sequence in which we make the displacement vary until the difference of the forces vanishes. From an initial displacement u 0 , we try to find a correction δ n u after n iterations that balances the linearized set of equations:

Kn-1 (u n-1 )δ n u = r(u 0 + δ n-1 u ) + f (2.40)
where K is the tangent stiffness matrix and r is the internal elastic force vector. At each iteration, both the matrix K and the vector r need to be computed, and the linear system needs to be solved. Since the convergence of the Newton-Raphson method is only valid for a displacement u 0 near the solution, large external loads must be applied by small increments and can require a large number of iterations to converge.

If we were to consider dynamics, in Equation (2.21), ü would be different from 0. Dynamic systems aim at studying the deformations that depend on the time, thus taking into account the acceleration of the system. We want to know all the transitory states of the movement. In this case, a temporal integration scheme has to be used before the linear system solver. For that, we must define a time-step h and a time integration scheme in order to define the new positions in terms of the old ones. There are two possible strategies: explicit temporal integration or implicit integration. We will not go any further in the time integration since we only deal with static problems throughout this manuscript.

Model parameterization

Material properties

The proper modeling of an object requires the identification of the constitutive law (presented in Section 2.1.3) and the related material properties. These parameters play a major role in the estimation of the displacement field. Therefore, their understanding and identification is crucial for the accuracy of the method.

However, identifying such parameters can be very difficult in practice, mostly for soft tissues. Indeed, the typical approach for estimating them relies in traction experiments where different known forces are applied on a sample of the material.

The sample elongation and the change in its cross section are measured in order to estimate the Young's modulus E and the Poisson's ratio ν. The Young's modulus defines the ratio of stress to strain. It is a measure of the stiffness of the material. The Poisson's ratio characterizes the compressibility of the material (change in the volume). While the Young's modulus can take any positive value, the Poisson's ratio is positive and smaller than 0.5. Incompressible materials (such as soft tissues) have a Poisson's ratio close to 0.5. The Lamé coefficients λ and µ are theoretical material-dependent quantities that are in practice obtained form the values of E and ν through the following relations:

λ = Eν (1 + ν)(1 -2ν) , µ = E 2(1 + ν) . (2.41)
Note that such definition is only valid for linear elasticity and for some cases of hyperelasticity (such as the Saint-Venant Kirchhoff model that is the only hyperelastic law used in this manuscript).

For biological tissues, the traction experiments are controversial because the biological samples are examined ex vivo. Therefore, the mechanical properties obtained may differ from the in vivo properties. Other methods such as elastography [START_REF] Sarvazyan | Shear wave elasticityimaging: A new ultrasonic technology of medical diagnostics[END_REF], Muthupillai et al., 1995, Xu et al., 2007] or Bayesian filtering [START_REF] Peterlík | Image-driven stochastic identification of boundary conditions for predictive simulation[END_REF] attempt to determine these parameters in vivo. Whereas in the first part of this manuscript we will consider average material properties, in the second part we will identify them with Bayesian filtering (see Chapter 5).

Boundary conditions

In order to guarantee the uniqueness of the solution of Equation ( 2 [START_REF] Nikolaev | Stochastic Correction of Boundary Conditions during Liver Surgery[END_REF]. For instance, the structures attaching the organs (such as ligaments or muscles) are difficult to locate and characterize based on medical images. Therefore, in the second part of this manuscript, we will propose a way to characterize them using Bayesian filtering (similarly to material properties).

Conclusion

In this chapter we have presented elastic biomechanical models for soft-tissue simulation in order to estimate the displacement field u that is the quantity of interest. Moreover, the finite-element method is depicted with an emphasis on the discretization choices considering the further use of convolutional neural networks for displacement field estimation. In addition, we introduced the model parameters and boundary conditions that play a major role in the accuracy of the prediction as we will see in Chapter 5.

Biomechanical Problem Formulation 22

Biomechanical Problem Formulation 23

Chapter 3

Artificial Neural Networks: Fundamentals

Components of an artificial neural network

Definition

An artificial neural network (ANN) is a sequence of equations that approximates the underlying relationships in a given data set. Such a system learns the patterns in the data by analyzing examples, generally without receiving any task-specific rules. An ANN is an application f , nonlinear with respect to its parameters θ that associates an output y to an input x such that y = f (x, θ). A neural network (NN) can be used for solving regression or classification problems and it can be trained with labeled (supervised learning) or unlabeled data (unsupervised learning). The universal approximation theorem [Cybenko, 1989, Hornik, 1991] states that a feed forward neural network made of artificial neurons can approximate any real-valued continuous function on compact subsets of R n .

In Figure 3.1 a simple NN is depicted. The input layer appears in orange and consists of five inputs. It receives the information supposed to explain the function to be approximated. Then, there are two hidden layers with two and five artificial neurons respectively. A hidden layer is an intermediate layer allowing to model the nonlinear process. In such a feed forward NN, the outputs of each layer are the inputs to the next layers. The output layer appears in green and is the NN prediction.

Artificial neuron and activation functions

Each artificial neuron consists of a function h j of the input x = (x 1 , ..., x d ), weighted by a vector of connection weights w j = (w j,1 , ..., w j,d ), completed by a neuron bias b j , and associated to an activation function φ such that y j = h j (x) = φ(< w j , x > +b j ).

(3.1) • The identity function φ(x) = x,

• The sigmoid function

sigmoid(x) = 1 1 + e -x ,
• The hyperbolic tangent tanh(x) = e x -e -x e x + e -x ,

• The rectified linear unit (ReLU) ReLu(x) = max(0, x),

• The hard threshold function

φ β (x) = 1 x≥β .
Historically, the Sigmoid was the most used activation function as it is differentiable and translates the input ranged in [-∞; +∞] to the range [0; 1]. However, the exponential function is computationally expensive and the Sigmoid has the problem of vanishing gradients (the same problem arises for the hyperbolic tangent). The problem of vanishing gradients is inherent to the backpropagation optimization. As we move backwards in the NN, the gradients tend to get smaller and smaller leading to difficulties in training the first layers. In particular, training the model with vanishing gradients can be very time consuming and lead to low accuracy in the prediction.

On the other hand, the ReLU activation function does not suffer from vanishing gradients. It is cheap to compute and it accelerates the convergence of the Stochastic Gradient Descent algorithm (that will be explained later). However, as the ReLU function and its derivative are equal to 0 for negative values, no information can be obtained in that case for the concerned neuron. To overcome this issue, a small positive bias can be added to the activation function to ensure that each neuron is active:

φ(x) = max(x, 0) + min(x, 0), (3.2)
where is either a fixed parameter set to a small positive value, or a parameter to estimate (e.g. to learn).

ReLU is used only within hidden layers. For the output layer, a different activation function can be used depending on the type of problems we are dealing with (regression or classification). In binary classification, the output is a prediction of P(Y = 1/X) since it ranges in [0, 1] so the Sigmoid activation function can be considered. However, for multi-class classification, the output layer has one neuron per class i, giving a prediction of P(Y = i/X). The sum of all these values has to be equal to 1. Hence one can use the Softmax activation function that reads as:

softmax(z) i = e (z i ) j e (z j ) . (3.3)
On the other hand, for regression problems one can use as activation of the output layer the hyperbolic tangent, a linear function or the identity function (meaning no activation function).

Artificial Neural Networks: Fundamentals 26

Specific types of neural networks

In this section we will present the two types of neural networks used in this thesis: the multi-layer perceptron with fully connected layers and convolutional neural networks.

Multi-layer perceptron

A multilayer perceptron (MLP) is a feedforward neural network composed by at least three layers of neurons where the output of a neuron of a layer becomes the input of a neuron of the next layer. In a MLP, each neuron of a layer is linked to all the neurons of the next layer but has no link with the neurons of the same layer. The mathematical formulation of a MLP with L hidden layers reads as follows:

We set a (0) (x) = x.

For k = 1, ..., L (hidden layers)

z (k) (x) = b (k) + W (k) a (k-1) (x), a (k) (x) = φ(z (k) (x)). (3.4) For k = L + 1 (output layer) z (k) (x) = b (k) + W (k) a (k-1) (x), a (k) (x) = ψ(z (k) (x)). (3.5)
where φ is the activation function of the hidden layers and ψ is the output layer activation function. Superscript (i) stands for the i-th layer. At each step, W (k) is a matrix with number of rows the number of neurons in the layer k and number of columns the number of neurons in the layer k -1.

The number of parameters of such a network is equal to the sum of the multiplications of the number of neurons between each consecutive layer, which in practice can be very high. Moreover, since the input of an MLP must be a vector, depending on the dimension of the problem, this approach might be inefficient when dealing with images or higher dimension inputs.

Convolutional neural networks

In this section we will present Convolutional neural networks (CNN). For the sake of simplicity we will explain them in 2 dimensions. The main advantage of CNNs is that they work directly with tensors of any dimension. As depicted in figure 3.2, a CNN is made of several types of layers: convolutional layers, pooling layers and fully connected layers. In the following, each type of layer will be detailed. CNN to classify handwritten digits. Source: https:// towardsdatascience.com Convolutional layer A convolutional layer is composed of 2 main parts: the object and the filters. Filters or kernels are matrices of learnable weights that extract the features from an input signal. The discrete convolution of two functions f and g reads as:

(f * g)(x) = t f (t)g(x + t).

(3.6)

For a 2D image I, the 2D-convolution with convolution kernel K is given by:

(K * I)(i, j) = m,n K(m, n)I(i + n, j + m). (3.7)
The convolution kernel is dragged on the image. At each position, we compute the convolution between the kernel and the patch of image over which the kernel is hovering. Then, the kernel slides by a number s of pixels (s being the stride) till it parses the complete image. To control the size of the output, zero padding can be added around the image (margin of zeros).

After the convolution operations, a ReLu activation function is generally added.

If K is a convolution kernel of size k × k and x is a patch of the image of size k × k pixels, the activation is obtained by sliding the k × k window and computing z(x) = φ(K * x + b) where b is a bias.

Several convolution layers can follow one another where the output of a convolution becomes the input of the next one.

Pooling layer Pooling layers reduce the spatial dimension of the convolved features by taking the mean or the maximum on patches of the input (mean-pooling or max-pooling). These layers also act on small patches so there is also a pooling stride. This step allows to extract the dominant features which are rotational and positional invariant.

Fully connected layer A CNN generally ends with one or various fully connected layers (e.g. perceptron layers) in order to learn nonlinear combinations of the highlevel features extracted by the convolutional layers. The output of the convolutions must be flattened into a column vector before being fed to the dense layers.

In the same manner as for MLP, depending on the problem, a final activation function can be added (Softmax for classification for example).

Important remark

CNNs have been widely used to extract features from images consisting of 2D or 3D matrices of values that are spatially connected. Note that the use of convolutional filters expects the input to be encoded in a grid structure as such filters have a spatial representation of the inputs. This is a key point to keep in mind for next chapters, where we will compute the deformations of an object using CNNs. Indeed, we cannot directly feed a CNN with the mesh of an object. Instead, we will first have to express the physical quantities of interest in a regular grid structure. We will explain in detail this procedure in chapters 4 and 7.

Training a neural network

Once the architecture of the network is set, the optimal parameters (weights and biases) must be found. To do that we perform a gradient descent minimization of a loss function that is chosen by the user. The general procedure for training a network consists of 7 main steps:

1. Initialize the weights with close to zero values;

2. Feed the network with an input;

3. Forward propagation: neurons are activated depending on their weights. Spread activations until the prediction is computed;

4. Compute the error between the prediction and the expected value using the loss function;

5. Backpropagation: spread the error back into the network. Update the weights such that the error is minimized and adjust the learning rate;

Artificial Neural Networks: Fundamentals 29 6. Repeat steps 1 to 5 and adjust weights after a batch of samples;

7. When the entire data set has been seen by the network, it is called an epoch. Repeat more epochs.

In next section we will present the different loss functions available and some variations of the gradient descent algorithm. Then we will explain the backpropagation algorithm.

Loss functions

To learn, the network must know the committed error over each batch. To do that, the gap between the prediction f (x, θ) and the expected value y is measured using a loss function L. There are several options for the loss function depending on the considered problem. A common way for finding the parameters of a network is to maximize the likelihood (or the log likelihood), which is equivalent to minimizing the loss function (which is the opposite of the log likelihood). For a distribution D of pairs (X, Y ), the expected loss function reads as

L(θ) = -E (X,Y )∼D (log(p θ (Y /X)) . (3.8) If the model is Gaussian, that is to say if p θ (Y /X = x) ∼ N (f (x, θ), I)
, maximizing the likelihood is equivalent to minimizing the quadratic loss

L(θ) = E (X,Y )∼D Y -f (X, θ) 2 .
(3.9)

In the case of binary classification, maximizing the log likelihood is equivalent to minimizing the cross-entropy. Setting f (X, θ)

) = p θ (Y = 1/X), L(θ) = -E (X,Y )∼D [Y log(f (X, θ)) + (1 -Y )log(1 -f (X, θ))] . (3.10)
This loss function is well adapted to the Sigmoid activation as the logarithm prevents small values of the gradient. Finally, for a multi-class classification problem, a generalization of the previous loss function to k classes is considered

L(θ) = -E (X,Y )∼D   k j=1 1 Y =j log(Y = j/X)   . (3.11)
Alternatively, one could encode the knowledge of the data domain in the cost function. For instance, in the case of modeling biomechanics with a NN, the cost function can correspond to the strong formulation of the elasticity equations [START_REF] Raissi | Physicsinformed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF]. We will not enter into detail for this type of loss functions as they are not consider throughout this manuscript.

Artificial Neural Networks: Fundamentals 30

Minimization with Stochastic gradient descent

To minimize a cost function, the gradient descent algorithm is one of the most popular methods. In the case of neural networks, the minimization is performed over the network parameters (weights and biases). For the sake of visual clarity, in the following the biases {b j } j=1,...r are contained in the weights {w j } j=1,...n .

The gradient descent canonical formula reads as: θ := θ -α.∇L(θ), (3.12)

where θ = [w 1 , w 2 , ..., w n ] T is the parameters vector, α is the learning rate and ∇L(θ) is the gradient of the cost function such that ∇L(θ) = ∂L ∂w 1 , ∂L ∂w 2 , ..., ∂L ∂wn . The learning rate α is an hyperparameter of the minimization. It determines the step size at each iteration while moving towards a minimum. It impacts the learning speed and the convergence. A large value of α might lead to fluctuation around an optimum thus making the optimization diverge (can result in underfitting). On the contrary, a too small value of α can slow down the convergence speed and converge to a local minimum. In this case, the network can adjust too well to the training data set and it is no longer able to generalize (overfitting). The optimal value of α cannot be analytically calculated for a given model on a given dataset. It is usually advised to adapt it during training but there is however no general rule. It is usually acknowledged that it should be initialized at a large enough value (0.1 for instance) and progressively be reduced during successive iterations of the training process.

Gradient descent needs the cost function to be convex. When this is not the case, taking the steepest slope is no longer enough (due to local minimums and saddle points). The risk is therefore to drastically slow down the backpropagation and to fall on a local minimum. The Stochastic gradient descent (SGD) algorithm addresses this issues. Its canonical formula is given by: θ := θ -α.∇L(θ; x (i) ; y (i) ).

(3.13)

Let {(x n , y n )} N n=1 be the training data set made of N samples. Instead of computing the gradient over all the samples of the data set, SGD computes the gradient over one randomly selected sample or subset of samples called a batch. At each step, m training samples are randomly selected without replacement and a combination of the m corresponding gradients is computed (usually the average). m is called the batch size. This process is repeated until the network has seen the entire training data set: this is called an epoch. The maximal number of epochs is usually fixed in advanced and the training stops when this value is reached or when an accuracy criterion is satisfied (early stoping). SGD is faster than standard (deterministic) gradient descent and provides more fluctuations of weights, which augments the probability of not being stuck in a non optimal local minimum.

Due to the high sensitivity of the SGD to the learning rate, variations of the algorithm have been proposed with adaptive learning rates. One of the most famous is the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with adaptive moment estimation. The method computes individual adaptive learning rates for different parameters from estimates of first and second moments of the gradients.

Backpropagation

The backpropagation of the gradient was proposed by [ Le Cun, 1988] as an efficient way to compute the gradient of a neural network through a method called chain rule. After each forward pass through a network, backpropagation performs a backward pass to adjust the network's parameters computing the gradient one layer at a time, iterating backward from the last layer to avoid redundant calculations of intermediate terms in the chain rule.

In this section we will explain how to compute the gradient of Equation (3.9) using the backpropagation algorithm. The empirical quadratic loss of Equation (3.9) is proportional to (3.15) where K is the size of the network's output. In a regression model, the output activation function ψ is generally the identity function, to be more general, we assume that

N i=1 R i (θ), (3.14) Artificial Neural Networks: Fundamentals 32 with R i (θ) = K k=1 (Y i,k -f k (X i , θ)) 2 ,
ψ(z 1 , ..., z K ) = (g 1 (z 1 ), ..., g K (z K )), (3.16)
where g 1 , ..., g K are functions from R to R. Let us compute the partial derivatives of R i with respect to the weights of the output layer. Recalling that

z (L+1) (x) = b (L+1) + W (L+1) a (L) (x), (3.17) we get ∂R i ∂W (L+1) k,m = -2(Y i,k -f k (X i , θ))g k (z (L+1) k (X i ))a (L) m (X i ). (3.18)
Differentiating now with respect to the weights of the previous layer

∂R i ∂W (L) m,l = -2 K k=1 (Y i,k -f k (X i , θ))g k (z (L+1) k (X i )) ∂z (L+1) k (X i ) ∂W (L) m,l . (3.19) with z (L+1) k (x) = j W (L+1) k,j a (L) j (x), (3.20) a (L) j (x) = φ b (L) j + < W (L) j , a (L-1) (x) > . (3.21)
This leads to ∂z

(L+1) k (x) ∂W (L) m,l = W (L+1) k,m φ b (L) m + < W (L) m , a (L-1) (x) > a (L-1) l (x). (3.22) Let us introduce the notations δ k,i = -2(Y i,k -f k (X i , θ))g k (z (L+1) k (X i )) s m,i = φ z (L) m (X i ) K k=1 W (L+1) k,m δ k,i .
(3.23)

Then we have

∂R i ∂W (L+1) k,m = δ k,i a (L) m (X i ), ∂R i ∂W (L) m,l = s m,i a (L-1) l (X i ). (3.24)
known as the backpropagation equations.

Artificial Neural Networks: Fundamentals 33

The values of the gradient are used to update the parameters in the gradient descent algorithm. At step r + 1, we have:

W (L+1,r+1) k,m = W (L+1,r) k,m -α r m i=1 ∂R i ∂W (L+1,r) k,m W (L,r+1) m,l = W (L,r) m,l -α r m i=1 ∂R i ∂W (L) m,l , (3.25)
where α r is the learning rate that satisfies α r → 0, r α r = ∞, r α 2 r < ∞, for example α r = 1/r. We use the backpropagation equations to compute the gradient by a two pass algorithm. In the forward pass, we fix the value of the current weights θ (r) = (W (1,r) , b (1,r) , ..., W (L+1,r) , b (L+1,r) ), and we compute the predicted values f (X i , θ (r) ) and all the intermediate values (z (k)

(X i ), a((k)(X i ) = φ(z (k) (X i ))) 1≤k≤L+1
that are stored. Using these values, we compute during the backward pass the quantities δ k,i , and s m,i and the partial derivatives given in Equation (3.23). We have computed the partial derivatives of R i only with respect to the weights of the output layer and the previous ones, but we can go on to compute the partial derivatives of R i with respect to the weights of the previous hidden layers.

In the backpropagation algorithm, each hidden layer gives and receives information from the neurons it is connected with. Hence, the algorithm is adapted for parallel computations. The computations of the partial derivatives involve the function φ , where φ is the activation functions. φ can generally be expressed in a simple way for classical activation functions. The backpropagation algorithm is also used for classification with the cross-entropy but we will not enter into details for this aspect.

Conclusion

The objective of this thesis is to train neural networks for estimating physical quantities such as displacement fields or forces applied to an object. For this reason, in this chapter we have presented the different network architectures that we will use, and defined their constitutive elements for a better understanding of the approach that will be presented in next chapters. As seen in Section 3.3.3, the backpropagation equations involve a lot of tedious calculations that are fortunately automatically handled by frameworks such as Tensorflow 1 and Pytorch 2 that are both used during this thesis. In order to train the networks, large amounts of high fidelity data are needed which is the main difficulty of the approach, in particular when considering medical applications where acquiring thousands of real samples can be delicate. In the following chapters, we will demonstrate how to use FE simulations as accurate data generators to overcome the lack of real patient data. As seen in Section 3.2.2, the generated data must be expressed in a "network interpretative manner", namely

Chapter 4

The U-Mesh Framework

Introduction and related work

As said in the general introduction, there are many applications in engineering where nonlinear deformations need to be simulated in real-time. Some important examples can be found in the field of medicine but also in other areas of mechanical engineering, such as training of complex industrial processes [START_REF] Amundarain | Virtual reality for aircraft engines maintainability[END_REF] or virtual prototyping [START_REF] Barbic | Six-dof haptic rendering of contact between geometrically complex reduced deformable models[END_REF] just to name a few. While there are different numerical strategies for solving equations associated with elastic materials, we only consider the finite-element method (FEM), for its accuracy and ability to simulate a large range of materials on potentially complex domains. However, obtaining real-time simulations with this method when considering nonlinear materials, becomes a real challenge, in particular if this has to be done on consumers level hardware rather than a high-end parallel computer.

In order to speed up FEM simulations, several techniques have been proposed. We review here only some of the main ideas that were proposed. First, since solving the system of equations resulting from the FEM discretization is usually the bottleneck of the computation, many works have focused on linear solvers. Domain decomposition methods are based on the "divide & conquer" paradigm. Such methods consist of splitting the global problem domain into smaller independent sub-domains, making the approach suitable for parallel computing. This allows to build efficient preconditioners even though additional computation is required to synchronize the solution between neighboring sub-domains. Under the right conditions, in particular if the number of processors of the computer matches the number of sub-domains, a superlinear speedup can be obtained [START_REF] Haferssas | An Additive Schwarz Method Type Theory for Lions's Algorithm and a Symmetrized Optimized Restricted Additive Schwarz Method[END_REF]. This is, however, impossible to achieve when considering problems (even if relatively small) which need to be solved in real-time on consumers level hardware. This is mainly due to the limited number of cores (only 10 cores on the latest Intel i9 processor) and communication costs which are significant compared to the expected computation times (about 50 ms per time step for an interactive simulation).

Another option for speeding up simulation times is to lower the computational complexity of the problem through a reduction of the model's degrees of freedom. Depending on the problem, and acceptable loss of accuracy, it is possible to obtain speedups of several orders of magnitude. Proper Orthogonal Decomposition (POD) is one of the main model order reduction methods. POD techniques compute off-line the solution to several complete models and extract the modes that describe best the solution to the complete problem. Based on a priori knowledge of the solution, it encodes the high dimensional problem in a smaller subspace defined by a truncated basis of singular vectors. The dimension of such basis determines the ratio between accuracy of the method and computation times (the smaller the basis, the larger the error). In the context of real-time simulation of nonlinear solids, several examples of POD have been proposed. [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF] proposed a POD method to simulate the palpation of the cornea. Haptic feedback rates were achieved, but with a relative error of about 20%. If more accurate solutions are needed, the number of modes used in the POD can be increased at the cost of higher computation times. Thus, to keep the method numerically efficient in the case of nonlinear materials, hyperreduction can be used in order to further reduce the computation times while reducing the error [Ryckelynck, 2005]. [START_REF] Goury | Fast, Generic, and Reliable Control and Simulation of Soft Robots Using Model Order Reduction[END_REF] applied the hyperreduced POD to control and simulate soft robots with very good accuracy in 25 ms per time step. However the POD may be in some cases insufficient to capture correctly the high degrees of nonlinearity that can be found for example in biological soft tissues as it relies on a linear combination of few basis vectors [START_REF] Bhattacharjee | A nonlinear manifold-based reduced order model for multiscale analysis of heterogeneous hyperelastic materials[END_REF]. To precisely account for nonlinearities it might be necessary to recompute the entire stiffness matrix which is burdensome and not always possible [Niroomandi et al., 2017]. Another model order reduction algorithm is the proper generalized decomposition (PGD), which, contrarily to POD, builds a reduced-order approximation without relying on the knowledge of the solution of the complete problem. PGD assumes that the solution of a multiparametric problem can be expressed as a sum of separable functions that are constructed by successive enrichment by invoking the weak form of the considered problem. An approach based on PGD is proposed in [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF] for the simulation of hyperelastic soft tissue deformation at high frequency. However, when the solution is non-separable, PGD offers no particular advantage over classical FEM techniques.

A last class of worth mentioning solutions consists of using the Graphics Processing Unit (GPU) as a particular type of parallel machine. Although each core of the GPU is very limited in its computational performance, the very high number of cores available (several thousand) makes it possible to obtain significant speedups on computationally heavy problems. For instance, NiftySim [START_REF] Johnsen | NiftySim: A GPU-based nonlinear finite element package for simulation of soft tissue biomechanics[END_REF]] is a GPUbased nonlinear finite element toolkit for the simulation of soft tissue biomechanics where speedups of 300x are obtained. SOFA [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF] is an Open-source

The U-Mesh Framework 38 Framework focused on real-time simulation of complex interactions with deformable structures, which provides GPU-compatible FEM codes [START_REF] Comas | Efficient nonlinear FEM for soft tissue modelling and its GPU implementation within the open source framework SOFA[END_REF].

Recently, machine learning started to revolutionize several fields (vision, language processing, image recognition, genomics) due to the continuously increasing amount of data available and the development of new algorithms and powerful GPUs. Deep learning, a class of machine learning methods based on learning data representations, as opposed to task-specific algorithms, has demonstrated strong abilities at extracting high-level representations of complex processes. With sufficient ground truth data, machine learning algorithms can map the input of a function to its output without any mathematical formulation of the problem, thus actuating like a black box. The term black box relates to the non-clarity in the way the model uses training data to reach particular conclusions. These ambiguities make it difficult to know the reasons why the system behaves the way it does. In this chapter we will put a lot of effort in clarifying such ambiguities by exploring connections of the chosen architecture with model order reduction techniques. To fulfill the need for ground truth data, since FEM can provide as much noise-free data as required, it seems interesting to train learning algorithms with such virtually generated data [START_REF] Lorente | A framework for modelling the biomechanical behaviour of the human liver during breathing in real time using machine learning[END_REF], Luo et al., 2018, Roewer-Despres et al., 2018, Tonutti et al., 2017, Fetene et al., 2018, Runge et al., 2017].

For example [START_REF] Lorente | A framework for modelling the biomechanical behaviour of the human liver during breathing in real time using machine learning[END_REF] proposed a machine learning approach for modeling the mechanical behavior of the liver during breathing in real-time. They trained several regression models using external displacement and elasticity parameters as input, and the FEM-based nodal displacements as output. Although they reached good accuracy their method is restricted to small displacements. [START_REF] Roewer-Despres | Towards finite element simulation using deep learning[END_REF] proposed a preliminary work using a deep-autoencoder to approximate the deformations of a nonlinear muscle actuated object. They showed that their method produces lower reconstruction errors than the equivalently sized PCA model. However the computational gain of their method is not clear and it was limited to a very simple model and coarse mesh. [START_REF] Tonutti | A machine learning approach for real-time modelling of tissue deformation in imageguided neurosurgery[END_REF]] treated a simple problem using two different networks, one to predict the magnitude of the displacement and the other one to predict its direction. However, it seems that these two networks require a specific training for each node of interest, which could be prohibitive for large meshes. Moreover their model is limited to small deformations on relatively simple shapes, with restricted input forces. For instance, the considered displacements never exceed 5 mm (for an organ of size 20 cm) and only 11 nodes of the mesh are excited. On the contrary our approach can handle complex and complete volume deformations of arbitrary shapes with one single network for any force application point. All the cited references propose to train neural networks with FEM generated data for various purposes. However, none of them is justifying the choice of the network architecture used to do this. In this chapter, we explore connections of

The U-Mesh Framework 39 the chosen architecture with model order reduction techniques in order to support interpretability and explainability to go beyond the black-box usage of neural networks.

The objective of our work, whose preliminary results are presented in this chapter, is to go beyond the state of the art on machine learning applied to computational mechanics. In particular we show that we can predict, in real-time, the shape of a nonlinear elastic structure with a very good accuracy using a deep network inspired by model order reduction methods. Our solution relies on a U-Net architecture trained on FEM-generated data sets; both aspects are presented in Section 4.2 while results and their comparison against a model order reduction algorithm are presented in Sections 4.3.

Displacement field estimation with deep neural networks

As mentioned previously, an important area of applications for real-time simulation of nonlinear material is in the field of computer-aided surgery where accuracy is also very important, but often left aside for the sake of rapidity. Achieving a better trade-off between accuracy and computation time is therefore mandatory to tackle more ambitious problems. From this chapter on, we propose a method that does not need such a compromise. It allows for extremely fast and accurate simulations by using an artificial neural network that partially encodes the stress-strain relation in a low-dimensional space. Such a network can learn the desired biomechanical model, and predict deformations at haptic feedback rates with very good accuracy. This section is divided in two main segments. First the selected network architecture is detailed, followed by our strategy to encode the stress-strain relationship and boundary conditions into the network and the data set generation used to train the network.

Network architecture

Formally, our network h is a parameterized function that accepts a 3 × n x × n y × n z tensor of input constraints C and produces a tensor of volume displacements U v of the same size as output. The domain Ω is sampled by a 3-dimensional grid of resolution n x × n y × n z . Practically speaking, the nodes of this grid match the nodes of the FEM mesh, although this is not required and an interpolation could be used instead. The tensor C represents the constraints applied over the surface boundary Γ of the domain. In particular, C can be seen as a traction t applied to Γ N or as an imposed surface displacement U s on Γ D . The tensor U v contains the

The U-Mesh Framework 40 volumetric displacement at the nodes of the mesh in response to the constraint C. In particular, each vector C[:, i, j, k] represents the constraint vector (c x , c y , c z ) applied over the node (i, j, k) of the grid. Similarly, the vector U v [:, i, j, k] represents the displacement (u x , u y , u z ) of each node (i, j, k) of the grid.

Our problem consists of finding the function h that produces the best estimation of the displacement field given prescribed constraints C (in this chapter we only deal with contact forces e.g. traction). This is performed by minimizing the expected error over the distribution D of pairs (C, U v ):

min θ E (C,Uv)∼D h(C) -U v 2 2 , (4.1)
where θ is the set of parameters of the network h. In practice, the expectation of Equation (4.1) is approximated by Monte-Carlo sampling with a training set

{(C n , U vn )} N n=1 of N samples: min θ 1 N N n=1 h(C n ) -U vn 2 2 . (4.2)
We build our training set by randomly applying forces on the mesh and running FEM simulations to produce corresponding displacements.

Let us characterize the architecture chosen for our network h. We propose to use the U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], a modified fully convolutional network initially built for precise medical image segmentation. Note that the use of convolutional filters expects the input to be encoded in a grid, due to their spatial representation of the inputs. As depicted in Figure 4.1, the network is similar to an auto-encoder, with an encoding path to transform the input space into a low-dimensional representation, and a decoding path to expand it back to the original size. Additional skip connections transfer detailed information along matching levels from the encoding path to the decoding path.

The encoding path consists of k sequences of two padded 3×3×3 convolutions (k = 4 in [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]) and a 2 × 2 × 2 max pooling operation (see Figure 4.1). Intuitively, each 3D convolution filter learns to isolate the different characteristics of the displacement field (orientation, direction, amplitude). At each step, each feature map doubles the number of channels and halves the spatial dimensions. We assume that the number of channels is directly related to the amount of detectable variations in the displacement field. In the bottom part there are two extra 3 × 3 × 3 convolutional layers leading to a (c × 2 k )-dimensional array. This feature space is intuitively similar to the Galerkin projection of the equations of motion onto the reduced space in POD, where the order of the singular vector truncation is equivalent to the number of neurons in the latent space. A difference however remains with

The U-Mesh Framework 41 The U-Mesh Framework 42 the presence of convolutional operations at each layer. In a symmetric manner, the decoding path consists of k sequences of an upsampling 2 × 2 × 2 transposed convolutions followed by two padded 3 × 3 × 3 convolutions. The features from the encoding path at the same stage are cropped and concatenated to the upsampled feature maps. At each step of the decoding path, each feature map halves the number of channels and doubles the spatial dimensions. There is a final 1×1×1 convolutional layer to transform the last feature map to the desired number of channels of the output (3 channels in our case).

x x 0 y 0 x z 0 step 0 3 c x x 0 y 0 x z 0 x x 1 y 1 x z 1 x x 1 y 1 x z 1 x c 2 1 step 1 step k x x k y k x z k x x k y k x z k x c 2 k-1 x c 2 k x x
The number of steps k and the number of channels c control the accuracy of the prediction just like the number of singular vectors in POD. Higher values of k and c lead to a more complex network suitable for difficult meshes at the expenses of longer computational times for both training and prediction, and higher requirements of training data. We tested several values for k and c in order to select the most appropriate values for our experiments depending on the desired accuracy and the eventual time restrictions.

Offline training data generation

Without lack of generality, we consider the boundary value problem of computing the deformation of a hyperelastic material under both Dirichlet and Neumann boundary conditions described in Chapter 2. Hexahedral (H8) elements are used to discretize the computational domain. This choice is not only motivated by the good convergence and stability of such elements: hexahedral elements are also required for our convolutional neural network (see Section 4.2.1).

In order to train such a network, we build a data set of pairs (C, U v ) obtained with the FEM. Once we are given a 3D mesh with its corresponding constitutive law, material properties and boundary conditions, we perform multiple simulations by applying random forces to nodes of the object. After each simulation, the pair of applied forces and obtained deformation is stored as an element of the data set. To speed up the generation of the training and testing data sets, the linearized system of Equation (2.40) is solved using an iterative preconditioned conjugate gradient method [Shewchuk, 1994].

The variability of the data relies on the force magnitude, its direction and its application point. In this chapter, a force is applied on a local surface area Γ N whose location varies such that the boundary of the computation domain is completely covered. The force direction is uniformly sampled on the unit sphere and the force magnitude is a uniform random value between 0 and 1. At each sample of the data set, one force is applied on a small region Γ N . There are Λ samples for each Γ N ,

The U-Mesh Framework 43 meaning that Λ different forces are applied per region. We can consider one force per sample or several forces applied simultaneously at different locations.

The network is trained minimizing Equation (4.2) with training data generated as explained above. The minimization is performed using the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] (see Chapter 3 for details).

Results on synthetic objects

In this section, we perform a model selection over the space of hyperparameters k and c in order to find the combination that leads to the best results in a cantilever beam. Then we apply our method, with the selected hyperparameters, to two benchmark examples: a cantilever beam and an L-shaped object under point loads. All our experiments are performed in a GeForce 1080 Ti using a batch size of 4 and 100, 000 iterations for training. We use a PyTorch implementation of the U-Net. We recall that the batch size is the number of samples that are given to the network at each iteration of the minimization process.

Validation metrics

To assess the efficiency of our method, we perform a statistical analysis of the mean norm error e over a testing data set {(C m , U vm )} M m=1 , which is built similarly as the training data set. Note that the training data set and the testing data set are disjoint. Let U vm be the ground truth displacement tensor for sample m generated using the FEM described in Section 2.2 of Chapter 2 and h(C m ) the U-Mesh prediction. The mean norm error between U vm and h(C m ) for sample m reads as:

e(U vm , h(C m )) = 1 n n i=1 |U v i m -h(C m ) i |. (4.3)
where n is the number of degrees of freedom of the mesh. We compute the average e and standard deviation σ(e) of such norm over the testing data set:

e = 1 M M m=1 e(h(C m ), U vm ), (4.4)
and

σ(e) = 1 M -1 M m=1 [e(h(C m ), U vm ) -e ] 2 .
(4.5)

The

U-Mesh applied to a cantilever beam

We consider a deformable beam of size 4 × 1 × 1 m 3 subjected to fixed boundary on one end. The beam follows the Saint-Venant-Kirchhoff behavior described in Section 2.1.3 of Chapter 2 with a Young's modulus E of 500 Pa and a Poisson's ratio of 0.4, and is discretized with 135 H8 elements. To generate the data set, 100 forces are applied on each of the 64 nodes of the upper face (that is Λ = 100), up to a total of 6, 400 samples. 80% of the samples are used for training (N = 5, 120) and the remaining 20% are kept for testing (M = 1, 280). Using this data set, we select the best architecture by training several U-Nets with different combinations of hyperparemeters k and c. In Table 4.1 are reported the training and prediction times as well as e and σ(e). As seen in this table, the higher the feature space size (FSS), the lower the errors. Prediction time is proportional to the depth of the network.

Choosing the best model is a tradeoff between network performance and speed, and we will show that the selected hyperparameters lead to good results also on different problems. Choosing k = 3 and c = 128 seems to be a good compromise for our needs. The selected parameters appear in bold text in Table 4.1. For the selected set of parameters, e over the data set is equal to 0.0007 ± 0.0006 m for a maximal deformation of 0.724 m. In Figure 4.3, we show the sample with maximal error. We perform a sensitivity analysis of the method to the amplitude of deformation.

The results are plotted in Figure 4.4. We perform a least squares line fit to find the relation between the maximal deformation and the mean norm error e. We can observe a very small sensitivity of the error e with respect to the deformation amplitude, and a very small error in the estimation of the displacement field in general. This shows an important characteristic of our method, in addition to its very limited computational cost.

In Table 4.2 are shown the prediction errors for 3 simultaneous input forces and their corresponding training and prediction times. In this scenario, the Young's modulus is set to 400 P a and three forces are applied simultaneously on three different regions

The U-Mesh Framework 45 Figure 4.5: Largest deformation for 3 simultaneous force applications (side and front view for the same sample). U-Mesh output is in blue, the reference is in green and the rest shape is in grey.

of the upper face of the beam. In order to avoid mechanical coupling, such regions must be far enough from each other. There are 12, 360 possible combinations of nodes fulfilling the above stated condition. Only one force (Λ = 1) is applied per viable combination up to a total of 12, 360 samples (N = 9, 888 samples for training and M = 2, 472 samples for testing). We can see that the errors and the time needed for the prediction are comparable to that needed for only one force application.

It is important to note that there are 2x more samples in this data set (than in previous ones) due to the large number of possible combinations. This explains the particularly low error in this case. The sample with maximal deformation (1.0035 m) is shown in Figure 4.5. At the tip of the beam, the relative l2 norm is equal to 1.5%.

So far we have seen that U-Mesh is able to predict the displacement field due to one or several forces with high accuracy and in an extremely short amount of time.

Nevertheless, FEM codes are also able to compute the solutions on such meshes in limited computation times, even for hyperelastic models. Hence, in order to put ahead the real contribution of our work, we test our method on a computationally expensive problem.

We consider the same beam as previously but this time discretized in 3, 267 H8 elements (see selected parameters, that is, 128 channels for the first layer and 3 steps. The metrics computed over the testing data set are shown in Table 4.3. The most interesting result is the low prediction time (only 3 ms). A very optimized version of the Saint-Venant-Kirchhoff FEM using a Pardiso solver [START_REF] Kourounis | Toward the next generation of multiperiod optimal power flow solvers[END_REF] that is among the most efficient solvers available, takes more than 300 ms to solve such simulation.

The speedup obtained with U-Mesh is of 100x. All the samples of the testing data set have an error bellow 0.0265 m and an average error of 0.0019 m for a maximal deformation of 1.011 m.

In the following paragraphs we will show that U-Mesh generalizes well on other geometries. In particular, we will see the performance of U-Mesh applied to an L-shaped object.

U-Mesh applied to an L-shaped object

We apply U-Mesh on an L-shape of size 28.424 × 10 × 40 m 3 discretized in 335 H8 elements. Since the U-Net requires a regular grid as input, the L-shaped object is embedded in a regular grid (with zero-padding). The Young's modulus is equal to 500 P a and the Poisson's ratio is 0.4. To build the data set, external forces ranging from 0 to 40 N are applied on the bottom face of the L (with Λ = 100) up to a

The U-Mesh Framework 48 total of 6, 000 samples (N = 4, 800 and M = 1, 200). We train a U-Net with 128 channels in the first layer and 3 steps. In Figure 4.8 are shown some samples of deformed L-shapes. The U-Mesh output is in blue and the reference solution is in green. In order to perceive the amount of deformation, the rest position is also shown (thin grey lines). The average and maximal errors are given in Table 4.4, and the prediction times are in the same range as for the beam scenario. The average error is equal to 0.00648 m where the maximal deformation is 8.9016 m. The slope of the regression in Figure 4.9 shows that the increase of the error with the deformation amplitude is controlled. It is worth noting that the outliers of this graph (such as the one marked in red) still correspond to small errors (see Figure 4.10).

Comparison of U-Mesh and POD

In this section we compare the predictions made by U-Mesh to the simulations computed on a reduced model using POD. We used the POD code available at https://github.com/SofaDefrost/ModelOrderReduction that works as a plugin of the SOFA framework [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF]. The POD consists of three phases.

First, an offline phase where all the potential movements of the beam are sampled and stored in the so-called snapshot. This offline phase is the equivalent to the data generation phase in U-Mesh and is also computationally intensive since it performs many fine simulations. In a second phase, the snapshot space is condensed in a reduced

The U-Mesh Framework 49 basis using singular value decomposition and a truncation [START_REF] Goury | Fast, Generic, and Reliable Control and Simulation of Soft Robots Using Model Order Reduction[END_REF]]. This phase "corresponds" to the training of the U-Net and is generally faster. Finally, the resulting reduced model allows for faster simulations since there are fewer degrees of freedom. We applied the simple POD and the hyperreduced-POD (HPOD) to the fine beam depicted previously (see Figure 4.7) in order to compare the performance of POD and U-Mesh.

We first compare the computation times for a given accuracy. The truncature error of the POD was set such that the mean norm error obtained with POD is similar to the one obtained with U-Mesh. To reach the desired accuracy in the considered deformation range, 3 modes were preserved. Computation times are reported in Table 4.5. With the selected number of modes, POD is about 6 times faster than the full FEM model whereas U-Mesh is more than 200 times faster than the full FEM model.

Let us now compare the relative errors of the two methods for a given computation time. The fastest version of the reduced model is the one considering only one deformation mode and using hyperreduction. As presented in Table 4.6, HPOD can compute deformations in 5 ms but with an error that is 14 times larger than the one produced by U-Mesh.

Conclusion

In this chapter we have proposed a U-Net architecture that can learn the relation function between an input force and an output deformation for regular geometries and make predictions with high accuracy in a record amount of time. The take-homemessage of this chapter is that for a given network architecture, the prediction time is nearly constant (and very short), regardless of the size of the problem. Furthermore the accuracy of the prediction, which depends on the quality and the size of the data set, is controllable since we generate this data. We believe that such an approach Despite the promising results of our method, with the current U-Mesh, it is not possible to make very accurate predictions when applying a force somewhere out of the sampled input domain of the training data set. This limitation is also inherent to POD. In the same manner, we are restricted to the geometry used to train the network. This is due to the fact that neural networks are not good at extrapolating, hence the importance of a good sampling of the input domain when generating the training data sets.

Yet, there are several directions to investigate to make this approach more broadly usable. In particular in the context of biomechanics, the geometries are non regular structures, the inputs can be different from forces and the object characteristics (such as material parameters and boundary conditions) can vary from one patient to another. In next chapters we will see how to generalize the U-Mesh framework to patient-specific modeling of soft tissues.

Chapter 5

Patient-specific modeling

Introduction

In part I of this thesis work, we have explored how to simulate objects in realtime without taking care about the parameterization of the models. However, when considering computer-assisted interventions, the accuracy of the simulations and consequently their parameterization are essential. Depending on the surgical intervention, the demanded precision is variable but a target registration error of 5 mm is usually accepted [START_REF] Ruiter | Model-based registration of X-ray mammograms and MR images of the female breast[END_REF]. For instance, in image guided biopsies (where our mechanical model can be used to enrich the medical images), the tip of the puncturing needle must be in contact with the tumors whose size can be relatively small (less than 2 cm). To guarantee this level of precision the modelisation should be adapted to the patient anatomy and to its mechanical properties. As a matter of fact, the reliability of a biomechanical model highly depends on its parameterization. Therefore, once the constitutive law is chosen, the patient-specific material properties and boundary conditions must be carefully set in order to achieve the desired accuracy [START_REF] Bosman | The role of ligaments: Patient-specific or scenario-specific[END_REF].

Our finite-element simulations depend on parameters of the constitutive model such as Young's modulus, Poisson's ratio and the boundary conditions (BCs) (see Sections 2.3.1 and 2.3.2). Since soft tissues may be represented with an incompressible material, Poisson's ratio can be safely set to a value close to 0.5. However, the value of Young's modulus E is more difficult to estimate for a given organ as it varies with the age of the patient or even with treated pathology. For instance, considering hepatic surgeries, a cirrhotic liver is significantly stiffer than average livers. Therefore, values from the literature do not directly match each patient. Besides that, the location and the elastic properties of the attachments of the organ play a major role in the accurate approximation of the displacement field. Such BCs are not visible in the preoperative images and it is difficult to estimate them intraoperatively as they are often out of the field of view of the surgery.

The elastic properties of materials can be identified by solving inverse prob-lems [START_REF] Zhao | Pointwise Identification of Elastic Properties in Nonlinear HyperelasticMembranes-Part II: Experimental Validation[END_REF], Gee et al., 2010, Lu et al., 2009, Sinkus et al., 2010] or using elastography techniques [START_REF] Sarvazyan | Shear wave elasticityimaging: A new ultrasonic technology of medical diagnostics[END_REF], Muthupillai et al., 1995, Xu et al., 2007] initially developed for diagnosis purposes. Some works have focused on the estimation of BCs intraoperatively such as [START_REF] Peterlík | Model-based identification of anatomical boundary conditions in living tissues[END_REF], [START_REF] Plantefeve | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF] and [START_REF] Johnsen | NiftySim: A GPU-based nonlinear finite element package for simulation of soft tissue biomechanics[END_REF] but these methods are difficult to use in practice as either additional intraoperative scanning is required or they are sensitive to anatomical variations. Moreover, when acquiring information intraoperatively, observational errors may occur, thus introducing uncertainty to the system. Alternative solutions accounting for such uncertainty rise from the use of Bayesian methods. For instance, authors in [START_REF] Nikolaev | Stochastic Correction of Boundary Conditions during Liver Surgery[END_REF], Peterlík et al., 2017], employed the reduced-order unscented Kalman filter (ROUKF) to model the BCs of a liver as stochastic parameters, leading to more accurate simulations of the deformations of the organ.

In this chapter we will briefly see how to obtain a patient-specific geometry of the patient's anatomy and then, similarly to works in [START_REF] Nikolaev | Stochastic Correction of Boundary Conditions during Liver Surgery[END_REF], Peterlík et al., 2017], we will describe a method to determine the patient-specific parameters of an organ using the ROUKF algorithm.

Patient-specific geometry

The first step towards patient-specific modeling is the patient-specific geometry of the organ. Generally, the 3D anatomical model of the organ is constructed from preoperative volumetric medical images such as computerized tomography (CT) scan or magnetic resonance imaging (MRI). Both acquisitions consist of a sequence of 2D slices that are merged into a 3D image using volume rendering. The color intensity of each voxel translates variation in material properties thus allowing to distinguish the different structures. Then, a process called segmentation is performed with a medical imaging software such as 3D Slicer1 or with deep learning algorithms [START_REF] Zhou | A review: Deep learning for medical image segmentation using multi-modality fusion[END_REF] in order to extract a volumetric surface of the organ and its internal structures (see Figure 5.1 for an example).

Model parameterization

We present here the procedure to estimate the value of Young's modulus and the BCs of an organ using observations of the target model using the ROUKF. To this end, each sought parameter p is described as a stochastic parameter associated to a Gaussian probability density function (PDF). Initially p ∼ N (µ 0 , σ 0 ) with µ 0 the mean value of p reported in the literature and σ 0 its standard deviation. The aim of the assimilation process is to iteratively reduce the standard deviation σ of p in order to find the most likely value for µ. To this end, the PDF of the parameter p is transformed based on observations. The transformation of the PDF is modelled using a ROUKF which can handle non-linear processes, and is computationally efficient [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF].

ROUKF: Overview of the algorithm

Once the FE model of the organ is built (from a preoperative CT scan for instance), the constraints imposed on the surface Γ D need to be identified in order to generate a deformation. In other words, the organ attachments representing Dirichlet boundary conditions and the traction or displacement imposed on the free part of the boundary need to be identified. In this work, the former can be either fixed (Section 9.2.1) or set as a stochastic set of parameters (Section 9.2.2) and the latter can be determined intraoperatively. During the intervention, points in the surface of the organ can be tracked in each video frame. Such points are called features and are separated into control features and observation features. The control features govern the deformation of the liver model (imposed displacement on Γ D or traction if a force sensor is available) and the observation features correspond to ground truth data (used in the filter correction phase to compute the Kalman gain). The control features can be selected close to the surgical tool and be used to prescribe displacements in the mechanical model.

An efficient implementation of a Bayesian inference method able to process nonlinear systems like our models is the unscented Kalman filter (UKF) [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF].

Compared to an extended Kalman Filter, it does not require to compute the Jacobian of the system, which would be prohibitive given the size of our problem. The

Patient-specific modeling 59 unknown data to be estimated (the stochastic state of the system) is described as a Gaussian distribution, which transformation through the nonlinear system is performed using an unscented transformation (see [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF] for details). The main idea is to parameterize the Gaussian distribution using a set of sigma points, which hold the mean and covariance information, but are easier to transfer through a nonlinear function. The general algorithm is described in Algorithm 2. It consists of a loop that contains two main steps. During the prediction step we form the new hypothesis about the estimated state, while during the correction step we correct it by comparing the predicted measurements with (noisy and partial) observations. Compute prediction phase: 7:

x σ * i = x i + √ P i I -generate sigma points 8:
for each sigma point k do 9:

x σ k i+1 = T (x σ k i ) get result from deformation step 10:
end for 11:

x i+1 = E x σ * i+1 -compute predicted state as mean of sigma points 12:

P i+1 = x σ * i+1 -x i+1 x σ * i+1 -x i+1 T + Q -compute predicted covariance 13:
Compute correction phase:

14: get q (o)
i+1 -observation features 15:

for each sigma point k do 16:

q (o) i+1 σ k = H( x σ k i+1 ) -get predicted observation 17:
end for 18:

P xq (o) = x σ * i+1 -x i+1 q (o) i+1 σ k -E q (o) i+1 σ * T -compute cross covariance 19: P q (o) = q (o) i+1 σ k -E q (o) i+1 σ * q (o) i+1 σ k -E q (o) i+1 σ * T + W -comp. obs. cov. 20: K i+1 = P xq (o) P -1 q (o) -

compute Kalman gain

21:

x i+1 = x i+1 + K i+1 q (o) i+1 -E q (o) i+1 σ *
-compute corrected state 22:

P i+1 = P i+1 -P xq (o) P -1 q (o) P T xq (o) -compute corrected covariance 23: end for
The prediction step can be very costly when using a model with many degrees of freedom, as it is the case when using a FEM method. Using the simplex method to generate the sigma points would require m + 1 simulations if m is the number of elements in the stochastic state vector (line 9 of the algorithm). With a mesh of n nodes and k stiffness parameters, this would mean 3n + k + 1 simulations. A simple FEM mesh of only a few hundred nodes would be too time-consuming for a Patient-specific modeling 60 clinical application, as it would require more than 300 simulations for each step of the assimilation process. To solve this issue, we use a ROUKF instead of the UKF. This method significantly reduces the computation cost since only k + 1 simulations (in the best case) are required. This approach was proposed in [START_REF] Peterlík | Image-driven stochastic identification of boundary conditions for predictive simulation[END_REF].

Let us assume there are k unknown parameters in our model, so k different parameters to estimate that can be either the elasticity of the material or the elasticity of the organ attachments. Since we are using the simplex version of the ROUKF, there are k + 1 sigma points meaning that k + 1 evaluations of the model are performed in each prediction step of the assimilation process. At each step of the assimilation, the control features are extracted from the actual video frame and mapped onto the FE model through barycentric coordinates, in order to prescribe displacements. At the first step, µ and σ are initialized to µ 0 and σ 0 for each parameter. Then, k +1 vectors of parameters are sampled and k + 1 simulations are performed. Each simulation corresponds to one of the sampled k + 1 values of the parameter and they can be done in parallel as they are independent. After the simulations for all sigma points are performed, the a priori expected value and covariance matrix are updated. This is called the prediction phase. Later, in a correction phase, the extracted observation features are compared to the model predicted positions to compute the innovation that is used to compute the Kalman gain. The a posteriori expected values and the covariance matrix are computed based on the Kalman gain.

Conclusion

In this chapter we have seen the importance of the patient-specific parameterization of the finite-element models. The ROUKF allows for real-time estimation of the sought parameters based on intraoperative observations. In chapters 7, 8 and 9 we will see how to introduce this knowledge into the U-Mesh for allowing patient-specific predictions of the displacement field of an organ. In next chapter, we will see how this data assimilation process can help in robotic interventions.

Chapter 6

Force Estimation in Robotized Injections

Introduction and related work

The work presented in this chapter was produced at the beginning of this PhD. It consists of a thematically isolated example of combining finite-element simulations with neural networks in the context of robotized intravitreal injections.

Intravitreal injections are among the most common surgical interventions in ophthalmology with more than 4 million injections worldwide in 2014 alone [START_REF] Ullrich | Assistive Device for Efficient Intravitreal Injections[END_REF]. This procedure is principally used in the treatment of diabetic maculopathy and for injecting vascular endothelial growth factor inhibitors in the treatment of age-related macular degeneration. Besides, we observe an increasing demand for such therapy due to the growing prevalence of diabetic patients and aging demographics. Intravitreal injections are mostly performed by doctors and the cost of such therapy has to be reduced. The increasing workload and the reduced reimbursement makes it difficult for hospitals to handle the situation. The time spent by the clinicians performing the injections has to be minimized while preserving the precision and the safety of the patient.

At the same time, robotic assistance in ophthalmology provides the ability to improve manipulation skills, along with shorter and safer surgeries [START_REF] Meenink | Robot-assisted vitreoretinal surgery[END_REF]. To this end, [START_REF] Ullrich | Assistive Device for Efficient Intravitreal Injections[END_REF] proposed a robotized intravitreal device capable of assisting injections into the vitreous cavity. However, designing such robotic systems requires to solve multiple challenges in terms of safety, cost and time efficiency and usability. The position accuracy and the orientation of the needle are of particular interest since the injection must be performed in a small region (the pars plana) of the eye. If the region is missed, damage of the eye lens or the retina might occur. In addition to accurate positioning, the ability to estimate or measure the force exerted by the needle on the sclera during the procedure could offer an important additional safety for the patient [START_REF] Jagtap | Applied force during vitreoretinal microsurgery with handheld instruments[END_REF].

As a matter of fact the knowledge of the force plays a central role in the control loop of surgical robots for patient safety during robotic assisted interventions [START_REF] Weber | Instrument flight to the inner ear[END_REF]. The force can be either measured through force sensors or estimated through vision-based methods. The use of force-based control algorithms allows for an improved human-machine interaction, more realistic sensory feedback and telepresence. Beyond this, force sensing or force estimation can facilitate the deployment of essential safety features [Haidegger et al., 2017]. A considerable amount of work relying on force sensors has previously been done, focusing on the development of miniaturized devices to ease their integration with actual systems. Force sensors usually need to meet several additional requirements, such as being water resistant, sterilizable and insensitive to changes in temperature [START_REF] Haouchine | Vision-based Force Feedback Estimation for Robot-assisted Surgery using Instrument-constrained Biomechanical 3D Maps[END_REF]. The major limitation of conventional sensors is thus the associated cost since most surgical tools are disposable [Haidegger et al., 2017]. To overcome this point, alternative solutions have been proposed, such as qualitative estimation of forces based on images [START_REF] Mura | Vision-based haptic feedback for capsule endoscopy navigation: a proof of concept[END_REF], Haouchine et al., 2018]. [START_REF] Mura | Vision-based haptic feedback for capsule endoscopy navigation: a proof of concept[END_REF]] introduced a vision-based haptic feedback system to assist the movement of an endoscopic device using 3D maps generated with a Shape-from-Shading method where the 3D shape of the surface is recovered from a 2D image of that surface. [START_REF] Haouchine | Vision-based Force Feedback Estimation for Robot-assisted Surgery using Instrument-constrained Biomechanical 3D Maps[END_REF] incorporated the use of a biomechanical model of the organ in addition to the 3D maps to estimate force feedback in robot-assisted surgery. This approach is, however, limited to the ability to evaluate tissue properties of the organ quantitatively [START_REF] Aviles | A recurrent neural network approach for 3d vision-based force estimation[END_REF].

Deep learning has already been suggested to improve existent characteristics of robotic assisted surgeries such as instrument segmentation and detection [START_REF] Pakhomov | Deep residual learning for instrument segmentation in robotic surgery[END_REF], as well as force estimation.

For instance, in [START_REF] Aviles | A recurrent neural network approach for 3d vision-based force estimation[END_REF], interaction forces in minimally invasive surgeries are estimated with recurrent neural networks using camera acquisitions combined with kinematic variables and deformation mappings. In a follow-up paper, [START_REF] Aviles | Sensorless force estimation using a neuro-vision-based approach for robotic-assisted surgery[END_REF]] used a neuro-vision based approach for estimating applied forces in the same context. In [START_REF] Aggarwal | Covfefe: A Computer Vision Approach For Estimating Force Exertion[END_REF], authors used two neural networks to classify muscle force exertion levels to prevent musculoskeletal disorders based on features extracted from video data and blood volume changes. However, in all these approaches an intermediate step to determine the surface deformation is required.

In this chapter, we present a method for estimating contact forces directly from an OCT image of the scleral deformation without the need for a specific image feature extraction method beforehand, as in [START_REF] Aviles | Sensorless force estimation using a neuro-vision-based approach for robotic-assisted surgery[END_REF], Aviles et al., 2014]. The technique relies on an image classifier for estimating force quantiles during robotized intravitreal injections. An imperative requirement for machine learning algorithms to work is the huge quantity of data to train on. Currently since intravitreal injec-tions are executed manually, there is no available information on the force induced by the needle. Hence, we propose to build a patient-specific biomechanical model of the sclera to generate a very large number of force-OCT image pairs, which are then used for training any supervised machine learning method. We will show that this approach allows us to avoid the need for large data sets of real OCT images. The simulations are parameterized to match experimental results and compensate for the absence of real data. The parameterization can be carried out with different approaches. For instance, simple optimization algorithms could be used to minimize the errors between the output of the simulation and the OCT image over the set of possible parameters. However, in the view of an online estimation of the parameters, we choose to use a Reduced-Order Unscented Kalman filter [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] to estimate the stiffness of the scleras. We then train a two-layer image classifier algorithm with the generated images and their corresponding forces. This solution allows a straightforward force estimation process to take place. This two-stage process makes it possible to classify force ranges with 93% accuracy while requiring very few experimental data, as demonstrated on several ex vivo porcine eyes undergoing robotically-controlled needle insertions.

Method

In this work, we first construct a numerical model of the interesting portion of the eye to synthesize images of the deformed sclera under needle-induced forces. To be as close as possible to the actual organ behavior, we propose to use Bayesian inference as a means of identifying material parameters. Then, a data generation process takes place in order to train a neural network.

Training data generation

Biomechanical model We simulate the deformation of the sclera under needleinduced forces by modeling the eye as an elastic half-sphere subject to Dirichlet boundary conditions (see Figure 6.1a and Chapter 2 for details). Since applied forces and resulting deformations remain small, we choose to describe the stress-strain relationship as linear, using Hooke's law from Equation (2.19). Dirichlet boundary conditions are added to prevent rigid body motion of the sclera, while a constant pressure is applied to the inner domain boundary to simulate the intraocular pressure (IOP), as illustrated in Figure 6.1a. The intraocular pressure plays an important role in the apparent stiffness of the eye and its variability is well studied, as high eye pressure can be an indication for glaucoma. It is worth-mentioning that the common IOP measurement devices are influenced by the stiffness of the eye. The value measured by the tonometer is not absolute but linearly related to the Young's Force Estimation in Robotized Injections 65 modulus [START_REF] Hamilton | Young's modulus in normal corneas and the effect on applanation tonometry[END_REF]. The external force due to IOP is given by f p = S×P where S is the surface area of the eye (in m 2 ) and P the intraocular pressure, given in P a. Note that f p is normal to the surface. To simulate the needle pushing on the sclera, we apply a local force f n in a small region of interest near the virtual needle tip. The external forces of the system are then formed by the IOP, the force induced by the needle and the gravitational force f g applied in the negative y direction. Finite Element simulation We solve the equation for the constitutive model using the finite element method presented in Chapter 2. The eyeball is considered as nearly spherical and is discretized using a quadrilateral surface mesh of a half sphere of radius 12 mm, discretized using the Catmull-Clark subdivision method. The obtained quadrilateral surface is extruded conforming to the scleral thickness in order to generate almost regular hexahedra (see Figure 6.1b). Scleral thickness plays an important role in the deformation of the sclera. Under a given force, a thick sclera is less deformed than a thin one, therefore it is important to take this thickness into account.

In our simulation, we use a mesh composed of 14, 643 hexahedral elements, resulting in about 5 seconds of computation time to obtain the static solution of the deformation. All the experiments from this chapter are run on a Dell Precision laptop equipped with an intel Core i7 2.90GHz, a Quadro M1200 Graphics Processing Unit and 16Go of RAM. Since this computation has to be repeated thousands of times to generate the training data set, we take advantage of the linearity of the model and pre-compute the inverse of K to speed up the generation of the training data. This leads to a substantial computational speed-up (x10).

Force Estimation in Robotized Injections 66 Model parameterization Our finite element simulations depend on parameters of the constitutive model and the geometry of the eye: Young's modulus, Poisson's ratio, intraocular pressure and scleral thickness. It is therefore essential to understand which parameters can be assumed constant, and which ones vary (and within what range) to properly generate training data.

Data from the literature [START_REF] Olsen | Porcine sclera: thickness and surface area[END_REF] report a thickness for corneal-scleral limbus in porcine eyes ranging from 630 µm to 1030 µm. We perform thickness measurements on porcine eyes and get a variation in thickness of more than 35%. Hence, we choose to consider the thickness as a parameter in the training, and we simulated scleras with five different thicknesses: 400 µm, 500 µm, 600 µm, 700 µm and 800 µm.

The IOP is also known to vary from patient to patient, and can be measured using a tonometer. So if we were to apply our method on a patient, we could use this information to parameterize the simulations used for the training. However, this is different for the porcine eyes used in our study. The IOP decreases gradually with postmortem time [START_REF] Balci | The importance of measuring intraocular pressure using a tonometer in order to estimate the postmortem interval[END_REF] and is halved only three hours after death. The IOP we measured ranged from 1 to 4 mmHg which is very small and does not affect much the deformation of the sclera. Therefore we considered the IOP constant and equal to 2 mmHg (266 P a) for our experiments and simulations. The Poisson's ratio ν is a constant value for the sclera according to the literature [START_REF] Asejczyk-Widlicka | The elasticity and rigidity of the outer coats of the eye[END_REF] and can be set to ν = 0.45. This leads to a nearly incompressible behavior during the deformation. With this in mind, we can reasonably assume that our training data set can be generated from numerical simulations in which the thickness, Poisson's ratio and IOP are known and constant.

The value of Young's modulus E, however, is more difficult to estimate as it varies depending on the porcine breed and experimental conditions. Therefore, values from the literature are not directly applicable. On the other hand, measuring it using a "classical" experimental biomechanics approach would be burdensome. For this reason, we propose to use a Bayesian approach to estimate the value of E using observations from our OCT images. Such an approach could also be used to perform data assimilation on actual patients.

The knowledge of the elasticity parameters is key to build a good model of the sclera. The goal of this work being the estimation of the force range based on the shape of the sclera, it is therefore very important to correctly estimate E to avoid introducing errors in the force prediction. Since the exact value of E is not a priori known (only its average value based on data from the literature), we describe it as a stochastic parameter associated to a Gaussian probability density function (see Section 5.3 in Chapter 2).

Force Estimation in Robotized Injections 67

Data set generation Once the model is correctly parameterized, we generate a synthetic data set {(f nk , I k )} N k=1 of N samples where f nk is the needle-induced force for sample k and I k is the corresponding simulated image of the deformation. The forces range from 0.0 N to 0.06 N and are applied at random locations and normally to the scleral surface for each of the different thicknesses. For each sample, a 2D cross-section of the complete 3D mesh is extracted in order to simulate the OCT image I k . The simulated images have to be informative of the scleral thickness and its deformation, hence a binary image suffices for our purpose. Therefore the images are post-processed with a contour detection algorithm (findContours OpenCV function) and binarized (threshold function).

Neural network for image classification

We look into modelling a function that can estimate the range of forces applied by a needle when observing a single OCT image and the deformation within it. To do this, we make use of an artificial neural network (NN) to provide a robust and reliable function capable of estimating forces applied to the sclera (see Chapter 3 for details). In our case, the input to the model is the cross-sectional OCT image that depicts anatomical information of the sclera, from its surface to 1 millimeter below. In order to provide high-frame rate imaging, we opt to use 2D OCT crosssections as the imaging modality over 3D volumetric OCT scans that are more common but slower to acquire. In our set up, we use the low-cost solution introduced in [START_REF] Apostolopoulos | Efficient OCT volume reconstruction from slitlamp microscopes[END_REF] to image 2D OCT cross-sections.

While we are interested in estimating the force applied by the needle, we choose to categorize the applied forces into three interval ranges. This reduces the need to be sensitive to exact force values, which not only is unnecessary in this instance, but also allows us to set up our inference task as a classification problem, whereby forces are grouped into ranges of clinical relevance.

To do this, our NN consists of two fully connected hidden layers and a classifier as output layer (see Figure 6.2). The input layer has 5600 neurons (corresponding to the size of the input images 140 × 40) while the hidden layers have 600 neurons. These optimal values are found through a grid search. We use ReLu activations throughout the network and a softmax activation at the output layer. The network was implemented using Tensorflow and the Keras Python library.

Given the above simulation model, we can train our NN with virtually an infinite amount of synthetically generated data. With the objective of stopping an incumbent needle from damaging the sclera, we use our simulation model to generate OCT images with forces between 0.0 N and 0.06 N applied on the Finite Element mesh of the sclera. This range was then divided into three ranges, or classes, consisting of:

Force Estimation in Robotized Injections 68 whereby class 0 corresponds to virtually no danger to the sclera, class 1 indicates a considerable force and class 2 is a dangerous force that should trigger a needle removal signal. The complexity of our problem makes it challenging to establish a clear cut between "no danger to the sclera" and "a needle removal is advised". For this reason we introduced an intermediate class which can be seen as a gray area between the two bounding classes. If the force is classified in this gray area, it becomes the clinician's responsibility to stop or continue the insertion.

We train our NN from scratch using a gradient descent and the cross-entropy loss. A 0.8 dropout factor between layers was used to help with generalization.

Experimental results

In this section, we first describe how data acquisition on ex vivo porcine eyes was performed. Then, we present the virtually generated data set, and associated training, with an emphasis on the parameterization of the biomechanical model. Finally, Force Estimation in Robotized Injections 69 we validate the neural network on unseen real data to demonstrate the high level of accuracy of our approach at classifying contact forces.

Experimental set up

To validate the predictive accuracy of the NN on real data, we acquired a collection of M samples {(f nk , I k )} M k=1 from ex vivo porcine eyes. Five porcine eyes were obtained from the local abattoir and transported to the experiment room while kept at low temperature. Experiments began within 3 hours of death, and we ensured they were completed less than 10 hours postmortem. During the tests, the eyes were moisturized with water and fixed with super glue on a 3D-printed holder to ensure fixed boundary conditions on the lower half of the eyeball. We choose super glue as it preserves the tightness of the eyeball. We then measured the IOP with a tonometer and obtained, for all the eyes, an IOP close to 2 mmHg (about 266 P a). The intraoperative pressure is lower than usual since it decreases dramatically after death [START_REF] Balci | The importance of measuring intraocular pressure using a tonometer in order to estimate the postmortem interval[END_REF]. However, for the results presented here, we did not inject any fluid during the experiments to compensate for this low IOP.

A surgical robotic arm that was designed for high-precision drilling during robotic cochlear implantation [START_REF] Weber | Instrument flight to the inner ear[END_REF] was used to measure the force applied by a needle. The robotic arm is fitted with a six-axis force/torque sensor (Mini40, ATI) at its wrist underneath a quick release for the end effector. After mounting a 22G Fine-Ject needle (0.7 x 30 mm) at the tip of the end effector, the robotic arm guided the needle forward along its axis. The needle was positioned as close as possible to the B-scan without intersecting it to bypass the generation of shadows in the OCT image. The margin between the needle and the B-scan was considered in the simulation. The robot was programmed to move towards the center of the eye along a path normal to the sclera. Contact and insertion forces, in the direction of motion, were continuously recorded during the insertion process. We also recorded the associated OCT images by storing B-scans over time. The custom-designed OCT device uses a 840 nm ± 40 nm wavelength light source, with an A-scan rate of 50 kHz for a resulting 2D image of resolution 512 × 512 pixels with 12 bits per pixel, corresponding to an area of 15 × 4 mm 2 . In these images, the black bands surrounding the OCT focus are removed leading to images of size 140 × 40 pixels.

In Figure 6.3, a fragment of the collected data including 5 complete acquisitions is shown. To each trial corresponds one OCT image captured at a punctual time and a continuous flow of force values. The forces are filtered so that the noise is reduced. As the images were collected at a lower frequency than the forces, the corresponding forces were averaged over an interval of two seconds around the imaging time. The vertical lines correspond to the imaging times. The relaxation properties of soft tissue explain the slight decrease of the force a few seconds after the force is applied.

Force Estimation in Robotized Injections 70 Table 6.1: Distribution of samples in the three force ranges and scleral thickness for each eye.

The computed average force is used as ground truth for the following supervised learning. For each eye, one position near the corneal-scleral limbus was selected and the needle moves forward 0.5 or 1 mm in the same direction. Note that the pierced samples are excluded from the data set. Overall, 45 valid trials were performed on different eyes. The spread of each acquisition among the force classes is depicted in Table 6.1.

Identification of the Young's modulus

According to the collected porcine data, the Young's modulus E of the biomechanical model is estimated using the ROUKF. The PDE of E is transformed based on observations taken on the OCT images acquired during the experiment on an eye

Force Estimation in Robotized Injections 71

(see Figure 6.4a). We consider the corresponding ground truth force (measured with the robot's force sensor) to run the simulations. According to [START_REF] Asejczyk-Widlicka | The elasticity and rigidity of the outer coats of the eye[END_REF], we set µ 0 and σ 0 to 0.49 M P a and 0.34 M P a respectively. We estimated the Young's modulus for each eye independently based on all OCT images where observation points can be extracted. We extract 10 observations on each OCT image for a given eye (see Figure 6.4a). The finite element mesh is registered to the OCT image using rigid registration. To each observation in the OCT (yellow dots) corresponds a predicted observation attached to the model (red dots). During the assimilation, the known force is applied and the observations obtained are used by the filter correction phase to compute the Kalman gain. We employ the simplex version of the ROUKF, and since we only want to estimate one parameter, two evaluations of the model are performed in each prediction phase of the assimilation process (i.e. there are two sigma points). The evolution of µ and σ with the iterations of the ROUKF for eye 1 are shown in Figure 6.5.

After convergence, we see that the final value of µ is 0.27 M P a, as reported in first column of Table 6.2. When applied to other porcine eyes, estimated values of E vary only slightly, as seen in Table 6.2. Therefore we consider the Young's modulus E to be constant across all our experiments, with an average value of 0.25 M P a. Using this values leads to very good visual agreement between the OCT images and the simulations for the three force ranges for all the eyes (see Figures 6.6 and 6.7) thus validating the assumption of linearity.

eye 1 eye 2 eye 3 eye 4 µ in MPa 0.27 0.22 0.24 0.29 σ in MPa 4.7e -3 6.7e -3 3.9e -3 4.9e -3 Table 6.2: Mean and standard deviation of the estimated Young's modulus in 4 different porcine eyes after performing a data assimilation with a ROUKF. The fifth eye is not included since the observations required in the assimilation process are not consistent across the OCT acquisitions for this eye.

Training data generation

Data set generation To run our simulations we used the Simulation Open Framework Architecture1 . We generated 3200 images of deformed scleras undergoing normal forces going from 0.0 N to 0.045 N for the smallest thickness. For each of the other four thicknesses, we generated 4000 images of deformed scleras where the forces vary from 0.0 N to 0.06 N at different random locations. For each thickness, the simulation took approximately half an hour. In Figure 6.7 are shown different examples of the output of the simulation (bottom) matching the OCT images (top). Overall, we created a data set having 19200 synthetic images within 2 hours and a half (see Figure 6.8b). In Figure 6.7 below we show some synthetic images.

Neural network training

The generated labelled data set is split such that 90% of the images is used for training and the remaining 10% is left for validation. All in all, 17280 images are used to train the artificial neural network that is validated on the other 1920 images. Figure 6.8a depicts the accuracy and the loss of the neural network on both training and validation data sets over each epoch. The validation accuracy curve displays 100% accuracy when classifying unseen synthetic images. This curve is above the training accuracy curve probably because of the high dropout applied during the training.

Validation on unseen real data

Our work aims at classifying force levels using only OCT data of the deformed scleras as input. All the OCT images collected during the experiments are filtered to obtain inputs similar to the synthetic ones (see Figures 10.2, 6.9b and 6.8b).

Once the OCT images have been processed to look like the simulated ones, we can use them as input to the NN and perform force predictions. For each OCT acquisition, the force measured by the robot is converted into a class label (0, 1 or 2) and is taken as the ground truth (target class). A class with label 0 corresponds to a minimal force for which there is no risk of damaging the sclera. A class with label 1 corresponds to forces ranging from 0.005N to 0.03N and indicates that the sclera is being considerably deformed. A force class with label 2 means that the sclera is potentially being damaged and that a withdrawal of the needle is advised.

The performance of the classification is reported in the confusion matrix in Table 6.3. Each row of the table gives the instances in a target class, and each column gives the instances in a predicted class. For each category, we highlighted the correct decisions in red and show that the overall accuracy of the classifier is very high (93%). For all the experimental data set, the lowest score of the NN was obtained for the force class 0 with 71% accuracy. For target class 2, the precision reached 100%. In Figure 6.10 are shown the raw force values (measured with the robot's force sensor) and the force range thresholds. The plot shows that the forces measured are uniformly spread throughout the force ranges. Note that the three misclassified samples are close to the upper bound of the force range and they are overestimated. We believe that misclassification of contact forces in the lowest ranges has a limited impact since the risk of damaging the sclera with such forces is almost null. On the other hand, it is essential that the forces of range 2 are predicted correctly, which appears to be the case on our (limited) data set.

Conclusion

This chapter introduces a method for improving the safety features of robotized intravitreal injections. We show that our vision-based method, which combines numerical simulation and neural networks, can accurately predict the level of force applied by a needle, using only 2D OCT images of the scleral deformation. By being real-time, this classification can lead to an immediate withdrawal of the needle once it reaches a certain alarm threshold. To cope with the issue of performing a large number of experiments to populate our training data sets, the NN is trained on synthetic images generated from simulations of sclera deformations. To automatically parameterize the simulations from the experimental data, we use a Kalman filter which performs data assimilation using the sequence of OCT images.

This work was carried out at the beginning of the PhD. At that time, we were inexperienced in deep learning algorithms so we selected the simplest network architecture for solving our classification problem. With the current experience, if we were to do this work again we would probably use a CNN instead. Indeed, CNNs are

Force Estimation in Robotized Injections 76 aimed to deal with images and they are usually faster to train than fully connected neural networks (due to the smaller number of learnable parameters). Moreover, CNNs are less sensitive to image framing and scaling and would be better suited if we were to add real OCT images in our data sets.

Moreover, in order to use our method in clinical practice (e.g. in vivo eyes) further simulation efforts should be performed. For instance, anatomically reasonable boundary conditions should be considered such as muscle fixation and eye movement should equally be included. It also seems essential to include the intraocular pressure as an additional input in the data set. Indeed, IOP is known to play a role in the deformation of the sclera and cornea of patients. To perform the prediction, a measure of the IOP would simply need to be performed on the patient using a tonometer. The reason we have not done it yet comes from the inaccuracy of the IOP measurement which is influenced by the stiffness of the eye. This coupling makes it more challenging to parameterize our simulations. Another improvement of the method would be to use adaptive force thresholds depending on the scleral thickness and stiffness as for thinner or softer scleras, the puncture force is lower. The simulated images could also be improved to match actual surgical scenarios better. In particular, adding the shadows generated by the needle in the OCT image would be an essential feature.

We also propose to address the sensitivity of the NN predictions to image framing and scaling by randomly cropping each simulated image, and augmenting the data set with these additional images. Note that in our data generation process, the optical interference phenomenon happening in OCT is not simulated. Instead, the synthetic and real images are binarized. This skeletonization is a way of filtering the complicated aspects of the OCT images that are a priori non-informative for our purpose. Training a network directly from the OCT images without the skeletonization would lead to a significantly more complex problem requiring the simulation of light propagation and the use of Convolutional Neural Networks.

Finally, the objective of this work was to estimate a force range, but it would also be pertinent to estimate the location and the angle of the force as a slight error in the latter might damage the retina or the eye lens. In this context we could augment the data set by including various needle insertion angles and locations, and use a regressor to predict the complete force vector.

Force Estimation in Robotized Injections 77

Contributions of this chapter:

• Simple fully connected neural network for force classification Chapter 7

The U-Mesh for Augmented Hepatic Surgery

Introduction

In Chapter 4, we have presented the performance of the U-Mesh framework on simple geometries, matching the required grid-like structure for the input to the network. In this chapter, we will see how to extend the method to any kind of geometry. In particular, we will first show the ability of the U-Mesh in predicting the deformations of a liver geometry. Second, we will apply it to augmented hepatic surgery by registering the preoperative images to the intraoperative ones.

In liver surgery, it is essential to locate the internal structures (such as tumors and blood vessels) that need to be preserved for the post-operative regeneration of the liver tissue. While the initial position of these structures is known from the preoperative images, their actual position during surgery is often hidden or uncertain. To guide the surgeon, augmented reality techniques are used to enrich visual information through fusion of intraoperative images and a preoperative 3D model of the patient's anatomy. This is usually done by overlaying a virtual representation of the liver built from preoperative images over intraoperative images or through augmented reality glasses. However, surgical manipulations and interactions with the surrounding anatomy can induce significant deformations to the patient's liver. As a consequence, the virtual model of the liver has to account for non-rigid transformations and produce its deformed state in real-time, which is difficult given the complexity of the physical systems needed for accurate biomechanical modeling. Therefore, in order to build an augmented view of a liver during surgery we need to perform an elastic registration of the preoperative model to the intraoperative images acquired with a 3D imaging device (see Figure 7.1). While in minimally-invasive surgeries a laparoscopic camera can be used to acquire a video of the abdominal cavity, in open surgeries an RGB-D sensor can capture the surface deformation of the tissues. From such images, a partial point cloud of the liver surface can be extracted using one of the methods listed in [START_REF] Petit | Environment-aware non-rigid registration in surgery using physics-based simulation[END_REF]. In this chapter we will see how the U-Mesh can produce an estimation of the complete displacement field given such a point cloud. This chapter is divided in two main segments. To start with, the U-Net learns to predict the displacement field of a virtual liver given an input contact force. In a follow-up, the U-Net is used in an augmented reality context, where the full volumetric displacement field needs to be estimated from a partial surface deformation.

In both cases the network is trained with FEM-generated data and the procedure is very similar to the one presented in Chapter 4 except from the meshing procedure.

Displacement field on regular grids

To train our network, we need to generate many samples consisting of a volumetric output displacement field U v , given an input constraint C (contact force or surface displacement). The simulation that generates this data needs to be as accurate as possible, yet computationally efficient to make it possible to generate data and train the network in less than a day (i.e. in the shortest amount of time available between preoperative data acquisition and surgery). As seen in Section 4.2, the input to the U-Net must be a regular grid. In the following we will present two options to fulfill this requirement.

First, the nodes of the finite-element mesh can be mapped into the nodes of the input grid. That way, the FE mesh and the U-Net input grid are independent and any type of elements can be used to mesh the object (tetrahedral elements, non regular hexahedral elements, linear or quadratic elements etc). In order to mesh complex geometries (in particular to handle the boundaries of an object), using 4-node tetrahedral elements can be more convenient than using 8-node hexahedral elements [START_REF] Wang | Back to elementstetrahedra vs. hexahedra[END_REF]. Therefore we propose to start discretizing a liver geometry with tetrahedral elements and map a regular grid onto the tetrahedral mesh. The regular grid follows the deformation of the FEM mesh (see Figure 7.

4b). Only the nodes

The U-Mesh for Augmented Hepatic Surgery 81 of the grid that are inside the liver volume are mapped and the outer nodes are zero-valued. Indeed, the mapping is only exact for inner grid nodes since we can use the FE shape functions which is not the case when solving the grid nodal values outside the geometry (ill-posed problem). The objective being to find the location of the internal structures, this simplification is not a limitation.

The second option consists of building the FE mesh such that it directly matches the U-Net grid. Here we will rely on 8-node hexahedral elements known for their better convergence, and lock-free behavior, especially with close to incompressible materials (such as the liver) and strong shear stresses [START_REF] Benzley | A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis[END_REF]. In addition, we combine this choice of elements with an Immersed-boundary method (IBM), which allows us to create an hexahedral mesh directly from the segmented CT or MRI image (see Figure 7.2). The surface mesh is embedded into a sparse hexahedral grid (Figure 7.5a), that is in turn embedded into a regular grid (Figure 7.5b) which matches the structure of the first layer of our U-Net. The sparse grid consists of rectangular cuboid cells. Cells that are overlapping the domain boundary are kept, therefore approximating the exact shape and volume of the object. The choice of discretizing the initial domain with a sparse and regular grid brings several benefits. The usual trilinear interpolation functions of an 8-nodes hexahedral element are reduced to a linear mapping, and, similar to 4-nodes tetrahedral elements, its jacobian remains constant inside the element. However, using such discretization requires particular care of the boundary elements. Volume integration of the displacement field inside these partially filled cells is carried on by recursively subdividing the cell into 8 sub-cells. The stiffness matrix inside a boundary cell is then accumulated from its sub-cells using the linear mapping, as in [START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF]. Since we are using a fine mesh of the domain, only one level of subdivision is enough to obtain an accurate approximation of the volume integral. Note that the finer the grid, the

The U-Mesh for Augmented Hepatic Surgery 82 smaller is the difference between the exact volume and the one represented by the sparse grid.

7.3

The U-Mesh applied to a liver

Results on a synthetic liver

A surface mesh is obtained from a preoperative CT scan of a human liver. The length of the liver is 0.2 m. We will first mesh the liver using T4 elements and then using H8 elements. In both cases the Young's modulus will be set to 5, 000 P a and the Poisson's ratio to 0.48.

Liver meshed with tetrahedral elements

In this paragraph the liver geometry is discretized into 4859 tetrahedral elements (1059 nodes in total). A 16 × 15 × 16 regular grid is mapped onto the tetrahedral mesh and follows the deformation of the FEM mesh (see Figure 7.4b). Dirichlet boundary conditions are then added by fixing some nodes in the area separating the two lobes to mimic the effect of the vascular tree and of the falciform ligament of the liver [START_REF] Abdel-Misih | Liver anatomy[END_REF].

For the training data generation, normal forces of random magnitudes are computed on the liver surface and mapped to the grid nodes using barycentric coordinates. Only one force is applied at each time step on a small region of the surface. We decided to limit the size of the data set to fit the time requirements of a clinical routine where sometimes only a few hours are available between the preoperative CT scans and the surgery. Hence a data set of only 2, 000 samples is generated in 135 min. N = 1, 600 samples are used to train the network in 149 min and M = 400 samples are left for validation.

Once the network is trained, we make predictions over the validation set. The results are reported on Table 7.1. The sample with maximal error is shown in Figure 7.4a. The outputs are predicted in only 3 ms. In Figure 7.3 are shown some samples of U-Mesh-deformed livers and their corresponding relative errors computed at one of the lobe tips. The output of U-Mesh is in green whereas the reference solution is in red. ponding relative errors computed at one of the lobe tips. Furthermore, the slope of the regression in Figure 7.8 shows that the increase of the error with the deformation amplitude is also controlled for this scenario.

Results in augmented surgery

To apply the U-Mesh in order to generate an augmented view of an hepatic surgery, it needs to learn to predict full volumetric displacements from partial surface point clouds that give information about the position of some points of the surface of the liver. These positions can be translated as prescribed constraints. We choose to use a FEM combined with an Immersed-boundary method to generate the data sets.

We can assume that during an open liver surgery, half of the surface of the liver is visible to the camera. As depicted in Figure 7.10a, 100 points are uniformly sampled in the visible part of the surface to mimic a point cloud. Then, 100 simultaneous forces of random magnitude and direction are applied to these points in order to generate nearly random displacements. The training data set consists of pairs of (U s , U v ) where the input to the network U s corresponds to the surface point cloud mapped onto the regular grid. For the same reasons stated in previous section, we limited the size of the data set to 2, 000 samples (N = 1, 600 for training and M = 400 for testing). The reader may refer to Figure 7.9 for examples of the generated

The U-Mesh for Augmented Hepatic Surgery 86 deformations. It is worth mentioning that no patient-specific parameterization of the biomechanical model is required since for homogeneous materials, the relation between the surface and the volumetric displacements is independent of the stiffness of the object [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF], and only depends on the Poisson's ratio (set to 0.49 as soft tissues can generally be described as incompressible).

Figure 7.9: Examples of generated deformations to train the network.

Ex vivo liver Once the network is trained, we assess our approach on ex vivo human liver data, on which ten markers are embedded to compute target registration errors (TRE). Note that the network is trained on the considered ex vivo liver geometry. During the experiments, surface data is obtained with an RGB-D sensor and ground truth data acquired at different stages of deformation using a CT scan. Before computing the displacement field, the point cloud needs to be cropped to the portion of the surgical image that contains the liver. This is done by segmenting the associated color image, similarly as in [START_REF] Petit | Real-time tracking of 3D elastic objects with an RGB-D sensor[END_REF]. The RGB-D point clouds can be interpolated onto the regular grid to obtain per-node displacements on the surface and can be given as input to the network that in turn predicts the volumetric displacement fields.

The marker predicted positions are compared to our ex vivo ground truth by computing TREs (see Figure 7.10b). The average TRE at the 10 markers is of 7.5 mm with a maximal value of 10.5 mm. The solution of the Saint Venant-Kirchhoff model, used to train our model, gives nearly the same error (which was expected) but for a computation time of 1550 ms.

In vivo liver With a PyTorch implementation of the U-Net running on a GeForce 1080 Ti, the network can predict the volumetric deformation of the liver in only 3 ms.

We can then apply this prediction each time the RGB-D camera generates a point cloud. The RGB-D point cloud is then interpolated onto the grid to obtain per-node displacements on the surface (i.e. U s ). Given this input, the network predicts the volumetric deformation, and the next point cloud can be processed. Each new RGB-D point cloud can be fed to the network, thus generating a continuous visualization of the internal structures of the organ. We show some frames of the obtained augmented reality view in Figure 7.11. The quality of the registration seems acceptable in the

The U-Mesh for Augmented Hepatic Surgery 87 frames on top. However, in the bottom frames, where occlusions appear, the quality is far from being satisfactory. We observe very similar behaviors when using the FEM (both with the Saint-Venant Kirchhoff or the co-rotational model). This means that the issue probably comes from the iterative closest point (ICP) algorithm, used for matching the points from the RGB-D point cloud to the surface model. As a matter of fact, our implementation of the ICP does not necessarily handle occlusions and when very large deformations happen between two consecutive frames, the matching fails. Another possible source of error is the model parameterization. Indeed, the liver model is only fixed at nodes that are close to the falciform ligament (whose location is approximate) whereas some other areas of the vascular tree might be stiff enough to consider them as fixed. In order to overcome this source of inaccuracy, in Chapter 9 we propose to estimate such fixations intraoperatively using Kalman filters.

Conclusion

In this chapter we have shown the ability of the U-Mesh framework in learning complex elastic deformations of any type of geometry. In particular, the U-Mesh learns the mechanical behavior of a liver and generates its deformed state about 500x faster than a reference FE solution. Since the network takes as input a regular sparse grid where displacements are imposed, we have shown an efficient FE Immersedboundary method based on the same hexahedral discretization from which thousands of deformed configurations are generated to train the network. Complex and accurate

The U-Mesh for Augmented Hepatic Surgery 88 Driven by surface displacement data, it makes this approach an ideal candidate for providing augmented reality during surgery.

The main limitation of our approach relies in the encoding of the RGB-D point cloud into the U-Net input grid. Each point of the acquired point cloud has to be associated to one (or possibly more) nodes on the hexahedral grid. The problem is that a generic point p 2 from the intraoperative RGB-D point cloud corresponds to a point p 1 on the preoperative liver surface (found using ICP), which is not necessarily a node of the grid (see figures 7.12a and 7.12b). However, our method requires that p 2 coordinates are mapped into grid coordinates. This problem does not have a unique solution. For example, we could set the value of p 2 to the grid node that is closest to p 1 (Figure 7.12a) or spread its value to the 8 nodes of the cell containing p 1 (Figure 7.12b). These two options introduce an approximation error that can be minimized by increasing the grid resolution.

In order to avoid such error, we decided to use the FE shape functions to obtain the nodal grid values. In particular, we computed the FE solution with the Saint-VenantKirchhoff model on the hexahedral grid by imposing the acquired surface displacement. Then, we stored the positions of the grid nodes associated to those cells containing the points of the RGB-D point-cloud to build a virtual point cloud that is in turn fed to the network (see Figure 7.12c). The main limitation of this approach is that we are also giving some volumetric information as input to the network (n 3 , n 4 , n 6 and n 7 are grid nodes inside the volume), which is not optimal. An alternative solution worth trying is depicted in Figure 7.12d. It consists of

The U-Mesh for Augmented Hepatic Surgery 89

encoding each point of the point cloud as a distance field from each of the grid nodes, similarly as works in [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration using a Data-Driven Biomechanical Model[END_REF]. The U-Mesh for Augmented Hepatic Surgery 90

Contributions of this chapter:

• The U-Mesh applied to non regular geometries • Immersed-boundary method to encode the displacement fields in a regular grid structure given as input to the network • The U-Mesh for augmented hepatic surgery with a simple set-up (one single RGB-D camera) • Publication: Mendizabal, A., Brunet, J.N., [START_REF] Brunet | Physics-based Deep Neural Network for Augmented Reality during Liver Surgery[END_REF]. Physics-based deep neural network for augmented reality during liver surgery. International Conference on Medical image computing and computer-assisted intervention, pp. 137-145

The U-Mesh for Augmented Hepatic Surgery 91

The U-Mesh for Augmented Hepatic Surgery 92

Chapter 8

The U-Mesh for Ultrasound-guided Breast Biopsy

Introduction and related work

To demonstrate the generality of our approach, we propose to apply it to a different organ and a different intervention: US-guided breast biopsy. Breast biopsy is the preferred technique to evaluate the malignancy of screening-detected suspicious lesions. To direct the needle towards the target, biopsy procedures are performed under image guidance, normally done with ultrasound (US) probes due to their ability to provide real-time visualization of both the needle and the internal structures [START_REF] O'flynn | Image-guided breast biopsy: state-of-the-art[END_REF]. However, proper needle placement with US remains a challenging task. First, malignant lesions cannot always be adequately visualized due to the poor image contrast of US. Furthermore, navigation towards complex 3D lesion geometries is commonly achieved using 2D freehand US (FUS) systems, which provide information in a lower-dimensional space [START_REF] Krücker | Clinical utility of real-time fusion guidance for biopsy and ablation[END_REF]. Since highly sensitive preoperative images (such as MRI or CT) can provide accurate positions of the lesions, finding a method to update these positions from real-time US images during an intervention would highly benefit current biopsy procedures. Several commercial and research platforms have implemented image fusion techniques that align preoperative and intraoperative data, exploiting rigid or affine registration methods [START_REF] Guo | Ultrasound imaging technologies for breast cancer detection and management: A review[END_REF]. However, when dealing with breast anatomy, large deformations arise due to compression forces applied by the US probe. To provide accurate probetissue coupling and acceptable image quality, an appropriate alignment procedure of the preoperative and US data is required.

Although FE models have been successfully employed for multimodal breast image registration, they have never been applied to registration between preoperative data and intraoperative US, due to difficulties in providing a prediction within real-time constraints [START_REF] Hipwell | A review of biomechanically informed breast image registration[END_REF]. This is especially true when considering large, non-linear deformations which involve hyperelastic objects, as it is the case for the breast.

As mentioned in Chapter 4, various techniques have been proposed to simplify the computational complexity of FEM in order to meet real-time compliance. The most optimized approach used to model breast biomechanics is the one proposed by [START_REF] Han | A nonlinear biomechanical model based registration method for aligning prone and supine MR breast images[END_REF], which relies on GPU-based TLED formulation. Despite the significant simulation speedup achieved, solving the FE system took around 30 s, which is still not compatible with real-time. Modelling methods that do not rely on continuum mechanics laws have also been used to approximate soft tissues behavior. Among these, the position-based dynamics (PBD) approach has been used to predict breast lesions displacement due to US probe pressure in real-time, providing comparable accuracy with FE models [START_REF] Tagliabue | Position-based modeling of lesion displacement in ultrasound-guided breast biopsy[END_REF]. However, not being based on real mechanical properties, such model requires an initial optimization of simulation parameters to obtain a realistic description of the deformation.

In this chapter we propose to use the U-Mesh framework to estimate breast tissue behavior in US-guided biopsy. Using FE simulations for model training in the context of MRI-US deformable image registration has already been proposed in [Hu et al., 2016], where the authors build a statistical model of prostate motion which can account for different properties and boundary conditions. In the case of the breast, the potentiality of employing machine learning techniques has been already shown in [START_REF] Martínez | [END_REF], where several tree-based methods have been employed to estimate breast deformation due to compression between biopsy plates. These methods have been trained on 10 different patient geometries with a very specific FE simulation, where the upper plate is displaced vertically towards the lower one.

In this chapter, we propose to train the U-Mesh to predict the deformation of internal breast tissues starting from the acquired surface displacements induced by the US probe. Our network can be seen as a patient-specific model. We train it on a single patient geometry before surgery, with a relatively small amount of training data. However, in contrast to the work of [START_REF] Martínez | [END_REF], FE simulations that compose the training set are generated with several random input displacements, making our approach able to generalize to different probe positions and compression extents.

Similarly to Chapter 7, we propose to use an Immersed-boundary method for generating patient-specific simulations. Results presented in Section 8.3 show the efficiency of the method when applied to both synthetic and phantom data. The contribution of this chapter consists of a novel method to generate a real-time capable soft tissue model to improve target visualization during needle-based procedures. The

The U-Mesh for Ultrasound-guided Breast Biopsy 94 position of lesions identified beforehand on preoperative images can be updated from intraoperative ultrasound data and visualized by the surgeon in real-time.

Method

This chapter consists of a data-driven method to estimate in real-time the displacement of the breast internal structures due to probe pressure during US scanning.

In our pipeline, we assume to have a patient-specific geometric model of the breast, obtained from preoperative imaging such as MRI, and to know the position and orientation of the US probe at each time, thanks to a spatial tracking system. If the tracking coordinate system and the coordinate system of preoperative imaging are registered, knowledge about the 3D pose and the geometry of the US probe directly allows to identify the contact surface between the breast and the probe. Since the US probe is represented as a rigid body, we can reasonably assume that when the anatomy is deformed by the probe during the image acquisition process, points on the breast surface below the US probe will be displaced to the same exact extent as the probe itself. As a consequence, our method can predict the displacements of all the points within the anatomy given as input the displacement of the surface nodes in contact with the US probe. The decision of relying on surface displacement inferred from the spatial tracking of the US probe instead of directly tracking surface deformations (through, for example, an RGB-D camera) was taken from the fact that probe-induced deformations are large but local, and the probe itself would occlude most of the deformed surface to the sensor, thus preventing an accurate estimation of the contact surface displacements.

Simulation of breast tissue

The training data set consists of pairs of (U s , U v ) where U s is the input partial surface displacement and U v is the volumetric displacement field. Even though the data generation process takes place in an offline phase, in order to generate enough training data with FE simulations within clinically acceptable times (the intervention can be performed on the day after preoperative scan is acquired), it is important to have simulations that are both accurate and computationally efficient.

We consider the boundary value problem of computing the deformation on a domain Ω under both Dirichlet and Neumann boundary conditions. Let Γ be the boundary of Ω (in our case, Γ corresponds to breast external surface, while Ω represents the entire breast volume). We assume that Dirichlet boundary conditions are applied to Γ D and are a priori known, whereas Neumann boundary conditions are applied to Γ N , a subset of Γ that represents probe-tissue contact area and changes depending on current US probe position. We choose to discretize the domain into 8-node hexahedral elements not only for their good convergence properties and lock-free behavior, but also because it is the required structure for the input to the network. To do that, the 3D breast geometry is embedded in a regular grid of hexahedral elements (see Figure 8.1) and we use an Immersed-boundary method to correctly approximate the volume of the object in the FE method computations (see Section 7.2 in Chapter 7 for details).

Training data generation

The input to the network corresponds to the displacement U s of the points belonging to the breast-probe contact area. The punctual displacements are spread to the nodes of the surrounding cuboid cell through a barycentric mapping and the corresponding volume displacement U v is obtained by the previously explained FE approach in response to U s . The data used to train the network must be representative of the application scenario and must allow the network to extract the pertinent features of the tissue behavior. In order to train our model to estimate breast volume deformation in response to pressure imposed with the US probe, we simulate several random probe-induced deformations using the following strategy:

• Select a random node p in the breast surface

• Select an oriented bounding box A centered in point p and normal to the breast surface, whose dimensions match those of the US probe lower surface, which represents current probe-tissue contact area

• Select all the surface points P falling within the box A

• Select as force direction d the normal to the surface at point p plus a random angle α (α ∈ -π 4 , π 4 )

The U-Mesh for Ultrasound-guided Breast Biopsy 96

• Apply the same force f of random magnitude (|f | ∈ [0.0, 0.8]) along direction d to the P selected points simultaneously

• Store the displacement at the set of points P (input to the network) and the displacement of all the points in the volume (output to the network)

• Repeat the procedure until N + M samples are generated

The choice of applying force f allowing some angle deviation from normal direction enables us to include in our dataset samples where the probe compression is not precisely normal to the surface, as it can be the case in freehand US acquisitions.

The maximal force magnitude (e.g. 0.8 N ) is set such that the amount of maximal deformation reproduced in the training dataset never exceeds too much that observed in real clinical settings. The described strategy is used to generate the set

{(U s n , U v n )} N n=1
of N samples which is used to train the network, and the set

{(U s n , U v n )} M n=1
of M samples which is left for validation. The training dataset is generated with the SOFA framework [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF] on a laptop equipped with an Intel i7-8750H processor and 16GB RAM.

Experiments and Results

The network presented in this work is used to predict US probe-induced deformations of a realistic multi-modality breast phantom (Model 073; CIRS, Norfolk, VA, USA). The 3D geometry model of the phantom surface and 10 inner lesions (diameter of 5-10mm) is obtained by segmenting the corresponding CT image, relying on ITK-SNAP and MeshLab frameworks [START_REF] Yushkevich | User-Guided 3D Active Contour Segmentation of Anatomical Structures: Significantly Improved Efficiency and Reliability[END_REF], Cignoni et al., 2008]. A Freehand Ultrasound System (FUS) based on a Telemed MicrUs US device (Telemed, Vilnius, Lithuania) equipped with a linear probe (model L12-5N40) is used to acquire US images of the 10 segmented lesions. The dimension of the probe surface is (5x1cm). For each lesion, we acquire US images in correspondence of four different input deformations. The MicronTracker Hx40 (ClaronNav, Toronto, Canada) optical tracking system is used to track US probe in space (Figure 8.2a). The overall probe spatial calibration error is below 1mm (±0.7147), estimated through the PLUS toolkit [START_REF] Lasso | PLUS: open-source toolkit for ultrasound-guided intervention systems[END_REF]. Landmark-based rigid registration is performed to refer the CT-extracted 3D model, the US probe and the US images to the same common coordinate system, exploiting 3D Slicer functionalities [START_REF] Fedorov | 3D Slicer as an image computing platform for the Quantitative Imaging Network[END_REF]. The registration process does not only enable us to extract the breast-probe contact area, as described in Section 8.2, but also to know in real-time the 3D position of any point belonging to the US image. In this way, it is possible to refer lesions position extracted from US images to the 3D space.

The U-Mesh for Ultrasound-guided Breast Biopsy 97

Predictions on synthetic data

Elastic properties of the physics model used to generate training data are set in accordance with the values estimated in [START_REF] Visentin | Iterative simulations to estimate the elastic properties from a series of mri images followed by mri-us validation[END_REF] for the same breast phantom considered in this study. However, as we are imposing surface displacements, the values of the elasticity parameters do not affect the displacement field inside the simulated volume as long as the ratio of the different stiffness values is maintained [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF], thus making the method reliable for any patient specificity. Dirichlet boundary conditions are imposed by constraining the motion of all the nodes belonging to the lowest phantom surface.

Using the method described in Sections 8.2.1 and 8.2.2, we discretized the breast phantom into 2174 hexahedral elements and we simulated several probe-induced displacements. Overall we generated N = 800 samples for training and M = 200 samples for testing. The U-Net is trained in a GeForce GTX 1080 Ti using a batch size of 4, 100000 iterations and the Adam optimizer. We used a Pytorch implementation of the U-Net. To assess the learning capability of the network, we perform a statistical analysis of the mean norm error e over the testing data set following Equation (4.3). We compute the average e, standard deviation σ(e) and maximal value of such norm over the testing data set. The obtained results are shown in Table 8.1. The maximal error is of only 0.266 mm and corresponds to the sample shown in Figure 8.3b. The most striking result is the small computation time required to make the predictions: only 3.14 ± 0.56 ms. In contrast, the FE method takes on average 407.7 ± 64 ms to produce the solution. Obviously, the resolution of the FE mesh could be reduced to accelerate the computations but at the cost of an accuracy loss.

Predictions on phantom data

In our experiments, we consider one lesion at a time and we reposition the US probe on the surface of the breast such that the lesion considered is visible on the US image. In order to validate our model, we manually extract lesions position from US image acquired at rest (i.e., without applying any deformation, when the probe is only slightly touching the surface) and we consider it as a landmark to track. We then impose four deformations of increasing extent for each lesion, and we compare the U-Net-predicted displacement with real displacements extracted from US images. The comparison is performed computing target registration error (TRE) between the predicted position of the lesion and its ground-truth position. The performance of our method is compared to that of the FE model used for data generation. In Table 8.2 are shown the target registration errors for each phantom lesion with respect to the applied deformation. The input deformations are classified into five ranges based on the probe displacements. Displacement ranges indicated as D15, D20 and D25 have a fixed length of 5 mm each and are centered at 15, 20 and 25 mm respectively. D10 and D30 contain the extreme cases under 12.5 mm or above 27.5 mm.

Values in Table 8.2 highlight that the average TRE for all the tumors and for all

The U-Mesh for Ultrasound-guided Breast Biopsy 99 the deformations is smaller than 6.194 mm which is comparable to the maximum value obtained with the FE method (6.080 mm). The average error increases with the deformation range just like in the FE method. There is no significant difference between the values of the two tables, meaning that in terms of accuracy, our method is comparable to the data generation method used to train it. In order to compute each deformation, the FE method needs about 407.7 ms whereas the U-Net predicts the deformation in only 3 ms.

The U-Mesh for Ultrasound-guided Breast Biopsy 100

Conclusion

In this chapter we have proposed to use the U-Mesh framework to learn the deformable behavior of the breast from numerical simulations based on the finite-element method, in order to bypass the high computational cost of the FEM. Our approach represents an interface between precise biomechanical FE modeling (not capable of real-time) and clinical applications requiring both high accuracy and very high speed. We have shown that our framework allows for extremely fast predictions of US probe-induced displacements of the breast during US scanning, achieving comparable accuracy to other existing methods. Therefore, it has the potential to be employed to update in real-time the estimated position of breast lesions identified on a preoperative scan on US images, enabling continuous visualization of the biopsy target, even when sonography fails to render it.

Although the FE model used to train our network does not perform in real-time, its prediction delay of less than 1 s might be considered already acceptable for our specific application. However, such good computational performance is achieved since in this preliminary evaluation we use a very simplistic model, that does not account for heterogeneity or complex boundary conditions happening in clinical cases. Usage of a more complex FE model will certainly cause an increase of computation load. On the contrary, an important feature of our approach is that the prediction time remains close to 3 ms regardless of the grid resolution and of the biomechanical model used for the data generation process. This means that increasing the complexity of the model used to generate the data set will not affect the prediction speed. Moreover, our pipeline allows the method to be insensitive to patient specific elastic properties as it imposes surface displacements. It is worth noting that for inhomogeneous objects, the displacement field still depends on the ratio of the different stiffnesses [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF]. Another advantage of our method is the easy meshing process. Any geometry can be embedded in a sparse grid and through the use of Immersed boundary simulations the deformations are correctly estimated.

The main limitation of our method remains the training process, which is burdensome and has to be repeated for every new geometry or application. However, we have shown that a limited amount of training data can be sufficient to train a U-Net such that it obtains accurate prediction within clinically acceptable times.

Chapter 9

Towards Patient-Specific Networks

Introduction

In Chapter 7, we have shown the potential of the U-Mesh for augmented hepatic surgery. However, we made the assumption that the boundary conditions (BCs), the material properties and even the constitutive law were known a priori. As seen in Chapter 5, this assumption can be reductive. The location and the elastic properties of the BCs are patient-specific and are not visible on preoperative images. Therefore, they cannot be identified in advance, which is necessary for the training data generation process. On the other hand, the identification of the elasticity parameters was not crucial in the previously considered scenario as we were dealing with surface displacements as input [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF]. Nevertheless, if we were to consider a scenario where the input to the network is a traction force, the material parameters do have an impact on the deformation and need to be correctly estimated. Furthermore, the material properties and BCs are relative to a constitutive law that can itself be inappropriate. Therefore, it is important to see how the U-Mesh performs in such situations and how it can be adapted.

This chapter is a proof of concept that explains how the U-Mesh could be integrated in a real clinical routine for patient-specific predictions of the displacement field. It consists of a set of possible research directions with some preliminary results or bibliographic references supporting their feasibility. The variety of parameters intervening in a simulation, may require several strategies in order to adapt the networks to patient-specificity. We will propose here two major solutions relying on two opposite paradigms.

The first one consists of training the network with all possible combinations of pertinent parameters, namely with exponentially large amounts of training data. Such parameters must be sampled and explicitly encoded in the input to the network (in particular the domain of possible geometries, elasticity parameters and even BCs).

The second option, more subtle, relies on transfer learning techniques where the net-work's weights are updated based on patient-specific data acquired intraoperatively. As a matter of fact, when more primary components of the simulation (such as the constitutive law) must be revised, it seems pertinent to explore this alternative as the learned underlying function must be modified. In an ideal scenario, we could use real patient data to retrain our model. An interesting point of this solution is that the parameters do not need to be explicitly identified as they are encoded in the data, meaning that the material properties, the BCs and the constitutive model can simultaneously differ from the ones learned preoperatively. However, as seen previously, collecting volumetric data on a patient without exposing him or her to radiation is in practice very difficult. To overcome this issue, we can generate new synthetic data (where the simulation has been parametrized using intraoperative observations) and retrain the neural network with these new simulations. The BCs and the material properties are relative to a constitutive law (that is difficult to characterize intraoperatively), meaning that their correct estimation can overcome the inaccuracy of the constitutive law (if any). Therefore, we will demonstrate how transfer learning can be used to adapt the network to patient-specific material parameters and BCs.

Once the network is trained as described in Chapter 7, we propose to estimate the boundary conditions and the relevant material parameters intraoperatively in order to update the network parameters with patient-specific FE simulations. To do that, we first use an image-driven stochastic assimilation method based on Kalman filtering to identify the BCs on the one hand, and the elasticity parameters on the other hand. Then, the pre-trained model is adapted to the patient specific properties through transfer learning. In this chapter we only show some preliminary results on a cantilever beam to showcase the potential of transfer learning in our applications.

This chapter is divided in three main segments. First of all we will see how to use Kalman filters to identify the BCs on the one hand, and the elasticity parameters on the other hand. Second of all, we will briefly explain how the identified parameters and other patient-specific properties (such as the geometry) can be encoded in the input to a generic network, in order to make patient-specific predictions. Lastly, we will demonstrate how the parameters identified in the first section can be used to update the network's weights using transfer learning.

Stochastic identification of patient-specific properties

We propose to use the ROUKF presented in Section 5.3, to estimate the value of Young's modulus and the BCs of a liver using observations of the target model. In this section, the ROUKF algorithm is first used to estimate the Young's modulus of In this paragraph, we aim at estimating the value of the Young's modulus of the liver using the ROUKF. This estimation is done using synthetic data, but a similar process can be followed for real data. We build a biomechanical model of the liver, with fixed boundary conditions (red points in Figure 9.1a) to mimic the effect of the falciform ligament and of the vascular tree. A force of fixed magnitude and varying amplitude is continuously applied to one of the liver lobes to generate observations (yellow points in Figure 9.1a). The amplitude of such force follows the sinusoidal function

1 2 × (1 -cos(2 × π × τ ))
where τ is a period. In this case, the control features defined in Section 5.3.1 correspond to the force applied (that is known). The Young's modulus is set to 5, 000P a in the reference simulation.

For the initialization of the ROUKF, we set µ 0 to 1, 000P a and σ 0 to 200P a. The state vector contains all the degrees of freedom of the mesh and the parameters to estimate (one parameter in our case). Hence, there are only 2 sigma points which allows a very fast assimilation process to take place as only two evaluations of the model need to be performed at each prediction phase. As depicted in Figure 9.1b, Towards Patient-Specific Networks 106 the value of the Young's modulus reaches rapidly a value close to the ground truth (at iteration 150, µ = 4948 and σ = 97). The value of the parameter converges to 4992 ± 15P a in 631 iterations (that is 500 seconds). Such assimilation process could take place before the surgery starts. Note that if the assimilation needs to be done in real-time, the simulations could be parallelized and simplified (we chose here a relatively high mesh resolution). Apart from material properties, the same approach can be employed to estimate the unknown attachments of the liver. We consider a scenario where an in vivo porcine liver is deformed with laparoscopic pincers. A video sequence of 7 s was acquired with a laparoscopic monocular camera inserted in the porcine abdomen inflated with gas. We assume there is a region Σ on which hidden boundary conditions are applied. Before the intervention, a CT scan was collected from the organ's geometry and a FE model was built following the pipeline described in Section 2.2. A tetrahedral mesh having 315 nodes was generated. The obtained model is fixed with elastic springs in a region Σ that is hidden to the laparoscopic camera view. There are 35 nodes in Σ meaning that 35 nodes were attached with springs to mimic the boundary conditions at these specific locations. The elasticity parameters of such springs are modelled as stochastic parameters (see Section 5.3.1). The elasticity values can range from 0 (no attachment) to high values (stiff attachment). A different elasticity is associated with each spring. Known surface displacements are prescribed on a small area of the visible surface to mimic the effect of the surgical tool based on the control features. The considered scenario falls within the category of displacement-zero traction problems, where the relation between surface and volume displacements is independent of the Young's modulus, for homogeneous materials [START_REF] Miller | On the prospect of patient-specific biomechanics without patient-specific properties of tissues[END_REF]. As a consequence, without lack of generality we set Y to the fixed value of 5, 000 P a, which is the average stiffness value for a liver. For the initialization step of the ROUKF, µ = 0 and σ = 0.01 for each parameter. Three assessment points are placed on the surface of the liver mesh in order to compute the prediction error between the observed data and the model. We ensure that the assessment points are different from the control and the observation features (see Figure 9.2a). In Figure 9.2b is shown the temporal evolution of the prediction error computed over the three assessment points. The error achieved with the stochastic simulation is compared to the error obtained with either fixed BCs (e.g. stiff attachment) or without BCs (e.g. spring stiffness set to zero). Results show that the stochastic simulation leads to smaller errors than the deterministic simulation.

Boundary conditions

The identified parameters can now be used to make patient-specific predictions by feeding them as input to the network (see Section 9.3). On the other hand, the obtained stochastic simulations can be used to generate new patient-specific data in order to retrain our network. To this end, in Section 9.4 we will explain how we can employ transfer learning methods to make the U-Net able to generalize to new unseen boundary conditions and elastic parameters.

Patient-specificity as input to the network

In the view of training a general network capable of predicting the displacement field of any liver geometry, with any type of BCs and material parameters, it is essential to teach the network to elucidate such variables. A solution might be to encode all these variables in the input to the network. That way, a universal network can be trained with a large variety of combinations of the sampled parameters. However, exhaustively sampling all the possible combinations would be extremely time consuming and probably useless to the network. Instead, a simple option consists of randomly sampling the parameter's space as done in [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration using a Data-Driven Biomechanical Model[END_REF], Pellicer-Valero et al., 2019]. Contrary to the approaches presented in chapters 4, 7 and 8 of this manuscript, the data generation and the training of a generic network would take more than a week as reported by [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration using a Data-Driven Biomechanical Model[END_REF]. Nevertheless, such expensive training is only performed once, and the trained model can be used for new liver geometries with new parameters as long as they are identified at prediction time. In this section we will see several strategies in order to encode patient variability in the input to the network.

Towards Patient-Specific Networks 108

Geometry and heterogeneity

As seen in Chapter 5, the first step towards patient-specific modeling of an organ is to account for its 3D geometry. The 3D surface of the organ can be segmented on a preoperative CT scan for instance. In order to encode such geometry surface in the input, several solutions have been proposed in the literature.

For instance, [START_REF] Pellicer-Valero | Real-time biomechanical modeling of the liver using Machine Learning models trained on Finite Element Method simulations[END_REF] provided a parameterization of the geometry by registering the liver surface to an average liver model. They performed a non-rigid registration using a modified Coherent Point Drift algorithm in order to express how much each node of the surface must be displaced (in each direction) in order to match the average liver surface. Then, Principal Component Analysis was used in order to reduce the size of the feature vectors to only 27 Principal Components. With this encoding procedure, they trained a feedforward neural network with tens of thousands of simulations on more than 100 liver geometries and obtained good accuracy on unseen liver geometries. Note that this work was only tested with synthetic data.

On the other hand, [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration using a Data-Driven Biomechanical Model[END_REF] proposed to express the preoperative volume geometry as a signed distance field. They generated a random 3D surface for each sample (e.g. 80, 000 different geometries) by extruding and remeshing an icosphere. They used a very similar network as ours so the obtained surface meshes must be encoded on a regular grid. To do that, at each grid node, the distance to the nearest surface point was computed. The points inside the organ volume take a minus sign thus leading to a signed distance field encoding the geometry. They trained the network with 460, 000 samples and tested it on real laparoscopic data. The obtained accuracy is not satisfactory yet.

The preliminary solution we tried consisted in encoding the geometry as a binary mask over the input grid (ones on the nodes of the sparse grid and zeros elsewhere). We trained our U-Mesh with ten different geometries generated from a beam to which some cells were removed. We ensured that each geometry reacted in a significantly different manner to a given external force. We generated 160, 000 training samples in total and trained the network for 100 epochs. Each sample ({F, M }; U v ) corresponds to an input tensor {F, M } of size n x × n y × n z × 4 and an output tensor U v of size n x × n y × n z × 3. F is the applied force encoded in the grid (3 values per node), M is the binary mask encoding the preoperative geometry (one value per node) and U v is the corresponding volumetric displacement. We obtained promising results but we did not extend these results to liver geometries yet.

The overall idea would be to train a general network with various liver geometries encoded using one of the above presented strategies and then, during surgery, make predictions with the patient's anatomy expressed in the input. Moreover, Towards Patient-Specific Networks 109 if we were to use a mask, we could encode material heterogeneity (assuming it is known) by setting the mask value to the normalized Young's moduli. Note that in [START_REF] Pellicer-Valero | Real-time biomechanical modeling of the liver using Machine Learning models trained on Finite Element Method simulations[END_REF], the usage of a binary mask was not recommended as in feedforward neural networks such mask should be flattened into a vector. This would mean that two voxels in the same location from two different livers represent the same feature which is not necessarily true.

Material parameters and boundary conditions

The other patient-specific parameters are the material properties and the BCs. Material properties such as Young's modulus and Poisson's ratio can be given as extra channels in the input. With the U-Mesh, this means giving an input tensor of size n x × n y × n z × 3 + i, i being the number of extra parameters. We trained the U-Mesh with variable Young's modulus lying in [0P a; 10, 000P a]. To do that, we sampled the domain and generated 16, 000 samples of random deformations for different values of E. We trained the network for 100 epochs and get small and constant errors when testing on different Young's modulus falling within the training range. An important requirement seems to be the normalization of such values over the sampled domain.

On another hand, fixed boundary conditions could also be given as a binary mask in the input, by setting to one the value of the nodes that are fixed and to zero the free nodes. However, we did not test this idea since having completely fixed nodes is a simplification of our FE simulations. Indeed, a better approximation of the boundaries consists of using springs to model organ's attachments and encoding those in the input seems non trivial. In the following section, we will see an alternative to deal with parameters that are difficult to explicitly encode in the input.

Transfer learning for patient-specific simulations

Let us assume that we have trained the network for predicting the deformation of a liver following the procedure described in Chapter 7. At this stage, the U-Mesh has learnt the chosen stress-strain relationship with average boundary conditions and material properties. In this context, the geometry is patient-specific as it is obtained from the patient CT scan. Now, using the intraoperatively estimated parameters (see Section 9.2) we can build a patient-specific FE model in order to generate a new small training data set. In this section we will see how to use this data set in order to update the network's weights through transfer learning.

As mentioned in Chapter 4, there exists a correlation between our method and model reduction techniques. There is an important body of work in this area, with a

Towards Patient-Specific Networks 110

well-established understanding of the process linking the fast (macro) model to the full (micro) model [De Angelo et al., 2019]. Such theory-driven approaches define how to generate reduced models with adapted parameters that characterize the full (micro-scale) model [START_REF] Boutin | Linear pantographic sheets: Asymptotic micro-macro models identification[END_REF]. Our deep learning approach does something similar by learning the key characteristics (deformation and parameters) of the full-scale model, but using a data-driven approach for this.

When applied in the context of surgery, both approaches share the same limitation.

The full model (micro-model) cannot always be correctly parameterized until the surgery has started, as some model parameters are not measurable in preoperative images. In this case, the use of transfer learning methods can offer a natural, datadriven solution for adapting the neural network to a particular patient. For methods based on reduced models, Bayesian approaches are probably a good alternative, as they can estimate material properties from a probability distribution and a priori knowledge of the parameter value.

As mentioned in the Section 5.3, boundary conditions have a significant impact on the accuracy of the predictions computed by biomechanical models. However, since they are hard to identify, we want to ensure the robustness of the U-Mesh to the variability of the BCs. We will show that a small amount of data is required to learn patient-specific BCs, when refining a network pre-trained with variable BCs from an appropriate distribution. This could help to significantly reduce the expensive cost of the offline data generation phase. Lastly, since real data can be sparse and noisy, we explore the behaviour of the U-Net when the input tensor C is highly sparse, and the effect of noise on the quality of the predictions.

Validation metrics

To assess the efficiency of our method, we perform a statistical analysis of the error over the testing data set {(C m , U vm )} M m=1 . Let U vm be the ground truth displacement tensor for sample m generated using the FEM described in Section 2.2 and h(C m ) the U-Mesh prediction. We define the mean Euclidean error e between U vm and h(C m ) for sample m as:

e(U vm , h(C m )) = 1 n n i=1 U v i m -h(C m ) i 2 (9.1)
where n is the number of nodes of the mesh. We compute the average e and standard deviation σ(e) of such norm over the testing data set. The mean Euclidean error represents the intuitive nodal distance, averaged over all the nodes of the mesh.

Towards Patient-Specific Networks 111 Table 9.1: Summary of the 7 strategies of interest. "BC4" stands for 4 adjacent cuboids in the middle of the hidden part of the beam. "BC8" stands for 8 adjacent cuboids fully covering the hidden part of the beam (see Figure 9.3).

Beam with hidden boundary conditions

In this section, we compare the accuracy of the U-Mesh either when trained from scratch with up to 16, 128 samples, or when pre-trained on various BCs and refined on the target BCs. We consider a deformable beam (size 4x1x1 m 3 , E = 300 Pa, ν = 0.4, 500 regular hexahedral elements) subject to fixed boundaries on a rectangular cuboid of its bottom part (see Figure 9.3). The beam follows the Saint-Venant-Kirchhoff behavior described in Section 2.1.3. We generate three different training data sets (N 1 = 16, 128; N 2 = 1, 209; N 3 = 100) and one testing data set (M = 4, 032), all drawn from the same distribution. We performed 10 trainings to compare 7 different strategies, summarised in Table 9.1. In strategies 1, 2, and 4, the U-Net is trained from scratch whereas in strategies 3, 5, 6 and 7, the U-Net is refined starting from a network pre-trained with 16,128 data with different boundary conditions (see Table 9.1).

In Table 9.2 are reported the validation metrics computed for each strategy on the same testing dataset (M = 4,032), as well as the index of the best iteration over 200,000 (with a saving step of 5,000). We see that strategy 3 performs better than strategy 2. More impressive yet are the strategies 5, 6 and 7 (especially 5, which, by refining, led to errors comparable to the one obtained with 12x more data without refining). Furthermore, they are substantially better than strategy 4 where no refinement was done. This is an example of a scenario where the U-Net cannot accurately learn a deformation model from scratch with very few (100) data, whereas it does learn an accurate model in a few thousands iterations using transfer learning in a few minutes (less than 15 minutes for scenario 3 and roughly 5 minutes for scenario 5 to reach best iteration).

The mild differences between the metrics obtained for strategies 5, 6 and 7 show that the data generated for pre-training must be reasonably distributed. Indeed, even though the network benefits from the diversity of BCs encountered in the pretraining stage, it is more efficient when these BCs are close enough to the target boundary conditions. Hence the need for a reasonable distribution.

So far we have seen that refining from an average model significantly reduces the quantity of data required to learn a deformation model. Results in Table 9.3 highlight the fact that it also speeds up the model convergence. Computing more metrics, we found that a good accuracy is reached approximately 20x faster when refining with 100 data than when starting from scratch with either 16,128 or 1,209 samples.

For completeness, we also investigated the case where the constitutive law changed between pre-training and refining stages. For pre-training, we modelled a beam with the linear Hooke's law, and for refining, we chose the Saint-Venant Kirchoff constitutive equation to model the deformations of the beam. In this scenario as well, we observed that transfer learning reduces the amount of data required to reach a given accuracy even when the base equations are complexified. We relaunched the training of 6 and 7 with a step of 50 iterations to ensure there was no significant overfitting before iteration 5,000, hence the values "1.5" and "1.15" (as a matter of fact, iterations 1,500 and 1,150 were actually slightly better than 5,000). Transfer learning situations are highlighted in red.

New boundary conditions and sparse data

As mentioned in previous paragraphs, real intraoperative data can be sparse and noisy. In this section, we show that the U-Net can still learn models when the training input tensors only contain a sparse view of the displacement U s imposed on the upper surface Γ D . In return, the accuracy is reduced and we show that transfer learning is barely relevant in such an adverse scenario.

We consider the same beam as described in Section 9.4.2, except that the beam is supported on both ends (fixed beam). In order to train the network, we built 2 data sets of sizes N 1 = 10, 080 and N 2 = 1, 008 (see Section 9.4.2). Here, every tensor C contains the values of an imposed surface displacement, on a randomly selected subdomain of the upper surface (in between 13 and 114 non-zero displacements in the testing data set, 67 in average) -see Figure 9.4. We trained the network either directly with N 1 or N 2 samples (strategies 1 and 2), or with N 2 data after a pre-training on a stiffer beam fixed at one end (Young's modulus of 500 Pa). The pre-training was done either with sparse data (strategy 4, same distribution as the refining data set), or dense data (strategy 3, full view of the imposed upper surface displacement U s ). In Table 9.4 are reported the validation metrics at best iteration. The average error with N 1 = 10, 080 is of only 3.03 mm and the maximal error is 57.1 mm -meaning less than 1.5 % of the length of the beam as maximal error. This shows that even though U-Net may learn much more accurate models with dense data, it still deals pretty well with sparse data when provided with a large enough training dataset. We should mention that we obtained very similar results by applying an additive white Gaussian noise of variance N = 10 -3 m on the testing dataset. With a variance N = 10 -2 m, the average mean Euclidean error e barely exceeds one centimeter. On another note, we see that there is no meaningful difference between the validation metrics of strategy 2 (1,008 data without refining) and strategy 4 (1,008 data with refining). Eventually, except when the data set is very small, we found that refining doesn't enhance accuracy in such a scenario. What is more, these results highlight the importance of pre-training the network with sparse data whenever the refining data is sparse. We further investigated the case where only very few data (N 3 = 100) are available, and found that it was not sufficient (with or without pre-training), although the refined model was more accurate. What remains valid is that in any scenario, the U-Mesh maintained a better accuracy with transfer learning in the first thousands of iterations. Reiterating these tests with sparse data without modifying the Young's modulus between pre-training and training stage corroborated these results.

Conclusion

In this chapter we have proposed several strategies to fulfill the real-time and precision requirements of patient-specific augmented reality. Based on a priori knowledge of the biomechanics of the organ, we select a constitutive model describing the relation between stresses and strains. Such relation is heavily affected by patient-specific properties such as boundary conditions and material characteristics. While obtaining these properties preoperatively may be troublesome, having information about them intraoperatively can be straightforward.

In our approach, the parameters of the preoperative finite-element or deep learning

Towards Patient-Specific Networks 115 models are updated based on intraoperative observations exploiting Bayesian filtering. The obtained parameters can be directly given as input to the network or can be used to retrain it with new data-driven simulations. In the first scenario, expensive data generation and training phases must take place in order to allow the network to generalize to any new liver. This brute force solution has led to interesting results [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration using a Data-Driven Biomechanical Model[END_REF], Pellicer-Valero et al., 2019] so it seems worth investigating. On the contrary, employing transfer learning to update a deep learning model is a more subtle solution, that would probably require further research efforts. When using transfer learning, not only the elasticity parameters and the boundary conditions can be changed simultaneously, but also the constitutive model itself. An interesting point of using deep neural networks is that the parameters do not need to be explicitly identified as they are encoded in the data (unless required as input as in Section 9.3). Hence, the network builds its own representation and through transfer learning, the weights of the network can be modified to match the targeted function. Moreover, we have reasons to believe that the U-Net learns local correlations in the displacement field rather than an overall model only. As a consequence, if the pre-trained model represents an average liver, transfer learning should not break the constitutive model learned previously. Note that the variability between livers can be high but it will always vary in a bounded range. For this reason, we believe that transfer learning is the key to an accurate and fast simulation of the deformations of a liver.

Note that when using Bayesian filtering, each parameter can only be modified individually. Indeed, in our pipeline using Kalman filters, in order to estimate the elasticity of the boundary conditions, the Young's modulus of the material needs to be fixed (and vice versa). A simultaneous estimation of both sets of parameters would be more complicated (yet possible), less precise (variance of the stochastic parameters will remain high) and would require very tedious fine-tuning of the filter. An alternative would be to use deep neural networks to estimate simultaneously the boundary conditions and the material parameters, but this might require a sequence of deformations as input.

Contributions of this chapter:

• Intraoperative identification of the organ's boundary conditions and stiffness for patient-specific finite-element simulations • Proof of concept: transfer learning to adapt an average network to patient specificity at the beginning of the surgery • Publication: Mendizabal, A., Tagliabue, E., Hoellinger, T., Brunet, J. Chapter 10

General conclusion

The objective of this thesis work was to come up with new strategies for computing real-time and accurate deformations of a soft structure. Without lack of generality, we considered the deformation of soft tissue in medical applications. In virtual laparoscopic training, surgeons practice on a completely virtual environment thanks to real-time simulations accounting for the interactions between surgical tools and the tissues. As an additional example, in the context of augmented surgery, in order to accurately locate the internal structures of an organ, a 3D model of the anatomy is overlaid on the surgical view by performing an elastic registration between preoperative and intraoperative images. Also in image-guided biopsy or robotized procedures, having a better control on the interactions of the robot with soft tissues can improve the outcome of the procedure and the safety to the patient.

Throughout this manuscript, we showed how combining finite-element simulations with machine learning techniques can help in robotized interventions and surgical guidance. For instance, in the context of robotic assisted interventions, where the knowledge of the force applied by the robot is a key feature for the safety of the patient, we proposed an image-based method for estimating such force without the need of any force sensor. We used a simple fully connected neural network for force classification, where the input to the network is an intraoperative image of the deformed tissue. The network was trained with synthetic data generated with a patient-specific finite-element model of the considered organ.

The corpus of this work relies on the U-Mesh framework: a method fulfilling the real-time and precision requirements of patient-specific augmented reality. Based on a U-Net architecture, it can learn the relation between a constraint and a deformation for various geometries and make predictions with high accuracy in a very short amount of time regardless of the size of the problem. To make predictions on irregular geometries, we proposed an Immersed-boundary method to encode the displacement fields in regular grids handled by the network. In particular, we showcased the potentiality of our method in open liver augmented surgery in humans, with very little equipment (one single RGB-D camera). Based on a preoperative 3D model of the organ, we used the U-Mesh as a non-rigid registration method in order to fit the intraoperative deformations. Our deep learning model is driven by a surface point cloud acquired with one single RGB-D camera, thus allowing for a markerless and radiation-free setup for building up an augmented reality. Our approach could be integrated in any operating room and has great potential in complex liver surgeries where small tumors must be resected. To demonstrate the generality of our pipeline, we also applied the U-Mesh in the context of US-guided breast biopsy. The objective was to track the internal tumors in real-time accounting for the large deformations induced by the US probe. To train our network, we simulated probe-induced compression using the Immersed-boundary method. The reached target registration errors (4 mm in average) was comparable to that obtained with finite-element numerical simulations not capable of real-time.

In all the presented scenarios, the U-Net is trained with synthetic data based on finite-element simulations. While some aspects of the simulation can be set preoperatively (such as the organ's geometry and average location of ligaments for instance), some others (such as material elasticity and boundary conditions) can only be identified intraoperatively and may require further effort depending on the demanded accuracy. To this end we proposed to estimate such patient-specific parameters based on intraoperative observations through Bayesian filtering. The obtained patient-specific simulations can be used to update our deep learning models employing transfer learning. In a similar direction, [Chakraborty, 2020] proposes to train a low-fidelity physics informed model, and refine it with few high-fidelity data exploiting transfer learning.

Despite the promising results of our method, there are some limitations worth mentioning. The accuracy of the U-Mesh is only guaranteed as far as the inputs are included in the range of the training data sets. In the same manner, we are restricted to the geometry used to train the network (unless encoded in the input). This is due to the fact that neural networks are not good at extrapolating, hence the importance of a good sampling of the input domain when generating the data sets.

Moreover, even if the choice of our network was motivated by its similarities with model order reduction techniques, we cannot claim that it was the best and only choice. Indeed, we have tried a simple fully-connected network and it seems to produce similar results in terms of accuracy while being about three times faster to train. Nevertheless, we observed that the U-Net has better extrapolation capacity (thus better abstraction of the problem) than the fully-connected networks. This means that even if it is not possible to make very accurate predictions when applying a force somewhere out of the sampled input domain of the training data set, the U-Mesh performs "less bad" than other networks.

General conclusion 119

Another limitation of the method is the expensive offline phase. The data generation can be extremely time consuming in particular when considering large meshes or when more complex input sequences are needed. It goes without saying that the larger the data sets, the longer the training. Hence it is important to build smart data generation strategies to cover all the force ranges without being exhaustive (to reduce data generation and training times). A possible alternative would be to include in the dataset only those deformation modes which represent significantly different deformations. Moreover, we performed a preliminary study on the sensitivity of the method to the Young's modulus variations in the training data set. We noticed that U-Mesh is more accurate when dealing with stiffer objects. This is an important result to keep in mind for an eventual smart data generation. Indeed, in a low deformation regime fewer data are needed to reach good accuracy. In a similar direction, we noticed that when using transfer learning, pretraining the network on a slightly softer object than the targeted average, accelerates the convergence.

Future research directions

In this thesis work, we have shown the potentiality of using physics-based deep neural networks in computer-assisted interventions. We have proved that the accuracy of the predictions satisfies the precision requirements of medical applications such as targeting small tumors. The large field of possibilities and the trendiness of our topic open space for new exciting research directions.

In particular, if we aim at integrating the U-Mesh in any operating room for augmented surgeries, it seems worth investigating how to generalize our network to any kind of geometry using more advanced strategies than those proposed in Section 9.3. Indeed, encoding the geometry seems to be a delicate yet important point, in particular when using regular grids. It is non trivial to express a discrete surface in a grid without introducing approximation errors. In the view of an optimal data generation process, we would like to explore statistical atlases to sample the geometry space in a meaningful manner as well as the organ attachments as done in [START_REF] Plantefève | [END_REF].

Even if the obtained target registration errors seem compatible with oncological surgery needs, we could further improve the finite-element modeling used to generate our training data set. Since the data generation takes place in an offline phase, we could enrich our models as much as needed. For example we could include incompressibility, account for heterogeneity, use quadratic hexahedral elements or simply refine the mesh. If, besides the displacement field, we also want to estimate the stress distribution within the organ, we could use more complex hyperelastic laws or expand our model to account for viscoelasticity. Handling models such as Ogden or Mooney-Rivlin would require no change to our method. Theoretically the architecture of the neural network is independent of the mechanical model but in practice, if the complexity of the model increases it might be relevant to deepen the network or enlarge the training data sets. On the other hand if we were to handle viscoelasticity, we would need to include a time term and the current state of the system in the network input as done by [START_REF] Meister | Towards Fast Biomechanical Modeling of Soft Tissue Using Neural Networks[END_REF].

In order to improve the robustness and accuracy of the method, we also plan to revise the surface point matching between the preoperative and the intraoperative surfaces. For now, we are using the iterative closest point algorithm but it might be relevant to also let the network guess the surface correspondences as done in [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration using a Data-Driven Biomechanical Model[END_REF]. However, an accurate initial rigid alignment of the surfaces seems crucial for such method to work, as the network produces a result in one iteration. Conversely, we could use the U-Mesh in an iterative process where at each iteration, a displacement can be imposed on a grid node based on the position of its closest point in the RGB-D point cloud. In a similar direction, we plan to integrate multiple RGB-D cameras in our set up, in order to improve the quality of the point cloud describing the surface deformation (larger field of view and less noise), or to use the second camera point cloud for validation purposes.

Another scenario worth investigating is the case of contacts between anatomical structures. Interaction between objects can be seen as external forces applied to their surfaces (the alternative option being to solve interactions through position constraints). Assuming we have two objects embedded in two U-Mesh grids, we can compute their motion until a contact is detected and then apply a simple penaltybased contact response. This contact response is a force applied on the surface of each object to cancel out their interpenetration. Using this force, we could then compute the deformation of each object following our method.

The work presented throughout this manuscript is an interface between precise biomechanical modeling (not capable of real time) and clinical applications requiring both high accuracy and very high speed. Even if there is still room for improvement, our pioneer method has proved successful in scenarios such as augmented surgery and image-guided biopsy. Pour modéliser la déformation des tissus mous, la méthode des éléments finis (EF) est la technique préférée de part sa capacité à résoudre des problèmes complexes, sa fiabilité et sa robustesse. Pour répondre aux critères de précision de la chirurgie guidée par l'image (erreur de recalage de l'ordre de 5 mm [START_REF] Ruiter | Model-based registration of X-ray mammograms and MR images of the female breast[END_REF]), le model doit prendre en compte les spécificités de chaque patient notamment d'un point de vue anatomique mais aussi biomécanique et physiologique. Par exemple, les fixations de l'organe ainsi que ses propriétés biomécaniques jouent un rôle très important dans le calcul de sa déformation. Par ailleurs, pour qu'il soit pertinant, le modèle doit faire des prédictions en temps-réel, ce qui est difficile de garantir en pratique surtout lorsque des déformations non-linéaires sont impliquées.

De nombreux travaux ont visé à repousser les limites de temps de calcul des méthodes EF. Très souvent, il s'agit d'un compromis entre temps de calcul et précision. Lorsqu'on cherche à modéliser des petits déplacements, les hypothèses de linéarité de la loi de comportement accélèrent considérablement les calculs. Par contre, lorsque les déformations sont plus importantes, ces simplifications ne sont plus représentatives de la réalité. Une solution serait alors d'utiliser les EF co-rotationnels qui prennent en compte les non-linéarités géomètriques en tempsréel [START_REF] Petit | Environment-aware non-rigid registration in surgery using physics-based simulation[END_REF], Haouchine et al., 2013a]. Cependant, cette approche ne s'étend pas à des lois de comportement plus sophistiquées comme les lois hyperlastiques pour lesquelles, d'autres alternatives ont été proposées. Par exemple, dans [START_REF] Marchessau | [END_REF], les auteurs utilisent la méthode dite MJED (Multiplicative Jacobian Energy Decomposition) afin de produire en temps-réel les deformations d'un foie hyperelastique, poreux et visqueux. Aussi, les auteurs de [START_REF] Miller | Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF] ont introduit les TLED (Total Lagrangian Explicit Dynamics) garantissant le temps-réel lorsqu'il sont associés à des schémas d'intégration explicite et des solveurs GPU [START_REF] Joldes | Real-time nonlinear finite element computations on GPU-Application to neurosurgical simulation[END_REF].

Plus récemment, l'emmergence des méthodes d'apprentissage automatique (machine learning (ML) en anglais) ont permis d'approximer une fonction quelconque, du moment qu'il y a suffisament de données d'entrainement, et ce, sans aucune connaissance a priori du problème. Aussi, ces méthodes font des prédictions en des temps très courts ce qui les rend très interessantes pour nos domaines d'application. En effet, le comportement biomécanique d'un tissus mous peut être implicitement encodé dans un model ML, pour prédire en temps-réel les déformations d'un organe en réponse à une traction de surface [START_REF] Morooka | Real-time nonlinear FEM with neural network for simulating soft organ model deformation[END_REF], Tonutti et al., 2017, Rechowicz et al., 2013, Mendizabal et al., 2019] ou bien à un déplacement imposé [START_REF] Pfeiffer | Learning soft tissue behavior of organs for surgical navigation with convolutional neural networks[END_REF], Brunet et al., 2019, Lorente et al., 2017, Mendizabal et al., 2019]. Cependant, la précision de la prédiction d'un réseau est fortement liée à la qualité et à la quantité des données d'entrainement. Dans l'idéal, il faudrait une infinité de données réelles et non bruitées du phénomène à approximer, ce qui est en pratique très difficile à assurer. C'est pour cela que les simulations EF peuvent être exploitées pour générer des bases d'entrainement représentatives de la réalité.

Parmis les méthodes de ML, l'utilisation des réseaux de neurones (RN) a considérablement augmenté. Par exemple, Tonutti et al. utilisent un RN pour prédire les déplacements de tumeurs cérébrales en réponse à des forces. De même dans [START_REF] Rechowicz | Development and validation methodology of the Nuss procedure surgical planner[END_REF] Le manuscript de thèse est divisé en 7 chapitres principaux. Les deux premiers pausent les bases théoriques des sujets que l'on considère, à savoir la formulation du problème bioméchanique et sa résolution par la méthode des EF, ainsi que les fondements théoriques des réseaux de neurones qui seront utilisés. Le chapitre 4 constitue un projet de classification de force lors d'interventions robotizées. Le chapitre 5 présente la plateforme U-Mesh qui est le coeur de cette thèse. Les chapitres 6 à 8 comportent les applications de U-Mesh dans des contextes chirurgicaux.

Estimation de force avec des réseaux de neurones

Les maladies ophtalmologiques nécessitant des injections intravitréales sont en hausse avec plus de 4 million d'injections dans le monde en 2014 [START_REF] Ullrich | Assistive Device for Efficient Intravitreal Injections[END_REF]. C'est pourquoi plusieurs travaux de recherche proposent de robotiser la procédure [START_REF] Ullrich | Assistive Device for Efficient Intravitreal Injections[END_REF] afin de libérer les cliniciens et permettre aux hopitaux d'assurer une demande toujours en hausse. Pour designer un tel système, il est nécessaire de développer un robot de précision pour assurer l'efficacité de l'injection ainsi que la sureté du patient. En effet, l'injection doit être faite dans une région très petite se situant entre la cornée et la sclère, sinon l'oeil risque d'être endommagé. Aussi il est important de connaître la force appliquée par l'aiguille sur la sclère afin de garantir la sécurité du patient [START_REF] Jagtap | Applied force during vitreoretinal microsurgery with handheld instruments[END_REF]. Les senseurs de forces étant couteux et compliqués à utiliser en pratique [START_REF] Haouchine | Vision-based Force Feedback Estimation for Robot-assisted Surgery using Instrument-constrained Biomechanical 3D Maps[END_REF], Haidegger et al., 2017], il est préférable d'estimer les forces en se basant sur des images intraopératoires [START_REF] Mura | Vision-based haptic feedback for capsule endoscopy navigation: a proof of concept[END_REF], Haouchine et al., 2018].

Dans ce chapitre on propose d'utiliser un réseau de neurones pleinement connectées (c'est à dire sans couche de convolution) pour classifier la force exercée par l'aiguille en se basant uniquement sur des images OCT de la sclère déformée. En effet pendant l'intervention, il est possible d'acquérir des images OCT 2D de la sclère (voir les images en noir et blanc de la figure 10.2) en capturant la déformation induite par l'aiguille. L'objectif étant d'identifier le moment où le robot excerce une force trop importante, on se propose de classifier les forces en trois rangs :

• rang 0 : peu ou pas de déformation,

• rang 1 : la sclère se déforme mais il n'y a pas de danger pour le moment,

• rang 2 : la sclère se déforme considérablement et on est sur le point de percer la sclère. Il est préférable de rétrograder l'aiguille.

Comme citer précédemment, la qualité de la prédiction d'un réseau de neurones dépend fortement des données d'entrainement. On propose donc de construire un modèle bioméchanique de la sclère afin de générer autant de données que nécessaire. Pour cela on modélise la sclère comme une demie-sphère avec une épaisseur non nulle, fixée aux extrémités (conditions aux limites de Dirichlet) et on applique une pression intraoculaire (IOP) (voire figure 10.1). Pour simuler la force exercée par l'aiguille, on applique une force locale en une zone petite autour de l'extrémité de l'aiguille.

On résout les équations obtenues avec la méthodes des EF. Pour cela, on discrétize la demie-sphère avec des élèments hexahédriques (H8) et on cherche l'équilibre des forces (internes et externes) avec un algorithme de Newton-Raphson. La paramétrization du modèle méchanique ainsi que la géométrie du modèle sont très importantes pour que les déformations soient représentatives de la réalité. Pour commencer, l'épaisseur de la sclère joue un rôle non négligeable : pour une force donnée, une sclère épaisse se déformera moins qu'une sclère fine. On va donc faire varier l'épaisseur de la sclère de 400 à 800 micromètres [START_REF] Olsen | Porcine sclera: thickness and surface area[END_REF]. D'un autre côté, la pression intraoculaire influence aussi la déformation. Elle peut être mesurée avec un tonomètre disponible en consultation ophtalmologique. Dans nos expériences, on utilise des yeux de cochons postmortem qui ont des pressions intraoculaires négligeables (2 mmHg). On peut donc la considérer comme constante dans nos simulations. Un autre paramètre important est la raideur de la sclère (le module d'Young) qui est variable selon les espèces de cochons et selon les conditions expérimentales. On utilise donc des filtres de Kalman pour estimer sa valeur afin de paramétrer nos simulations.

Notre base de données d'apprentissage consiste en des couples (F, I) avec F , le rang de force appliquée et I, l'image de la déformation résultante. L'image est obtenue en prenant une section 2D du modèle EF 3D. La figure 10.2 montre des exemples d'images OCT réelles vs images simulées pour les trois rangs de force. Une fois les données générées, on peut entrainer le RN. Il s'agit d'un réseau de classification avec deux couches cachées de 600 neurones chacune et une couche de sortie avec 3 neurones (un neurone par classe ou rang de force). L'entrée du réseau correspond à l'image OCT (réelle ou simulée) et la sortie (ou label) correspond au rang de force appliquée.

Figure 10.2: Images réelles OCT de la sclère ainsi que les simulations correspondantes pour les trois rangs de force.

Nous avons testé notre méthode sur des yeux de cochons avec 93% de précision dans les prédictions. Néanmoins, pour inclure cette méthode dans une procédure clinique, il faudrait impérativement inclure une valeur de pression intra-oculaire spécifique au patient, ce qui rendrait le modèle bioméchanique plus compliqué. De plus, la mesure de cette pression est souvent biaisée par la raideur de l'oeil. Il faudrait donc identifier ces deux paramètres conjointement lors de l'assimilation avec les filtres de Kalman.

Plateforme U-Mesh: estimation du champ de déplacement avec des réseaux de neurones

Il existe de nombreuses applications où il est nécessaire de calculer les déformations de structures non-linéaires en temps-réel. Comme vu en introduction, la méthode des EF est la plus utilisée pour résoudre ce type de problèmes. Cependant, dès lors que l'on considère des matériaux non-linéaires, il s'avère difficile d'assurer la contrainte de temps réel. Pour accélérer les temps de calculs, plusieurs méthodes ont étaient proposées comme la décomposition de domaine [START_REF] Haferssas | An Additive Schwarz Method Type Theory for Lions's Algorithm and a Symmetrized Optimized Restricted Additive Schwarz Method[END_REF], le calcul parallèle sur GPU [START_REF] Allard | SOFA -an open source framework for medical simulation[END_REF], Johnsen et al., 2015] et la réduction de modèle [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF], Ryckelynck, 2005, Goury and Duriez, 2018]. Plus récemment, les méthodes d'apprentissage automatique ont été proposée pour modéliser la méchanique des objets [START_REF] Lorente | A framework for modelling the biomechanical behaviour of the human liver during breathing in real time using machine learning[END_REF], Tonutti et al., 2017, Morooka et al., 2008] pour de petites déformations et un ensemble de forces restreint. Dans ce chapitre on présente la plateforme U-Mesh qui produit la déformation de géométries simples pour de grandes déformations avec un nombre variable de forces en entrée. On a choisie une architecture de type U-Net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], pour ses similarités avec la réduction de modèle (voire figure 10.3). Il s'agit d'un réseau de type encodeur-décodeur avec des filtres de convolutions pour extraire les traits d'intérêts des entrées du réseau.

Pour entrainer un tel réseau, on génère une base d'apprentissage avec les EF, constituée de couples (F, U ) où F est l'entrée du réseau et correspond à la force appliquée, et U est le déplacement volumique résultant (sortie du réseau). Pour cela, on construit un modèle EF de l'objet traité et on applique des forces aléatoires à la surface de l'objet. À chaque simulation, on enregistre le couple (F, U ). Un détail important pour que la méthode fonctionne, réside en la nécessité d'encoder les forces et les déplacements dans des grilles régulières. En effet, le réseau U-Net nécessite une structure de grille pour que les filtres de convolutions puissent extraire des informations spatiales cohérentes. On utilise donc des grilles héxahédriques pour nos simulations EF. Ainsi, les noeuds de la grille contiennent le déplacement ou la force des noeuds du maillage EF.

On a testé la méthode sur des géométries de références en méchanique, satisfaisant la structure de grille régulière. Tout d'abord nous avons modéliser le comportement d'une poutre à section carrée, fixée en une de ses extrémités. Une ou plusieurs forces sont appliquées aux noeuds de la face supérieure et données au réseau en entrée. Nous avons aussi modéliser le comportement d'un L. Sur la figure 10.4, apparaissent en vert les solutions EF et en bleu les prédictions U-Mesh. La différence est quasi imperceptible à l'oeil nu. En effet, l'échantillon avec l'erreur maximale (image de gauche) correspond à une norme l2 relative de 5.1% à l'extrémité libre de la poutre.

Le résultat le plus marquant est que les prédictions sont faites en 3 ms seulement.

U-Mesh pour la chirurgie augmentée du foie

Comme cité précédemment, pour construire une réalité augmentée d'un organe, il est nécessaire de faire un recalage non rigide entre les images préopératoires (scan CT et modèle EF) et les images intraopératoires (nuage de points et vue de la surface de l'organe). Dans ce chapitre nous allons voir comment faire ce recallage élastique en utilisant U-Mesh, pour prédire la déformation volumique de l'organe à partir d'un nuage de points surfacique clairsemé.

Comme vu dans le chapitre précédent, pour que U-Mesh puisse faire des prédictions, il est impératif d'encoder le champ de déplacement dans une grille régulière. Pour ce faire, plusieurs options s'offrent à nous. La première consiste à mailler la géométrie de l'objet avec n'importequel type d'éléments (tétrahèdres par exemple) pour y réaliser les calculs EF. Ensuite, le maillage peut être plongé dans une grille régulière et

U-Mesh pour la biopsie du sein

Lorsqu'il y a des lésions suspectes dans une mammographie, il est préférable de réaliser une biopsie afin de vérifier leur malignité. Pour guider l'aiguille vers la lésion, la biopsie est souvent guidée par des images ultrasons (US) [START_REF] O'flynn | Image-guided breast biopsy: state-of-the-art[END_REF]. Cependant, il n'est pas toujours facile de suivre les lésions sur ces images. En effet, elles ne sont pas toujours visibles suite au contraste réduit des US, et de plus, la navigation 3D est faite en se basant sur des informations 2D. Puisque nous disposons du scan CT préopératoire avec la localisation exacte des lésions, trouver une méthode qui mettrait à jour la positions des lésions pendant l'intervention serait très bénéfique. Comme pour la chirurgie augmentée du foie, on peut faire un recalage élastique entre le modèle préopératoire et l'image US en utilisant la méthode des EF. Là aussi, le sein est un organe hautement déformable et donc des équations non-linéaires entrent en jeu, ce qui rend les calculs incompatibles avec le temps-réel. Dans [START_REF] Han | A nonlinear biomechanical model based registration method for aligning prone and supine MR breast images[END_REF], les auteurs proposent d'utiliser une formulation TLED basée sur GPU permettant de calculer les déformations du sein en 30 ms, ce qui ne satisfait pas encore la contrainte de temps-réel.

Dans ce chapitre on propose d'utiliser la plateforme U-Mesh pour prédire la déformation du sein en réponse aux forces de compression induite par la sonde US. Dans ce cas, l'entrée du réseau comporte les déplacements des noeuds en contact avec la sonde US. La méthode est validée sur des données issues d'un fantome du sein. Une fois que les données sont générées et que le réseau est entrainé (∼ 5 heures), on prédit des déplacements sur des données synthétiques. On obtient des erreurs nodales moyennes de 0.05 mm ce qui montre bien encore une fois que le réseau a parfaitement appris le modèle avec lequel il a été entrainé. Dans un second temps, nous avons validé notre approche sur des données d'un fantome en sillicone contenant 10 lésions. On impose des déformations avec une sonde US qui est suivie dans l'espace. On connaît donc sa position et son orientation à chaque instant. Ainsi, on peut extraire l'emplacement des noeuds à la surface du sein en contact avec la sonde.

Pour chacune des 10 lésions, on place la sonde de telle sorte que la lésion considérée soit visible sur l'image US. Dans un premier temps, la sonde touche à peine le sein (donc aucune déformation n'est induite). Cette première étape permet de calibrer le système afin de référencer le modèle et la sonde dans un même repère. Ensuite on impose 4 déformations d'amplitudes croissantes. Ainsi on obtient 40 acquisitions avec des déformations variables. Lorsqu'on compare la prédiction U-Mesh avec la vérité terrain (obtenue à partir de l'image US), on obtient un TRE de 4.2 mm en moyenne. Encore une fois, on constate que U-Mesh fait aussi bien que la méthode EF utilisée pour l'entrainer.

En résumé nous avons construit une plateforme faisant l'interface entre des simulations EF précises mais couteuses en temps de calcul, et des applications nécessitant précision et rapidité simultanément.

Conclusion

Dans cette thèse nous avons démontré comment l'utilisation d'algorithmes d'apprentissage automatique en combinaison avec des méthodes numériques classiques, pouvait améliorer les interventions assistées par ordinateur. En particulier, nous avons présenté la méthode U-Mesh basée sur des réseaux de neurones pronfonds, qui modélise le comportement biomécanique d'un organe, tout en respectant les spécificités de chaque patient et en satisfaisant la contrainte de temps-réel inhérente à la chirurgie guidée par l'image.

Dans une phase préopératoire, on construit un modèle éléments-finis de l'organe considéré avec des paramètres moyens, afin de générer une base d'apprentissage pour entrainer le réseau. Ensuite, pendant la chirurgie, les propriétés mécaniques des organes sont identifiées avec de l'inférence Bayesienne et de nouvelles données patient-spécifiques peuvent être générées. Les poids des réseaux préalablement entrainés, peuvent être mis à jour par transfert de connaissance en utilisant les nouvelles données simulées.

•

  ∇f : gradient of f• ∂f ∂x : partial derivative of f with respect to x• I: identity matrix in R 3×3 • tr(A): trace of A • div(f ): divergence of f • u: first derivative of u • ü: second derivative of u • p ∼ N (µ,σ): variable p follows a normal distribution of mean µ and standard deviation σ Chapter 2 • f ext : external forces • f int : internal forces • m: body mass • a: acceleration • ρ: mass density • v: velocity • u: displacement field • Ω: computational domain • Γ: boundary of Ω • Γ D : subpart of Γ with Dirichlet boundary conditions vi • Γ N : subpart of Γ with Neumann boundary conditions • p: initial position of a particle in Ω • p: deformed position of p • φ: transformation undergone by the solid • F: deformation gradient in R 3×3

•

  Abbreviations• ANN: artificial neural network • AR: augmented reality • BC: boundary condition • CAI: computer-assisted interventions • CAS: computer-aided surgery

•

  NN: neural network • OCT: optical coherence tomography • PBD: position based dynamics • PCA: principal-component analysis • PDF: probability density function x • PGD: proper generalized decomposition • POD: proper orthogonal decomposition • ReLU: rectified linear unit activation function • ROUKF: reduced-order unscented Kalman filter • SGD: stochastic gradient descent • T4: 4-node tetrahedral elements • TRE: target registration error • US: ultrasound • VR: virtual reality xi Chapter 1

Figure 2

 2 Figure 2.1: Infinitesimal volume tetrahedron. The stress vectors t 1 , t 2 , t 3 and t n are linked to the normals -x 1 , -x 2 , -x 3 and n.

Figure 2 . 2 :

 22 Figure 2.2: Linearization of the relation between stress and strain. Source: www. simsolid.com

Figure 3

 3 Figure 3.1: Multi-layer perceptron with five inputs, five outputs and two hidden layers of neurons. Source: https://towardsdatascience.com

Figure

  Figure 3.2:

Figure 4 . 1 :Figure 4

 414 Figure 4.1: General network architecture for an object with a resolution of x × y × z nodes, c channels in the first layer and k steps. Note that x 0 × y 0 × z 0 is the padded grid.

Figure 4

 4 Figure 4.3: Front and side views of the maximal error sample obtained with a network having 3 steps and 128 channels. U-Mesh output is in blue and the reference is in green. The red arrow represents the force applied. Note that the difference between the two meshes is very small. The relative l2 norm at the tip of the beam is of 5.1% and the deformation amplitude is 0.45 m

  Figure 4.7). There are 336 nodes on the upper face and Λ is set to 100. Overall the data set has 33, 600 samples (N = 26, 880 samples for training and M = 6, 720 samples for testing). The U-Net is trained with the previouslyThe U-Mesh Framework 47 3: Error measures for a beam having 3267 H8 elements and one simultaneous force application.

Figure 4

 4 Figure 4.6: Sensitivity of e to the deformation amplitude for a 3, 267 H8-elements beam. The point cloud represents the values of e for randomly selected samples of the testing data set. The regression line of equation y = 0.0083 × x shows the low sensitivity of the U-Mesh to the deformation range. The average computation time is 3 ms.

Figure 4

 4 Figure4.7: Beam mesh with 3267 H8 elements deformed with U-Mesh, front and side views. U-Mesh output is in blue, the reference is in green and the rest shape is in grey. The computation time is equal to 2.9 ms for this sample. The relative l2 norm at the tip of the beam is 1.6% and the deformation amplitude is 0.4 m

Figure 4

 4 Figure 4.8: Samples of deformed L-shapes. In green is shown the reference solution, in blue the output of U-Mesh and in grey the rest shape.

Figure 4

 4 Figure 4.9: Sensitivity of e to the deformation amplitude for the L-shape. The point cloud represents the e all the samples of the testing data set. The regression line of equation y = 0.002046×x shows the low sensitivity of the U-Mesh to the deformation range. Maximal error sample highlighted in red and shown in Figure 4.10.

Figure 4 .

 4 Figure 4.10: Maximal error sample (e = 0.0327 m for a nodal deformation reaching 3.08 m). Reference solution is in green and U-Mesh prediction is in blue.

Figure 5

 5 Figure 5.1: Sequence of 2D CT images of the liver and resulting surface mesh.

Algorithm 2 :

 2 Main steps of unscented Kalman filter 1: Initialize data: 2: set x 1 -model positions and unknown parameters 3: set T = T (x 1 ) -finite element model 4: set I, P 1 , Q, W -initial filter parameters 5: for each simulation step i do 6:

  (a) Simplified model of the sclera (b) Deformed hexahedral mesh. Note that the colors of the hexahedra are just for visibility

Figure 6

 6 Figure 6.2: Artificial neural network used for image classification. The input is an image of size 140 × 40 pixels resulting in an array of size 5600.

Figure 6 . 3 :

 63 Figure 6.3: Continuous acquisition of the robot force sensor and corresponding punctual OCT images for one eye. The vertical lines give the imaging times.

Figure 6

 6 Figure 6.4: (a) Observation points in the real OCT image. (b) Predicted observations (in red) and real observations (in yellow) for the first estimation of E. (c) Predicted observations (in red) and real observations (in yellow) for the final estimation of E.

Figure 6

 6 Figure 6.5: Variation of µ and σ for the Young's modulus estimation using the ROUKF for 24 OCT images. The value of the parameter converges to 0.27 M P a.

Figure 6

 6 Figure 6.6: Real OCT images of the eye used for the data assimilation and matching simulated images with E = 0.25 M P a for the three force ranges.

Figure 6

 6 Figure 6.7: Real OCT images of the sclera and their matching simulated images with E = 0.25 M P a for the three different force ranges.

  (a) Loss and Accuracy curves for training and validation sets (b) Fragment of the training data set generated by our numerical simulation (a) Unprocessed OCT image (b) Filtered OCT image

Figure 6 .

 6 Figure 6.10: Filtered force values corresponding to the 45 OCT acquisitions. The three misclassified samples are overestimated. Note that we only kept the samples under 0.06N which is the point of rupture of the thinnest scleras. To correctly account for the elements above that point of rupture, the puncturing process should be included in our model which would complicate the simulation without bringing any gain.

Figure 7

 7 Figure 7.1: Augmented reality pipeline: preoperative internal structures are mapped in real-time onto the live image of the organ using a FEM model.

Figure 7 . 2 :

 72 Figure 7.2: Example of sparse grid discretization generated from a preoperative CT. Combined with an Immersed boundary method, it allows the use of regular hexahedral meshes to simulate the deformation of the organ, and can be used with our CNN.

  Figure 7.3: Various liver samples from the testing data set and corresponding relative errors computed on the tip of the deformed lobe for the tetrahedral topology. The rest shape of the liver is shown in grey.

Figure 7

 7 Figure 7.4: (a) Maximal average nodal error. (b) Regular grid mapped onto a FE tetrahedral mesh.

( a )

 a Figure 7.6: Various liver samples from the testing data set and corresponding relative errors computed on the tip of the deformed lobe. The rest shape of the liver is shown in grey.

Figure 7 . 7 :

 77 Figure 7.7: Sample with maximal maximal nodal error (0.0153 m) for the hexahedral mesh. The reference solution is shown in red and the U-Mesh predictions is in green. The rest shape is shown in grey.

Figure 7 . 8 :

 78 Figure 7.8: Sensitivity of e to the deformation amplitude for the liver. The point cloud represents the error e of all the samples of the testing data set. The regression line of equation y = 0.0021 × x shows the low sensitivity of the U-Mesh to the deformation range.

Figure 7

 7 Figure 7.10: (a)Virtual point cloud on the visible surface of the liver to generate random displacements (top view). (b)The U-Mesh prediction appears in red and the intraoperative CT scan is in gray (side view).

Figure 7 .

 7 Figure 7.11: Augmented reality frames obtained with U-Mesh predictions. The vascular tree appears in blue, the tumors are in gray and the parenchyma is in brown.

Figure 7

 7 Figure 7.12: Overview of the possible approaches to map a point-to-point displacement in grid coordinates. Blue curve: preoperative surface model. Red curve: intraoperative surface (i.e. the RGB-D point cloud). Purple curve: surface solution computed with the FEM. p 2 is the position of a point in the RGB-D point cloud and p 1 is the position of its corresponding point in the undeformed surface mesh (correspondence found with the ICP algorithm). (a) Set p 2 to the grid node that is the closest to p 1 . (b) Set p 2 to the 8 nodes of the cell containing p 1 . (c) Store the positions of the deformed grid nodes obtained with FEM for each point of the RGB-D point cloud, and create a virtual point cloud {n i } i=1...n in grid coordinates. (d) For each node in the grid, compute its distance with the closest point in the RGB-D point cloud.

Figure 8 . 1 :

 81 Figure 8.1: Breast surface mesh obtained from a preoperative CT scan Immersed in a hexahedral grid for FEM computations.

Figure 8

 8 Figure 8.2: (a) Experimental setup. From left to right: monitor showing real-time US images; CIRS breast phantom during FUS acquisition; optical tracking system that allows to map the real positions of the CIRS breast phantom and the US probe to the preoperative geometry model. (b) External surface and inner lesions of the CIRS breast phantom.

Figure 8

 8 Figure 8.3: (a) Sample with maximal deformation (79.09 mm). (b) Sample with maximal mean norm error (0.266 mm). The green mesh is the U-Net prediction and the red mesh is the FEM solution. The initial rest shape is shown in grey. (c) U-Net prediction on phantom data.

TowardsFigure 9

 9 Figure 9.1: (a) Liver simulation mesh made of 11, 000 tetrahedral elements. The red points highlight the fixed points (Dirichlet boundary conditions), the yellow points correspond to the observed features, and the green arrows illustrate the direction of the applied forces. (b) Variation of µ and σ for Young's modulus estimation using the ROUKF. The value of the parameter converges to 4992 ± 15P a in 500 seconds.

Figure 9

 9 Figure 9.2: (a) The first frame of the video sequence with features. (b) Temporal evolution of the assessment error computed for each assessment point.

Figure 9

 9 Figure9.3: Cuboid-like boundary conditions on which the U-Net is pre-trained in strategies 3 and 5. In strategy 6, there are four more cuboids so that the lower part of the beam is fully covered. In strategy 7, the U-Net is pre-trained on BC1.

Figure 9 . 4 :

 94 Figure 9.4: Randomly visible sub-domains of the upper face of the beam (in yellow).

  N., Nikolaev, S., Cotin, S. (2020) Data-driven simulation for augmented surgery. Advanced Structured Materials, vol 132, Springer Towards Patient-Specific Networks 116 Towards Patient-Specific Networks 117

Figure 10 . 1 :

 101 Figure 10.1: Modéle simplifié de la sclère.

Figure 10

 10 Figure 10.3: Architecture de réseau de la plateforme U-Mesh.

Figure 10

 10 Figure 10.4: Comparaison de la prédiction U-Mesh (en vert) et de la solution EF (en bleu) pour une poutre encastrée et une lettre L.

Figure 10 . 6 :

 106 Figure 10.6: Nuage de points RGB-D donné en entrée au réseau et prédiction résultante.

Figure 10

 10 Figure 10.7: Prédiction U-Mesh en réponse à la compression induite par la sonde US.

  D is a constant. Also surface forces f s are applied on a subset Γ N of the boundary Γ (the so-called Neumann conditions) such that σ

	throughout this manuscript, we can set the second derivative of the displacement to
	zero. Hence, Equation (2.21) reads as:
	-div(σ) = f ext .	(2.22)
	To ensure the uniqueness of the solution to this equation, we must impose Dirichlet
	boundary conditions to function u on a subset Γ D of the boundary of the domain Ω
	such that
	u(x, t) = u D , ∀x ∈ Γ D ,	(2.23)
	where u
	Biomechanical Problem Formulation 15

The Equation (2.21) is known as the strong formulation. It represents the fundamental principle of dynamics at each particle of the body. As we will only consider static problems

  D be the subpart of Γ on which Dirichlet conditions are applied and Γ N the subpart on which Neumann conditions are applied (such that Γ D ∩ Γ N = ∅). Dirichlet conditions specify the value of the function u on the domain boundary, while the Neumann conditions impose a traction at the boundary. In the first part of this manuscript, we assume the Dirichlet conditions on Γ D are known a priori, while Neumann boundary conditions on Γ N can change at any time step. In the context of soft tissue simulation, the correct estimation of boundary conditions is a non-trivial problem

	.22), boundary
	conditions must be imposed on the unknown function u. The boundary conditions
	describe the effect of the exterior environment on the considered object. There are
	two types of boundary conditions: Dirichlet and Neumann. We recall that Γ is the
	boundary of the domain Ω. Let Γ

  Table 4.1: Error measures for a beam having 135 H8 elements. Rows are sorted in decreasing e . The selected architecture appears in bold.

	c	k	FSS	e	σ(e)	pred_t train_t
				in m	in m	in ms	in min
	64	2	256	0.0028	0.0008	2	24
	16	4	256	0.0012	0.0009	3.2	40
	32	4	512	0.0009	0.0008	3.2	40
	64	3	512	0.0007	0.0007	2.5	80
	64	4 1024	0.0007	0.0007	3.3	95
	128 3 1024 0.0007 0.0006	2.5	35
	64	5 2048	0.0006	0.0005	4	426
	128 4 2048	0.0006	0.0005	3.45	320

U-Mesh Framework 44

Table 4

 4

		Prediction time	e
		in s	in m
	FEM	0.740	0.000
	HPOD	0.005	0.084
	U-Mesh	0.003	0.006

.6: Comparison at similar computation time: the mean norm errors of the full FEM model, of the HPOD and of the U-Mesh are given, for a computation time in the range of the millisecond. Only one mode is kept in the HPOD.

The U-Mesh Framework 51

  • Image-based force estimation in robotized procedures (no force sensor needed)• Training data generated with a simplified finite-element model of the eye • Identification of the stiffness of the sclera using Kalman filters • 93% accuracy in the force classification on real porcine data

	• Publications:
	-Mendizabal, A., Fountoukidou, T., Hermann, J., Sznitman, R., Cotin,
	S., 2018. A Combined Simulation and Machine Learning Ap-
	proach for Image-Based Force Classification During Robotized
	Intravitreal Injections. In International Conference on Medical Image
	Computing and Computer-Assisted Intervention, pp. 12-20.
	-Mendizabal, A., Sznitman, R., Cotin, S., 2019. Force classification
	during robotic interventions through simulation-trained neural
	networks. International journal of computer assisted radiology and sur-
	gery, 14(9), pp. 1601-1610.

Force Estimation in Robotized Injections 78

Force Estimation in Robotized Injections 79

Table 8 .

 8 1: Error measures over the testing data set for a breast having 2174 H8 elements, with maximal nodal deformation of 79.09 mm.

	The U-Mesh for Ultrasound-guided Breast Biopsy 98

Table 8 .

 8 2: Target registration errors in millimeters for different tumors and different deformation ranges in the breast phantom. The first table is for the proposed method, while the second table reports results obtained with the FE model used for data generation. Not-acquired data is reported as (-).

			U-Net predictions
	TumorID	D10	D15	D20	D25	D30	Mean	STD
	1	-	1.936 2.002 1.506 3.053	2.124	0.569
	2	3.211	2.905 4.068	-	4.137	3.580	0.534
	3	2.032	-	4.709 7.134 10.90	6.194	3.262
	4	0.505	2.225 5.313 5.903	-	3.486	2.217
	5	0.932	2.768 3.454	-	4.893	3.012	1.425
	6	3.923	6.349 5.625	-	6.724	5.655	1.075
	7	3.454	3.864 4.543 6.710	-	4.643	1.255
	8	2.422	3.261 4.320 5.136	-	3.785	1.030
	9	-	3.928 4.214 4.578 4.858	4.394	0.353
	10	5.529	3.272 3.940 4.846	-	4.397	0.860
	Mean	2.751	3.390 4.219 5.116 5.761
	STD	1.638	1.294 1.007 1.854 2.788
				FE method
	TumorID	D10	D15	D20	D25	D30	Mean	STD
	1	-	1.326 2.151 2.075 3.759	2.328	0.887
	2	1.956	2.738 3.945	-	4.025	3.166	0.865
	3	1.595	-	4.748 7.044 10.932	6.080	3.404
	4	0.755	1.991 4.544 5.120	-	3.103	1.795
	5	1.029	2.863 3.330	-	4.541	2.941	1.262
	6	2.579	3.409 2.871	-	2.337	2.799	0.400
	7	2.605	3.219 4.095 6.750	-	4.167	1.582
	8	2.695	2.748 4.321 5.411	-	3.794	1.139
	9	-	2.745 2.497 2.510 4.193	2.986	0.704
	10	2.916	2.542 3.015 3.868	-	3.085	0.485
	Mean	2.016	2.620 3.552 4.682 4.964
	STD	0.765	0.593 0.856 1.803 2.757

Table 9 .

 9 2: Error measures over all seven scenarios. Best iterations are given in thousands. Transfer learning situations are highlighted in red, first and second best results in green and blue.

	Towards Patient-Specific Networks 112

Part III Brief summary in french 123 Titre: Combinaison entre simulation numérique et apprentissage auto- matique -Applications à la formation, la planification préopératoire et l'assistance peropératoire 10.1 Introduction Dans

  un contexte de chirurgie augmentée, il est nécessaire de réaliser un recalage élastique des données préopératoires sur la vue intraopératoire de l'organe. L'objectif est de superposer en temps réel un modèle virtuel à la vue du champ opératoire, afin de visualiser les structures internes de l'organe (tumeurs, réseau vasculaire etc) qui ne sont pas naturellement visibles pour le chirurgien. Le geste chirurgical est ainsi mieux maîtrisé et l'intervention peut être planifiée de façon plus adaptée. Faire ce recalage dans le cas des tissus humains et avec des données intraopératoires très clairsemées, relève plusieurs défis. Par exemple, dans un contexte de chirurgie hépatique, seulement 30% de la surface du foie est visible à cause de la taille de l'incision ou le champ réduit de la caméra laparoscopique[START_REF] Plantefeve | Patient-specific biomechanical modeling for guidance during minimally-invasive hepatic surgery[END_REF].

	General conclusion 121 General conclusion 122

  la déformation de la cage thoracique est prédite par un RN. Dans ces deux papiers, seul le déplacement surfacique est estimé.Morooka et al. propose les RN pour estimer la déformation volumique d'un foie. Il utilise une ACP (analyse en composantes principales) pour compresser les modes de déformation des sorties du réseau dans le but d'accélérer l'entrainement. Leur model produit des résultats satisfaisants mais a uniquement était testé sur des données simulées.Afin de répondre aux contraintes de précision et de temps-réel, on suggère de combiner des simulations EF avec des RN pour modéliser le comportement biomécanique des tissus humains. Cette idée s'intègre très bien dans un contexte de chirurgie guidée par l'image, puisqu'il est ainsi possible de collecter beaucoup d'information lors d'une intervention, et ainsi d'apprendre continuellement à améliorer la modélisation. En particulier, nous proposons la plateforme U-Mesh qui permet de prédire en tempsréel les déformations d'un organe hautement déformable comme le foie avec une précision permettant de guider les chirurgiens au cours des interventions où il est essentiel de suivre cette déformation (par exemple lors d'une biopsie). Les réseaux sont entrainés avec une très grande quantité de données générées avec la méthode des EF. Pour garantir des modèles patient-spécifiques, on combine simulation et inférence bayésienne afin d'identifier les paramètres des modèles EF.

Biomechanical Problem Formulation 16

Biomechanical Problem Formulation 21

Artificial Neural Networks: Fundamentals 25

Artificial Neural Networks: Fundamentals 27

Artificial Neural Networks: Fundamentals 28

Artificial Neural Networks: Fundamentals 31

The U-Mesh Framework 50

The U-Mesh Framework 52

www.slicer.org Patient-specific modeling 58

Patient-specific modeling 61Patient-specific modeling 62

Force Estimation in Robotized Injections 64

Force Estimation in Robotized Injections 72

https://www.sofa-framework.org/ Force Estimation in Robotized Injections 74

Force Estimation in Robotized Injections 75

The U-Mesh for Augmented Hepatic Surgery 85

The U-Mesh for Ultrasound-guided Breast Biopsy 95

The U-Mesh for Ultrasound-guided Breast Biopsy 101

General conclusion 120

Acknowledgements

"They did not know it was impossible so they did it" quote by Mark Twain.

Contributions of this chapter:

• The U-Mesh framework: nonlinear dimensionality reduction for displacement field estimation trained with controlled finite-element simulated data • Predictions on two benchmark examples (cantilever beam and L-shaped structure) with very high accuracy in only 3 ms • Constant computation time regardless of the size of the problem • Publication: [START_REF] Mendizabal | Simulation of hyperelastic materials in real-time using deep learning[END_REF]. Simulation of hyperelastic materials in real-time using deep learning. Medical Image Analyses. Volume 59, 101569 10.1016/j.media.2019.101569

The U-Mesh Framework 53

The U-Mesh Framework 54

Part II

Patient-specific real-time simulations

55

Contributions of this chapter:

• U-Mesh applied to simulate breast deformation • Data generation strategy adapted to the studied phenomena (simulation of probe-induced displacements) • Accurate estimation of target location (∼ 4 mm error) in real-time • Publication: [START_REF] Mendizabal | Physics-based Deep Neural Network for Real-Time Lesion Tracking in Ultrasound-guided Breast Biopsy[END_REF]

Résumé

Il existe une multitude de domaines en ingénierie nécessitant le calcul de déformations non-linéaires en temps réel, notamment dans le domaine de la médecine, avec les simulateurs pour l'entrainement des chirurgiens ou bien la chirurgie guidée par l'image. Dans un contexte de chirurgie augmentée, il est nécessaire de réaliser un recalage élastique des données préopératoires sur la vue intraopératoire de l'organe. L'objectif est de superposer en temps réel, un modèle virtuel spécifique à chaque patient, à la vue du champ opératoire dans le but de visualiser les structures internes de l'organe (tumeurs ou artères par exemple). Afin d'obtenir une précision adaptée, il est nécessaire de construire un modèle bioméchanique personnalisé de l'organe. Pour ce faire, la méthode des éléments-finis est la technique préférée afin de prédire la déformation des tissus mous. Cependant la complexité des calculs impliqués (notamment dans le cas de déformations nonlinéaires) rend cette méthode incompatible avec les exigences de temps réel et de précision, inhérentes au domaine d'application visé.

Pour répondre à ces contraintes, on propose de combiner des simulations éléments-finis avec des réseaux de neurones profonds pour modéliser le comportement biomécanique des tissus humains. En particulier on présente la plateforme U-Mesh permettant de prédire en temps-réel les déformations d'un organe comme le foie avec une précision adaptée à la chirurgie augmentée.

Mots clefs : Simulation temps-réel ; Simulation patient-spécifique ; Réseau de neurones pronfond ; Méthode des éléments-finis ; Chirurgie augmentée

Résumé en anglais

Many engineering applications require accurate numerical simulations of non-linear structures in real-time. Some important examples can be found in the field of medicine, in order to develop surgical training systems, or in the field of surgical navigation where augmented reality can bring significant improvements to the clinical gesture. To guarantee the accuracy of the simulations, patient-specific modeling must be pursued by taking into account personalized material parameters and boundary conditions. In the context of augmented surgery for instance, it is essential to perform an elastic registration between the preoperative and the intraoperative images. To this end, a patient-specific biomechanical model must be built to produce real-time finite-element simulations of the deformed organ. This is in practice very difficult to achieve as the problems to be solved are highly complex, in particular when non-linear deformations are considered.

In this work, we propose a method combining finite-element simulations and deep neural networks in order to satisfy the rapidity and accuracy requirements of medical applications. In particular, we present the U-Mesh framework, capable of predicting in real-time the shape of a highly deformable organ like the liver in order to guide surgeons during interventions where following the organ's deformation is crucial for the surgery to be successful.

Keywords: Real-time; Deep neural networks; Finite element method; Data-driven simulation; Augmented surgery