Jérémie Allard

Christian Duriez

Cotin, who gave me the passion for this fascinating job. I also think of Stéphane Bordas, who inspired me to push myself and reach my

Introduction

T he end of the 1980s marked a turning point since which Minimally Invasive Surgery (MIS) has been progressively recognized as a safe and effective approach to meet surgical needs while decreasing the rates of complications and traumas for patients. In less than 20 years, it has become the gold standard for many practitioners and has completely revolutionized how surgery is approached in developed countries. On the other hand, despite the numerous benefits associated with patients, the paradigm shift from open surgery towards MIS has significantly raised the technical level required to perform surgery. One of the main reasons lies in the fact that MIS adds distance between surgeons and patients, thus preventing direct visual contact with organs, with poor access conditions and limited tactile sensations. As a result, the practice of surgery has evolved into an extraordinarily multidisciplinary and complex science that is now increasingly assisted by tools and devices to improve the quality of care.

Imaging systems have become central for diagnosis, but also for assisting practitioners during the interventions (Figure 1.1). These reasons explain why sciences such as physics, biochemistry, electronics, and computer science today contribute to the quality of patient care.

Currently, the training of medical personnel is carried out first by observing other more experienced practitioners. Although practice remains the best solution to learn coordination between eyes, hands, and instruments, it mainly remains possible on human subjects. Indeed, during students' residency, practice on animals, cadavers, or samples is usually conducted, while posing economic and ethical issues. Later in their curriculum, students are trained with cadavers or samples (porcine or bovine. . . ) and then with anesthetized animals. However, since the geometries and material properties of the tissues are different and rapidly change after death [START_REF] Abolhassani | [END_REF], it may lead to several inaccuracies with living patients as well as a lack of representative materials for proper training. As a result, the full breadth of training represents more than ten years due to the difficulty accessing internal structures and anticipating the large panel of possible complications. These limitations in training time and the difficulty to access surgical practice may lead to an inadequate preparation in a context where a series of worldwide directives aims at optimizing the use of medical resources.

Interactive medical simulations

Medical simulators have multiple benefits over traditional training methods, such as the possibility to train physicians in various scenarios, including common and rare cases, or to change the properties of the tissues in a repeatable manner [START_REF] Taschereau | [END_REF]. In addition, trainees do not require the direct supervision of trained clinicians, thus saving time for the medical staff without compromising patients' safety. However, medical simulators must provide a realistic and immersive experience and have enough flexibility to allow users to make mistakes. When it happens, it shall not be the end of the scenario. Indeed, the novice must also learn to recover from the mistakes since errors made at one point in the procedure will impact the following steps. Physics-based simulations can answer these requirements to provide high-fidelity scenarios and simulate surgeons' interventions in a virtual environment.

Following the development of MIS, computed-based simulators have received a considerable interest both in academics and industry Bro-Nielsen and Cotin (1996); Müller and Gross (2004); Courtecuisse et al. (2014aCourtecuisse et al. ( , 2015)); [START_REF] Nguyen | [END_REF]. Indeed, MIS allows reproducing realistic conditions in a virtual simulator: the rendering of synthetic images is now close to authentic medical images, and because interactions with organs (only possible through instruments) are limited they can be rendered with haptic devices. In addition, advanced numerical simulation allows for free and theoretically unlimited user interactions with virtual organs with many variabilities in terms of scenarios to match patient-specific conditions (geometry, accessibility of the tumors, physiological motions of the patient, and pathologies). The physics-based simulation could therefore open new avenues leading to improved training and more general care as well as risk reduction for patients [START_REF] Maier-Hein | [END_REF].

Computer Assisted Interventions

Medical imaging is now one of the essential aspects of in order to perform routine daily surgery. It represents not only a support for diagnostics but also an actual operative instrument during particular therapeutic procedures. There are many types of medical imaging modalities based on different technologies and techniques (Ultrasound, X-ray, MRI, fluoroscopy,...), each presenting benefits and drawbacks Seibel (1997). For instance, tomography techniques (CT and MRI) provide detailed anatomical representations of highly contrasted organs. CT scans are fast, but they involve X-rays absorption, whereas MRI increases the length of the procedure and requires non-ferromagnetic equipment.

Pre-operative images are usually acquired for diagnosis and surgical planning. Once the surgical plan has been created, the next step is to guide the surgery. However, from a physiological point of view, organs are not static, and even the most straightforward breathing motion can induce shape deformations. These modifications may invalidate pre-operative planning since the location of internal structures may significantly vary. Identifying internal structures and planning during the operation requires expertise and interpretation of intraoperative data and a mental effort of three-dimensional reconstruction to coincide with the pre-operative images. As a result, even expert surgeons request for new technologies for guidance and assistance during the intervention.

Nowadays, several levels of assistance are considered by the scientific and medical communities:

Visual assistance: Digital tools are more and more used for intraoperative assistance. The goals are, for instance, to display through Augmented Reality (AR) the internal structures (vessels, tumors. . . ) on top of the intraoperative images that are often sparse, incomplete, and of poor quality. An essential advantage of biomechanical models lies in their ability to predict the behavior of structures, providing this way a physics-based extrapolation, not just geometric, in areas where few or no intraoperative data are available.

Gesture Assistance: Robots were introduced in the operating rooms to overcome limitations related to human factors and thus assist surgeons in tasks requiring high accuracy. The main objectives are accuracy, speed, repeatability, reduction of efforts, remote action for hostile environments (Xray), and filtering of movements and tremors. To this end, an increasing number of surgical robots have been developed and applied to a wide range of surgical applications Troccaz (2013). However, structure deformations remain significant obstacles to access advanced algorithms and automation of complex tasks.

Percutaneous procedures

Percutaneous image-guided interventions are an essential example of the evolution of the surgical practice. Needle-based interventions are among the least invasive surgical approaches to access deep internal structures into organ volumes without damaging surrounding tissues. Percutaneous procedures have become the gold standard for numerous interventions such as biopsy, vertebroplasty, or radiofrequency ablation (RFA) that uses heat at the tip of a needle to destroy cancer cells. Unlike traditional invasive approaches, needles only affect a localized area around the needle shaft, reducing this way the occurrence of traumas and risks of postoperative complications Cowan et al. (2011). However, contrary to what one may think, needle-based interventions can be exceedingly complex (Figure 1.2). Indeed, the effectiveness of the treatment is highly dependent on the accuracy of the needle positioning (≈ 3 mm De [START_REF] Jong | [END_REF]), which can be particularly challenging when needles are manipulated from outside the patient using medical images (X-ray, fluoroscopy, or ultrasound) offering poor visibility of internal structures. Therefore, the possibility for patients to benefit from these alternative treatments depends first on the medical recommendations (concerns about the age, the extent, and the localization of the disease) and the practitioner's experience and preferences concerning needle-based approaches.

Human factors, organ deformations, needle deflection, and image limitations may lead to needle misplacement and raise significantly the technical level needed to perform these surgical acts. In addition, the needle's path may cross sensitive tissues such as nerves, bones, arteries, or organs that may lead to complications if damaged Chentanez et al. (2009). Concerning hepatic surgery, tumors can be located deeply inside the liver with poor access conditions (for instance, in posterior segments VI or VII), thus reducing the ability to control the needle tip close to the target. In addition, the liver is a very soft organ that tends to deform in contact with the diaphragm during respiratory motion resulting in non-uniform displacements (up to 5 cm).

The previous difficulties explain why interventional radiologists usually need several attempts to reach the targeted zone or need to operate under apnea Poon et al. (2004), which increases intervention time and risks for the patient.

Contributions and medical applications

The main application followed in this manuscript is to develop an autonomous but supervised robotic system for needle steering in percutaneous liver surgery. This project is motivated by the fact that FE simulators are now considered clinically relevant tools for training, formation, and surgical assistance. Thus, we believe that a sufficient knowledge of the mechanical behavior of structures can be extracted from FE models to guide a robot and reach a target, pre-identified at planning stages. Our approach consists in developing control models to correct and predict the deformation of structures and adapt the behavior of the robot to reach an objective. Based on a tomographic image, we acquire the geometry of the organ and plan the trajectory for the needle insertion. Intraoperatively, the robotic system automatically performs the needle insertion. FE models, including needle/tissue interactions and collisions with organs, are used to predict the behavior of internal structures. To limit errors, FE models are permanently registered with intraoperative data extracted from medical images. Input commands of the robot are computed automatically in real-time with an inverse simulation using the registered FE models. In addition, for safety reasons, the user shares the control between the automatic method and remote controllers.

Both the needle and the tissue undergo deformations (due to physiological motion or interactions with the needle), and the targets, trajectories, and obstacles may move simultaneously. Furthermore, the robot cannot be considered independent from its environment once the needle is inserted into the organ. Therefore, the robot must compensate or even generate deformations of tissues to align the shaft of the needle with the desired trajectory. For this purpose, we proposed to derive the robotic commands from 3 essential steps that lie at the intersection of several scientific domains (Figure 1.3):

1. Numerical models: The first step is to develop numerical methods for dynamic simulations of soft tissues, instruments, and their interactions in the virtual environment. The simulations must guarantee at the same time stability and accuracy while maintaining low computational expense to enable interactivity. Specific models accounting for complex phenomena (such as friction or heterogeneity) are necessary to increase the predictive nature of the simulation. Advanced numerical solvers exploiting time/space coherency and parallelism are introduced to solve the equations of motion, maintaining interactive computation time without scarifying the accuracy.

2. Non-rigid registration: An important question that will be addressed is to combine FE simulations with image-based constraints while allowing sufficient freedom to perform accurate predictions. Indeed, as an approximation, FE models cannot be used in an open-loop and must be corrected with intraoperative data (image, sensors). The advantages of relying on FE models for this task are multiple. FE models can be used to regularize the data extracted from images, being this way less sensitive to noise and tracking errors. In addition, FE models provide accurate solutions to interpolate the whole volume displacement of the organs, including the displacement of internal structures such as tumors or vessels.

3. Robotic control: The simulation will be used to locally predict, in the virtual environment, the behavior of the structure under specific motions of the robot. The critical scientific challenge is to provide alternative solutions to control a robot based on inverse Finite Elements (iFE) simulations. In this case, inverse simulations provide a way to anticipate the behavior of soft structures allowing in this way to adapt the robotic commands much faster than waiting for visual correction from images. However, when the robot interacts with soft tissues, it produces a highly nonlinear system, and the validity domain of the simulation is limited. In order to maintain the model's consistency, simulations are recomputed at a high frequency in the control loop of the robot. Thus, a significant challenge is to propose numerical solutions to solve iFE simulations sufficiently fast to enforce the stability of the control loop.

The research project introduced in this document is part of the logical evolution of surgery to improve safety, recovery time and offer less aggressive treatments for patients. Objectives aim first at improving the accuracy, stability, and realism of numerical simulations while maintaining a performance compatible with user interactivity. Based on these advanced and fast FE models, we introduce essential contributions to assist surgeons during interventions, providing in this way both visual and gesture assistance.

Outline

The remainder of this document is organized as follows. Chapter 2 is dedicated to the description of numerical methods for fast simulation of deformable structures with contacts. Chapter 3 concerns the non-rigid registration of soft tissues with intraoperative data. Chapter 4 introduces our solution for automatic needle steering using inverse Finite Element simulations. Each chapter is organized with the following homogeneous structure: After a short introduction of the context, we briefly review the related works from the literature. A generic methodology is then introduced, describing the general background and scientific challenges. Then, we introduce our contributions on the domain. The chapters end with the perspectives on the topic and a conclusion.

Finally, Chapter 5 gives an overview of the research projects funded and other related activities such as technological transfers, supervision, and teaching activities.

Introduction

Medical simulations are becoming recognized as a method for continuous professional development for medical personnel. Compared to the traditional formation, trainees can dive into complex scenarios or critical situations that otherwise are hard, expensive, and risky to replicate in real life. It means that trainees will be ready to go through the whole set of situations and get much more experience. Recent advances in the field open the possibility to use the simulations not only for training but also for pre-and intraoperative support (Figure 2.1). Indeed, virtual simulations of anatomical structures could provide strong support for the diagnostic phase, but also during the intervention in combination with augmented reality systems (see chapter 3).

For all the above applications, the obvious requirement is to simulate the mechanical response of organs with high accuracy. However, it is also important to note that significant computation time constraints always exist. Without meeting this requirement, the simulator cannot be used in an interactive context, regardless of the quality and accuracy of models. Therefore, the underlying contributions of this chapter are a set of advanced numerical techniques, including modeling, constraints definition, and high parallelization strategies to meet these requirements.

The main requirement to build useful surgical simulators is to reproduce the mechanical response of organs. For that, surgical simulators involve at least the following major challenges:

1. The simulation of complex interactions is necessary. The user must have the possibility to interact with virtual models and ideally to get haptic feedback from the scene.

2. It is necessary to consider the deformable environment to obtain a realistic global behavior of the tissues. In this case, boundary conditions are also fundamental because stiff interactions with the tool (controlled by the user) and complex interactions between deformable structures must be solved simultaneously.

3. The response must be computed in real-time to allow user interactions with the virtual body: pushing, prodding, palpation and cutting. . . It always requires a balance between accuracy/speed and stability.

4. Significant numerical issues must be overcome, as most organs have a heterogeneous stiffness (which leads to ill-conditioned problems) and are composed of multiple tissue types.

On top of the above requirements, it is also essential to enforce the stability of the simulation, in particular in response to user interactions. Indeed, although many solutions in the literature only focus on the accuracy of FE solutions, stability is crucial in a medical context.

Related Works

Biomechanical simulations with user interactions involve many challenges such as real-time computation of the deformation of soft tissues, collision detection, contact modeling, topological modifications, and haptic feedback (see [START_REF] Nealen | [END_REF]; Konofagou (2012) for a broad survey).

Simulation of deformable bodies

The first methods proposed to simulate the deformation of soft tissues in real-time relied on mass-spring systems, e.g., [START_REF] Kühnapfel | [END_REF]. Although such discrete methods are simple to implement and fast, they are challenging to parameterize with material properties such as Young's modulus. Moreover, they introduce anisotropy through the choice of the mesh giving rise to stability and accuracy issues (node flipping, difficulty to preserve the volume).

Finite element methods (FEM) provide high biomechanical realism Kirby and Logg (2012), mainly because the complex nonlinear behavior of soft tissue is directly accounted for through constitutive relations. Real-time computations were first achieved for linear elastic material models (see Bro-Nielsen and Cotin (1996); Cotin et al. (1999) or James and Pai (1999)). In linear elasticity, precomputations (offline) can be used to accelerate the online simulations. However, the small strain assumption is incorrect in practice and produces erroneous results when the solids undergo large deformations. The use of pre-computed solutions for highly nonlinear problems is intensively pursued, for instance, in Niroomandi et al. (2008) for hyperelasticity.

The co-rotational method, originated in continuum mechanics Freund (1970), was introduced by Felippa (2000) within the field of numerical methods. In this formulation, the stiffness of each element is assumed linear within the local frame described by its rotated state. This approach simulates geometric nonlinearity (i.e., large displacements and rotations) while maintaining low algorithmic complexity. However, the formulation is restricted to small strain ranges, excluding important physical aspects such as volume preservation during deformations. The formulation is still widely used in the computer graphics community. For instance, Frâncu et al. (2019) introduced a novel method of simulating incompressible materials undergoing large deformation without locking artifacts.

Later, other methods were proposed to simulate both geometric and material nonlinearity in real-time. The Total Lagrangian explicit dynamics (TLED) Miller et al. (2007) allows modeling non-linearity, and fast GPU-based implementations were proposed [START_REF] Comas | [END_REF]; Joldes et al. (2009). However, the main limitation is that the method is based on an explicit time integration scheme, entailing small time steps to keep the computation stable, especially for stiff materials. Marchesseau et al. (2010) proposed the Multiplicative Jacobian Energy Decomposition (MJED), a general algorithm to implement hyperelastic materials with implicit time integration schemes. However, although some validations of the behavior against real organs have been conducted, these models remain complex and expensive, and the simulation of realistic boundary conditions such as interactions between deformable organs and surgical instruments is still an issue.

Meshless, Position-Based dynamics (PBD), and Neural Networks are other strategies to model soft tissues in real-time. A detailed review of this topic goes far beyond the scope of this document, but a survey can be found in [START_REF] Zhang | [END_REF].

Interactive simulation and numerical solver 2.2.2.1 Time discretization

In the context of interactive simulations, an important choice is the time integration scheme. Explicit integration schemes are widely adopted Miller et al. (2007); Taylor et al. (2008) in surgical simulations. In this case, the solution only involves the (diagonalized) mass matrix leading to very fast, simple to implement, and parallelizable solutions [START_REF] Comas | [END_REF]. The equations of motion are decoupled, and each degree of freedom can be solved independently, making the solution process very fast and inherently well suited to parallelization [START_REF] Comas | [END_REF]. The major drawback of explicit dynamics is the need to satisfy the Courant-Friedrichs-Lewy stability condition, which forces a strict upper bound on the time step used for integration. Explicit methods are consequently particularly well suited to very soft tissues such as the brain Joldes et al. (2009), but very small-time steps (which prevent realtime computations) must be chosen for stiffer structures. Moreover, explicit simulations do not guarantee that the residual vector is minimized at each step, so the external and internal forces are balanced.

On the contrary, implicit methods are unconditionally stable 1 for any time step and arbitrary stiff materials Baraff (1996). Implicit schemes provide better control of the residual vector and hence that the external and internal forces are balanced at the end of the time steps. Of course, the advantages of such methods come at the cost of having to solve a set of linear equations at each time step. Nevertheless, this document shows that implicit integration schemes can offer a reasonable tradeoff between robustness, stability, convergence, and computation time, particularly when combined with a GPU implementation.

Solving the nonlinear equations

The nonlinear problem obtained by implicit schemes is usually solved using an iterative Newton-Raphson method. Each iteration of the Newton method consists of solving a linear problem whose solution reduces the error between internal and external forces. The mechanical matrices depend on the position of FE meshes (and potentially velocities) except for the linear model. The linear problems must then be recomputed for each simulation step and ideally for each Newton iteration. Therefore, the simulation performances are directly linked to the efficiency of the solver, which explains why earlier studies mainly focused on sparse linear algebra. The linear equations can either be solved by direct solvers or iterative solvers.

Direct solvers: Direct solvers provide the exact solution by computing a factorization (for instnace, the Cholesky factorization [START_REF] Barbič | Real-time subspace integration for St.Venant-Kirchhoff deformable models[END_REF]) or a decomposition (QR decomposition), or eventually, the actual inverse of the system matrix Bro-Nielsen and Cotin (1996) (though not recommended for large matrices). As proposed in George (1973), the nested dissection ordering has been widely used in direct solvers by exploring the parallelism of subproblems while reducing the fill-in of the matrix. Partitioning and reordering are usually implemented through external tools such as METIS Karypis and Kumar (1998).

Explicit assembly of global matrices is necessary for direct solvers to compute the factorization or decomposition of the system. There are several ways to construct sparse matrices; the most popular method is first to collect triplets (the row/column index and the value); then compress the triplets in a sparse format. A very efficient implementation is provided in the Eigen library. Hiemstra et al. (2019) proposed a row by row assembling method for isogeometric linear elasticity problems. Several approaches proposed to assemble the matrix directly on the GPU and minimize memory transfers Dziekonski et al. (2012); Zayer et al. (2017); Fu et al. (2018).

The solving phase can then be performed with so-called forward substitution, using the two triangular systems (in the case of the Cholesky factorization). Several efficient libraries exist on the CPU (Pardiso, MUMPS, Taucs) and GPU (cuSPARSE, MAGMA, AmgX). In addition, the solving stage can be improved by partitioning and reordering the system Herholz and Alexa (2018).

Despite the stability of direct solvers, the complete factorization or decomposition of large matrices is usually too time-consuming to be recomputed at each time step, and it is complicated to parallelize. Specific optimizations, inspired by the co-rotational model, have been proposed to incrementally update the sparse Cholesky factorization Hecht et al. (2012) but this approach does not extend to other material laws or element types. Other approaches aim at reducing the number of degrees of freedom of the model, either using condensation on surface nodes Bro-Nielsen and Cotin (1996), or reducedcoordinate models [START_REF] Barbič | Real-time subspace integration for St.Venant-Kirchhoff deformable models[END_REF].

Iterative solvers: In the interactive context, iterative methods are usually preferred because they offer the possibility to compute an approximated solution and limit the number of iterations to better control the time spent during the solving process. The most popular method is the Conjugate Gradient (CG) algorithm Saad (2003), because of the fast convergence and its simple implementation. Parallel solutions both on CPU Parker and O'Brien (2009); [START_REF] Hermann | [END_REF] and GPU Bolz et al. (2003); Buatois et al. (2009); [START_REF] Allard | [END_REF] have been proposed.

Iterative solvers are usually faster than direct methods and require less memory storage, but they converge slowly for ill-conditioned problems, i.e., when the ratio of the largest and smallest eigenvalues is large. This is the case when solving linear systems of equations associated with heterogeneous structures.

Preconditioner: Another intense area of research aims to improve the CG algorithm's performance using preconditioners to speed up its convergence.

In the context of interactive simulation, [START_REF] Baraff | [END_REF] proposed to use a diagonal inverse, often called a Jacobi preconditioner. More advanced preconditioners such as Cholesky factorization have also been studied Hauth et al. (2003). In the real-time context, one substantial limitation of this technique is the computational overhead added to the simulation: first, during the computation of the preconditioner itself, and second, during its use at each iteration of the CG (see Saad (2003) for details). Thus, the practical usefulness of preconditioners depends on the ability to find a balance between the computational overheads and the time saved by decreasing the number of CG iterations.

Simulation of the interactions

A key requirement for realistic surgical simulators is to treat contact between soft-soft, soft-stiff and stiff-stiff objects. The simulation of interactions is a vast topic. Pioneer work in the domain was proposed by Jean and Moreau (1992).

A standard solution for contact simulation consists of using a penalty method, which modifies the variational principle and solves the contact condition approximately. A contact force is added at each contact point proportional to the interpenetration and a user-defined penalty coefficient. The higher the coefficient is chosen, the better the constraints are satisfied, but the worse the condition number is affected. The choice of the penalty coefficient is often problem-dependent and depends strongly on the stiffness ratio between the contacting objects.

Lagrange multipliers or augmented Lagrangian techniques Hughes (2008); Renard (2013); Jean (1999); Jourdan et al. (1998); Peter Wriggers (1999) are usually preferred over penalty methods to treat contact constraints accurately and robustly. Methods used to solve contact equations with Lagrange multipliers can be classified into two categories: Quadratic Programming (QP) methods and Complementary Problem methods that could be linear (LCP) or nonlinear (NLCP).

QP methods define the constraints directly into the mechanical system. They can be used to treat the inequality of the contact Redon et al. (2002); Pauly et al. (2004), and also to simulate friction using a discretized pyramidal cone Kaufman et al. (2008). However, these publications address the case of rigid bodies in contact, where the number of degrees of freedom (DOF) is smaller than the number of contact constraints. Indeed, the number of equations resulting from QP methods is the same order as the number of DOFs of the interacting objects. Therefore, these methods are difficult to adapt for the real-time simulation of interactions between FE meshes of large dimensions.

Various numerical methods can be used to solve LCPs in physics-based an-imations Erleben (2013). Direct methods such as pivoting methods give an exact solution, but they are not computationally efficient. In contrast, iterative methods have been more widely applied in large-scale simulations, especially for real-time computations. Being simple to implement, projected Gauss-Seidel (PGS) [START_REF] Duriez | [END_REF], [START_REF] Courtecuisse | [END_REF], Macklin et al. (2016)) could handle the friction response with the Coulomb's friction cone combined in the LCP formulation. However, the algorithm is not efficient for ill-conditioned problems due to the slow convergence. Macklin et al. (2019) proposed using a Newton method to solve the nonsmooth functions that are reformulated from complementarity problems. A complementarity preconditioner accelerates the convergence of a conjugate residual (CR) algorithm for the constraint resolution. Although very robust and efficient for interactive simulations, the method is introduced in a context with the number of DOFs is smaller than the number of constraints. LCP can be used to simulate frictionless contact between deformable models Duriez and Andriot (2003) in real-time, whereas NLCP methods can be used to simulate friction contact with the exact friction cone [START_REF] Duriez | [END_REF].

An essential advantage of (N)LCP methods is that the number of constraint equations is proportional to the number of contacts, which is often much smaller than the number of DOF in the context of deformable models. However, the main difficulty is related to the fact that the solution process involves the compliance matrix, which is the inverse of a large system composed of the mass, the damping, and the stiffness of the deformable objects. Duriez et al. (2004) proposed to precompute the compliance matrix, but the solution is limited to linear elastic deformation. Otaduy et al. (2009), proposed to compute the compliance matrix using additional Gauss-Seidel iterations on the deformable models, but the method was not presented in a real-time context. Saupin et al. (2008b) proposed a method, named compliance warping, which is dedicated to co-rotational models. It consists of pre-computing the compliance matrix from the rest position and updating it using a local estimation of the nodal rotations. However, this approximation can become inaccurate for large deformations, and the method is limited to relatively coarse meshes. A prediction-correction scheme is introduced in [START_REF] Peterlík | Constraint-based haptic rendering of multirate compliant mechanisms[END_REF] to mitigate the inaccuracies, and the method extends the formulation to haptic feedback with generic constraints between the deformable models.

Simulation of topological modifications

When cutting through a finite element mesh, discontinuities in the displacement field must be introduced. Methods can be classified into three parts, namely the mesh-based methods, the meshless methods, and hybrid methods [START_REF] Wang | [END_REF].

Mesh-based methods simulate the deformations using FE and the cutting process with topological changes. Among those changes, the easiest one is the element deletion, which consists in deleting an element that crosses the cutting path Forest et al. (2005). However, element deletion induces volume loss in the mesh and an approximative cutting surface. Other methods use element duplication for each element crossing the cutting path Molino et al. (2005), but the volume of the mesh is still modified. Some methods were proposed to alleviate these issues by re-meshing locally to fit the cutting path [START_REF] Paulus | [END_REF]. A significant difficulty is to preserve the quality and the density of the mesh during the subdivision process to avoid distorted elements, which lead to convergence difficulties during the simulation. Several approaches were proposed to maintain a relatively good mesh quality Ganovelli and OSullivan (2001); Bielser et al. (2003); Molino et al. (2005); Sifakis et al. (2007).

Meshless methods use particles and their interactions with neighbors to compute the deformation. In [START_REF] Cheng | An interactive meshless cutting model for nonlinear viscoelastic soft tissue in surgical simulators[END_REF], the cutting process is simulated when the forces increase above a threshold separating the two parts around the cutting plane. The main advantage is the reduced computation time, but meshless methods lack geometrical information that leads to less efficient collision detection.

Hybrid methods combine the advantages of both using a meshed surface and particles to compute the deformations. In Si et al. (2018), the authors reconstruct the surface mesh using an interactive triangulation algorithm during the cut. Alternatives to generating finite element meshes, in particular, to handle topological changes, are extended finite element methods XFEM Moës et al. (1999); [START_REF] Koschier | [END_REF], that will enrich the shape functions of the elements to introduce discontinuities. These methods allow cuts, material interfaces, and domain boundaries to be described independently of a background mesh, which may also be progressively adapted using a posteriori error estimators Bordas and Duflot (2007); Duflot and Bordas (2008); Bordas et al. (2008); Ródenas et al. (2008) or local heuristics Menk and Bordas (2011). However, such methods are not developed in the real-time context, and the simulation of collisions and interactions remains unsolved, mainly because the geometry of the discontinuities is known implicitly.

Haptic rendering

Haptic rendering enables physical interaction with simulated objects in a virtual environment. Haptic devices control an object in the virtual envi-ronment and create an experience of touch by applying forces, vibrations, or motions to the user. The manipulated objects interact with each other generating reaction forces affecting the physics-based models, and more importantly, the controller of the haptic interface. Somehow, it represents a way to evaluate the level of accuracy of physics-based models in an interactive simulation.

Haptic simulations raise additional constraints. The computation of the reaction forces is limited by 1) the need for high-frequency refresh rates (between 500Hz -1KHz), 2) the stability of the device control law, and 3) the fidelity of the haptic rendering as perceived by the user. Several approaches have been proposed to deal with these constraints (see [START_REF] Susa | Haptic rendering based on finite element simulation of vibration[END_REF] for a broad survey).

In the pioneering works Bro-Nielsen and Cotin (1996); Cotin et al. (1999), positions are applied as bilateral constraints (equality conditions) and solved by the Lagrange Multipliers method. The approach was initially based on a superposition principle and was further extended Picinbono and Lombardo (1999) using a force extrapolation method. Later, Duriez and Andriot (2003) introduced Signorini's model for contact handling in the field of haptic rendering. In [START_REF] Duriez | [END_REF], an extension to friction contact response is proposed. The computation performance is obtained through the use of a pre-computed and condensed compliance matrix. However, models are limited to small displacement deformations, which is not realistic for soft-tissue models.

An intermediate representation of the constraints provided by the virtual environment was proposed in Adachi et al. (1995), Mark et al. (1996) or [START_REF] Garre | [END_REF]. This representation separates the haptic rendering from the physics engine, usually too expensive for haptic rate computation. In Jacobs et al. (2010) the simplified model relies on the linearization of nonlinear deformable objects that runs at low rates. In Mahvash and Hayward (2004) and [START_REF] Sedef | [END_REF] stable haptic rendering is efficiently achieved by calculating response forces from pre-computed data, e.g., by interpolations performed directly in the haptic loop. Saupin et al. (2008a) uses an asynchronous strategy between simulation and haptic rendering, sharing the full compliance matrix in the contact space. However, a precomputation of the compliance is used, preventing any use of hyperelastic models. The method was extended in [START_REF] Peterlík | Real-time visio-haptic interaction with static soft tissue models having geometric and material nonlinearity[END_REF] for needle insertion simulations, where a beam model computed at haptic rate is coupled with the pre-computed compliance matrix.

Methodology

Our methodology is based on the following general background of deformable simulations with contact using an implicit integration time.

FE models and constitutive law

We first introduce the deformation models and constitutive laws. where D e corresponds to the stress-strain matrix parametrized by the Young's modulus E and the Poisson's ratio ν:

D e = E (1 + ν)(1 -2ν)         (1 -ν) ν ν 0 0 0 ν (1 -ν) ν 0 0 0 ν ν (1 -ν) 0 0 0 0 0 0 1-2ν 2 0 0 0 0 0 0 1-2ν 2 0 0 0 0 0 0 1-2ν 2         (2.
2) C e is the strain-displacement matrix which is employed to compute strains at any point inside the element, based on its nodal displacements:

C e = 1 6V         b 1 0 0 b 2 0 0 b 3 0 0 b 4 0 0 0 c 1 0 0 c 2 0 0 c 3 0 0 c 4 0 0 0 d 1 0 0 d 2 0 0 d 3 0 0 d 4 c 1 b 1 0 c 2 b 2 0 c 3 b 3 0 c 4 b 4 0 d 1 0 b 1 d 2 0 b 2 d 3 0 b 3 d 4 0 b 4 0 d 1 c 1 0 d 2 c 2 0 d 3 c 4 0 d 4 c 4         (2.3)
where V is the volume of the tetrahedral element. The method of cofactors can be used to invert the following 4 x 4 matrix and obtain the above coefficients:

    1 x 1 y 1 z 1 1 x 2 y 2 z 2 1 x 3 y 3 z 3 1 x 4 y 4 z 4     -1 =     a 1 a 2 a 3 a 4 b 1 b 2 b 3 b 4 c 1 c 2 c 3 c 4 d 1 d 2 d 3 d 4     (2.4)
where x i , y i , z i are the nodal coordinates of the tetrahedra in the global coordinate system.

All the above matrices are constant and can be pre-computed since they depend only on the position of the mesh at the initial state. However, when large rotations occur, the linear approximation of the strain function leads to well-known inflation artifacts Müller and Gross (2004).

Co-rotational formulation

The corotational formulation is widely used in computer graphics (Felippa, 2000) for its stability and the fast computation time possible thanks to precomputations, while the model is not restricted to small displacements. For this purpose, the popular corotational method computes displacements in a rotated local coordinate system (Figure 2.2). In each element, the local frame is represented by a rigid rotation matrix R e , and its transpose is used to align the deformed tetrahedron with its reference, undeformed shape. The forces f e are computed in the rotated coordinate system, then transformed back to the world coordinates:

Re

f e = R e K e R T
e u e (2.5) where u e = pp is the displacement vector between p and p, being respectively the current and initial positions.

The rotation matrix R e associated to the element e is a 12×12 block-diagonal matrix computed as follows:

R e =     Re Re Re Re     (2.6)
where Re is a 3×3 rotation matrix. Several solutions have been proposed to compute the rotation matrix. For instance, it can be obtained from a polar decomposition of the nodal positions as proposed in Nesme et al. (2005).

MJED method for Hyperelastic materials

Any hyperelastic material is entirely determined by its strain energy function, which describes the amount of energy necessary to deform. This strain energy function often involves the invariant of the right Cauchy-deformation tensor defined as C = ∇Φ T ∇Φ, where Φ is the deformation function from the configuration at rest towards the deformed one. Most of the strain energy of hyperelastic isotropic materials is defined from 2 invariants of the deformation tensor I 1 = trC, I 2 = ((trC) 2 -trC 2 )/2 and the Jacobian J = det∇Φ.

The MJED method Marchesseau et al. (2010) proposes to decouple the invariant of C from the expression of J to avoid assembling the stiffness matrix and complex derivative expressions. The rationale behind the MJED method is that the derivative of J with respect to the nodal position is trivial for most elements (especially for linear tetrahedra), whereas its derivative with respect to C is not trivial. Since the converse is valid for the two invariants, the MJED method optimizes the force and the stiffness computation by separating both terms.

Implicit time integration

Equations used to model the dynamic behavior of bodies can be written within a synthetic formulation, given by Newton's second law:

M(p) a = G(t) -F (p, v) (2.7)
where p ∈ R n is the vector of generalized degrees of freedom (here, the mesh node positions), M(p) is the inertia function from which derives the Mass matrix M (assumed constant in this document). The function G gathers external forces at time t. F is a nonlinear function providing internal forces of the simulated object depending on its current position p. The internal forces are derived from the physics-based deformable model described above (i.e., corotational formulation).

Let's consider the time interval [t i , t f ] whose length is h = t f -t i . We have:

M(v f -v i ) = t f t i (G(t) -F(p, v)) dt p f = p i + t f t i v dt (2.8)
To evaluate integrals

t f t i (G(t) -F(p, v, t)) dt and t f
t i v dt we chose the following implicit Euler integration scheme:

M(v f -v i ) = h (g f -F(p f , v f )) p f = p i + hv f (2.9)
where g f is the value of function G at time t f .

We apply a Taylor series expansion and make a first-order approximation of the nonlinear function F:

F (p i + ∆p, v i + ∆v) ≃ f i + δF δp ∆p + δF δv ∆v (2.10)
This linearization corresponds to the first iteration of the Newton-Raphson algorithm. Therefore, the problem is solved with a single iteration, assuming a temporal coherency of the mechanical behavior. It may lead to minor numerical errors in the dynamic behavior, but these errors tend to decrease at equilibrium or with null velocity.

Replacing (2.10) in (2.9) and using ∆p = p f -p i = hv f and ∆v = v f -v i , we obtain:

(M + hB + h 2 K)∆v = h(g f -f i ) -h 2 Kv i (2.11)
where M is the mass matrix, f i and v i are respectively the internal forces and velocities at the beginning of the time step. K = δF δp is the global stiffness matrix and B = δF δv is the damping matrix, used to dissipate energy. Following [START_REF] Baraff | [END_REF], the damping matrix can be approximated with Rayleigh damping B ≃ αM + βK, where α is the rayleigh stiffness and β the rayleigh damping. Finally, equation (2.11) can be rewritten:

((1 + hα)M + h(h + β)K) A ∆v x = h(g f -f i ) -h 2 Kv i b (2.12)
Equation ( 2.12) provides a linear problem Ax = b that must be solved at each time step. Indeed, the left-hand side is a global system matrix A mainly composed of the mass and stiffness matrices. Using the co-rotational formulation the stiffness matrix A is obtained with:

K = (G e R e K e R T e G e T ) (2.13)
where G e is the globalization matrix transferring the local stiffness matrix K e of an element e to the global stiffness matrix K. The rotation matrices are not constant and must be recomputed at each simulation step. Similarly, although the hyperelastic MJED formulation does not rely on rotation matrices, the local stiffness matrices K e are not constant and depend on the current position/deformation of the model. Therefore, in both cases, K depends on the position of FE models. As a result, the linear system must be solved at each step. In an interactive context, the equation (2.12) is usually solved with the iterative Conjugate Gradient.

Preconditioner

The condition number κ of the matrix A measures how much the output value of the function can change for a slight change in the input argument.

For heterogeneous objects or ill-structured meshes, the condition number κ is often high, which raises convergence issues for the conjugate gradient algorithm. A common technique is to use a preconditioner to reduce the condition number, ensuring a faster convergence of the algorithm. By definition, a preconditioner approximates the system matrix A, which is less costly to invert. Solving equation (2.12) with a preconditioner P can be written:

P -1 Ax = P -1 b, such that κ(P -1 A) < κ(A) (2.14)
Several preconditioners can be used, from simple diagonal matrices [START_REF] Baraff | [END_REF] to accurate but costly Cholesky factorization. However, the performance improvement remains limited since the preconditioner itself is expensive to compute.

Time-stepping and collision detection

Collision response on mechanical objects leads to discontinuities in the velocities of the colliding points. For such discontinuous events, the acceleration is not defined: the problem belongs to the field of non-smooth mechanics.

To integrate the mechanics and the non-smooth events due to contact over time, we use a time-stepping method: The time step is fixed, and there is no limitation on the number of discontinuities that could take place during a time step (M.∼Anitescu et al. (1999)), but low-order integration schemes should be used. Although it could lead to excessive dissipation if the time step is too significant, it provides stable simulations. It is particularly relevant for interactive simulations involving contact with virtual devices controlled by an operator.

Collision detection

Before enforcing the contact between soft tissues or with surgical instruments (rigid or deformable), one needs to detect them (see Teschner et al. (2005) for a survey). A common strategy is to use proximity queries providing the minimal distances between mesh (even concave meshes). This approach has the advantage of "anticipating" the contacts before they appear. Specific structures such as hierarchical bounding volume may be used to accelerate the computation time of this task, such as bounding boxes Gottschalk et al. (1996) or spheres Hubbard (1995).

Collision detection returns either penetration depth or a minimal distance between meshes or pairs of primitives (between which distance constraints can be formulated to solve the contacts).

Contact mapping

Collision detection is a time-consuming operation. As introduced in [START_REF] Faure | [END_REF], to save computation time, the collision surface of the models is not necessarily the same as the underlying FE mesh used to compute the deformation. In this case, the constraint forces can be transferred from the collision models to the FE model using a mapping function J . This relation is used to map the positions p fe of the deformable model to the positions p col of the collision models:

p col = J (p fe ) (2.15)
The most used mapping is to couple the collision surface (edges, triangles, or quads) with the tetrahedral mesh supporting the FE computations. In this case, a barycentric mapping is used to link the different representations. At the beginning of the simulation, each vertex of the collision models is associated with the "closest" tetrahedron in the FE mesh. During the simulation, both the collision mesh and the FE model are deformed simultaneously. Therefore, this relation remains constant over the simulation:

p col = J col p fe (2.16)
where J col is the matrix of the mapping computed from the barycentric coordinates of collision points with respect to the associated tetrahedron's positions. In addition, the mapping also provides a way to transfer any given force f col applied on the collision models to an equivalent force applied on the FE nodes. The transferred forces f fe are then weighted with the same barycentric coefficients using:

f fe = J T col f col (2.17)
Similarly, collision detection may provide contact information at an arbitrary location on the collision surface (for instance, in the middle of a collision triangle). Another mapping relation is, therefore, necessary to transfer contact forces to equivalent nodal forces on the collision mesh (Figure 2.3). To simplify the solving process this relation is considered constant during the resolution of each time step. As previously, a mapping matrix J CD can be defined to transfer the contact forces f prox on the collision mesh: Finally, for a given object an arbitrary number of representation may be used and the mapping relations may be serially applied to finally compute equivalent forces (or displacements) on the DOFs of the FE mesh f fe = J T col . . . J T CD f prox . This operation is fast since mapping matrices are usually very sparse, and most of them can be pre-computed. Therefore, for simplicity, we now omit this relation considering that the constraints are directly applied to the mechanical model.

f col = J T CD f prox (2.18)

Constraint-based simulation

To avoid interpenetrations or to enforce specific behaviors, the simulation of two colliding objects 1 and 2 is now subject to a set of constraints. For this purpose, the simulation of the deformable bodies are coupled with the following equation:

Λ(p f,1 , p f,2 ) = H 1 (p f,1 ) -H 2 (p f,2 ) (2.19)
where Λ is a nonlinear function providing the violation of the constraint with respect to positions p f,1 and p f,2 of the deformable objects at the end of the time steps, and H 1 , H 2 provide the value of the constraints. In the case of collisions, Λ corresponds to the shape of the interpenetration (red area in figure 2.4) and must be canceled at the end of each time step. This function is highly nonlinear since constraint forces must be applied to cancel the violation. However, applying such forces will result in displacements, rotations, or deformations of deformable bodies that may generate new collision areas or invalidate the previous ones.

The constraint equations are linearized using a first-order Taylor expansion:

H(p f ) = H(p i + hv f ) ≈ H(p i ) + h ∂H(p i ) ∂p v i + h ∂H(p i ) ∂p ∆v (2.20)
where dH(p i ) dp is the linearization of constraints at the begining of the time step (using positions p i ). In practice, the evaluation and the linearization of the constraints equations are complex. For this purpose, collision detection is performed, providing a discrete set of proximities between both objects (green lines in figure 2.4). The number of proximities usually depends on the resolution of the collision mesh (that may be adapted thanks to the mappings) and the collision detection method itself (for instance, by filtering proximities afterward).

Once the collision information is available, all the constraint equations are evaluated along with the directions of the proximities assumed constant for the time step. This leads to several simplifications: First of all, H ≈ dH(p i ) dp , known as the Jacobian of the constraints, can be defined2 , providing the constraint directions (blue and orange arraows in figure 2.4). The dimension of H is c × n where c is the number of proximities, and n is the number of DOFs. Similariliy, the violation of the contraints H(p i ) ≈ Hp i is evaluated along the same constraint directions. With these simplification and replacing (2.20) in (2.19) the violation of the constraint at the end of the step can be rewritten as:

δ f = (H 1 p i,1 + hH 1 v i,1 ) + (H 2 p i,2 + hH 2 v i,2 ) δ i +h H 1 ∆v 1 +h H 2 ∆v 2 (2.21)
where δ f is the violation of the constraints at the end of the time step and

H 1 = dH 1 (p i,1 ) dp and H 2 = - dH 2 (p i,2 ) dp .
In addition, to enforce non-penetrating objects, constraint forces λ (also called Lagrangian Multipliers) are applied along the direction of the constraint. The governing equation (2.7) is modified as follows:

M(p) a = G(t) -F (p, v) + hH T λ (2.22)
After the linearization, equation (2.12) is modified as follows:

Ax = b + hH T λ (2.23)
Gathering equations (2.23) and (2.21) for the two interacting objects results in the following Karush-Kuhn-Tucker (KKT) system:

     A 1 x 1 + hH T 1 λ = b 1 A 2 x 2 + hH T 2 λ = b 2 hH 1 x 1 + hH 2 x 2 = ∆δ (2.24)
(2.25) (2.26) where ∆δ = δ f -δ i 2.3.4.1 Contact and friction models Contact equations are subject to Signorini's law (Figure 2.5a). It expresses that there is a complementarity3 relation between the contact force λ n and the distance δ n , along the direction n of the contact:

0 ≤ δ n ⊥ λ n ≥ 0 (2.27)
This model has several physical justifications including non-interpenetration and no sticking force. Moreover the contact force vanishes if the points are not strictly in contact. However, using Signorini's law, the contact force creates a frictionless response. Coulomb's friction law describes the macroscopic behavior of friction in the tangent contact space (Figure 2.5c). In this law, the reaction force is included in a cone whose height and direction are given by the normal force. If the reaction force is strictly included inside the cone, objects stick together;

otherwise, the reaction force is on the cone's border, and objects are slipping along the tangential direction. In this last case, the friction force must be directed along the direction of motion:

δT = 0 ⇒ ∥λ T ∥ < µ ∥f n ∥ (stick) δT ̸ = 0 ⇒ λ T = -µ ∥λ n ∥ δT ∥ δT ∥ = -µ ∥λ n ∥ T (slip) (2.28)
where µ is the friction parameter, and T is the direction of motion in the tangential plane to the contact normal n.

During 3D slipping motion (also called dynamic friction), the tangential direction is unknown. Indeed, the tangential velocity depends on the tangential force that has to be found during the solving process. It creates a non-linearity equation in addition to the complementarity state stick/slip.

Signorini's law and Coulomb's law are also valid in a multi-contact case. This model is defined for each potential contact provided by the contact detection algorithm. However, to solve these laws at every contact point, we have to consider the coupling that exists between these contact points (Figure 2.5b). This coupling comes from the intrinsic mechanical behavior of deformable objects.

Constraint solving and mechanical coupling

When objects are distant, the mechanical matrices A 1 and A 2 are not coupled and can be solved independently (Figure 2.6a). However, coupling terms must be added to take the constraints into account.

A 2 A 1 A 1 A 2 x 1 b 1 x 2 b 2 (a) No contact A 2 A 1 A 1 A 2 x 1 b 1 x 2 b 2 (b) Penalty A 2 A 1 A 1 A 2 H 1 T H 1 x 1 b 1 x 2 b 2 𝜆 𝛿 H 2 T H 2 (c) Augmented Lagrangian
Figure 2.6: When contacts occur, penalty-based methods 2.6b add off-diagonal terms coupling the mechanical matrices. In this case, the global coupled system must be solved at once. Instead the Augmented Lagrangian method 2.6c adds additional equations to couple the models but mechanical matrices are not modified.

Compared to penalty-based methods, a significant advantage of relying on Lagrange multipliers is that collision events never modify the mechanical matrices. It is therefore possible to rely on the schür complement method to solve the equation (2.24), (2.25) and (2.26) with the 5 following steps:

Step 1: interacting objects are solved independently while setting λ = 0. A set of independent linear systems of equations Ax = b must then be solved for each object. The Conjugate Gradient algorithm is generally proposed to solve this problem since A is large, sparse, symmetric, and positive definite. It provides

x free 1 = A 1 -1 b 1 and x free 2 = A 2
-1 b 2 called the free motion for each object, corresponding to the velocities update without considering any constraint nor collision. After integration, we obtain p free 1 and p free 2 :

p free 1 = p i,1 + hv i,1 + hx free 1 (2.29) p free 2 = p i,2 + hv i,2 + hx free 2 (2.30)
Step 2: the constraint laws are linearized during the time step. The linearization corresponds to the collision detection that provides the proximities of potential contacts and their associated normals between the position p1 and p2 of the collision meshes. Constraint equations are then transferred to the mechanical DOFs using the mapping relation (2.17), providing the Jacobian of the constraint H 1 and H 2 that are assumed constant during the time step (see [START_REF] Faure | [END_REF] for details). The collision detection is performed at the beginning of the time step, whereas the violation δ free is instead computed using the free motion p free 1 and p free 2 obtained at the previous step.

Step 3: The violation of the constraint is evaluated, and constraints are defined according to the collision proximities. With Euler equation one may note that:

∆δ -H 1 x free 1 -H 2 x free 2 = δ f -H 1 p free 1 -H 2 p free 2 (2.31)
Replacing (2.24) and (2.25) in (2.26) gives:

H 1 A -1 1 H T 1 + H 2 A -1 2 H T 2 W λ = H 1 p free 1 -H 2 p free 2 δ f -δ free (2.32)
where W is known as the compliance matrix, or the Delasus operator.

In the above equation, both λ and δ are unknown. For Signorini's law (equation (2.27)), this equation describes an LCP (Linear Complementarity Problem). If it is combined with Coulomb's law (equation (2.28)), we obtain an NLCP (Nonlinear Complementarity Problem). An important difficulty lies in the fat that this equation involves evaluating the inverse of large matrices A 1 and A 2 (same dimension as the number of DOFs) to compute the compliance matrix W.

Step 4: We obtain the value of λ using a Gauss-Seidel algorithm dedicated to the NLCP created by contact and friction equations. Considering a contact α, among m instantaneous contacts, one can write the behavior of the model in contact space:

δ α -W αα λ α unknown = α-1 β=1 W αβ λ β + m β=α+1 W αβ λ β frozen + δ free α (2.33)
where W α,β is a compliance matrix modeling the mechanical coupling between contact points α and β. For each contact α, this method solves the contact equations by considering the other contact points (α ̸ = β) as "frozen".

The new value of λ α is given by solving Signorini's law and Coulomb's law on this contact (see [START_REF] Duriez | [END_REF] for details of implementation and performance).

Step 5: When the value of λ is available for all the contacts, the corrective motion is computed:

p 1,t+h = p free 1 + h∆v cor 1 with ∆v cor 1 = A -1 1 H T 1 λ p 2,t+h = p free 2 + h∆v cor 2 with ∆v cor 2 = A -1 2 H T 2 λ
(2.34)

We finally obtain p 1,t+h and p 2,t+h , the positions of objects 1 and 2 that fulfills the contact and friction laws.

Compliance and mechanical coupling

Equation (2.32) requires the computation of A -1 , which is a large matrix (same dimension as the number of DOF) and changes at each time step. Although computing this inverse in real-time is only possible for coarse models, the resulting operator W plays a critical role in enforcing the constraints.

In the following, the term "mechanical coupling" describes the coupling between contact constraints applied on two subsets of the boundary of a deformable body. This coupling occurs through the deformation of the body itself. Indeed, even if the contact points are only defined on the boundary of the deformable bodies, they are all influenced by each other through the stiffness of the material.

Consider for example the figure 2.7. Although none of the deformed shapes show any interpenetration, the behavior computed in the figure 2.7b is not acceptable since stiff and soft parts are deformed the same way. It shows that the contact force distribution is closely related to the underlying heterogeneity of the material, which is represented by W. With an explicit integration scheme, the matrix W corresponds to the diagonal mass matrix. Similarly, if penalty methods are used, the force distribution would mainly depend on the geometrical interpenetration but not on the heterogeneity. In both cases, it would lead to unrealistic configurations such as in figure 2.7b, at least during the dynamic of the deformation. Saupin et al. (2008b) proposed the compliance warping technique, which consists in updating the pre-computed A 0 -1 with the nodal rotations, but this solution remains inaccurate for large deformations and requires storing a large dense matrix which makes the method unsuitable for fine meshes.

Contributions

In the previous section, we identified two critical difficulties to achieve realtime computations. (i) to compute the solution of the linear system defined in equation (2.12), in particular for heterogeneous objects. (ii) to compute the compliance matrix W in equation (2.32). These two problems are addressed in this section.

2.4.1 Fast linear solver for implicit integration 2.4.1.1 Implicit FEM Solver on GPU for Interactive Deformation Simulation In [START_REF] Allard | [END_REF], we introduced a set of methods to implement an implicit Finite Element solver on the GPU. Unlike previous GPU-based sparse solvers, we avoid the assembly of the system matrix and parallelize the matrix operations directly on the original object mesh. It considerably reduces the number of operations required, and more importantly, the consumed bandwidth, enabling the method to be fast enough for highly complex interactive stiff body simulations.

While an explicit integration scheme can be implemented with very few tasks (figure 2.8b), a lot more are required for implicit integration (figure 2.8c).

Yet, most tasks are standard linear vector algebra that can easily be implemented either on the CPU or GPU, except for Force and dForce (force evaluation and stiffness matrix / displacement vector product, i.e. lines 1 and 8 in the figure 2.9). These two tasks are the most expensive. The task Force is executed only once per simulation step, but the dForce is required for each iteration of the CG. We proposed a GPU-based implementation of these two tasks resulting in a complete GPU parallelization of the CG, without any need for transfers between CPU/GPU except to check convergence at each iteration (line 16 in the figure 2.9). For this purpose, we avoid constructing the large system matrix explicitly; instead, we compute matrix-vector products for each element in the original unstructured finite element mesh. The implementation of this step is based on the following techniques: 13 δ i = dot(r,r); Parallelization strategy: While the actual computations depend on the specific FEM formulation, the parallelization of tasks Force and dForce only depends on the mesh. Several elements are connected to each vertex, causing writes conflicts if we process them in parallel. Classically, FEM solvers rely on mesh partitioning, but this method cannot extract enough parallelism for GPUs. Instead, we rely on a two-step process, where values from each element are first computed independently, then each vertex gathers contributions from connected elements.

1 f 0 = f(p,v); // Force 2 b = f0 + h K v; 3 a = 0; 4 // Solve (M-h 2 K) a = b 5 d = r = b; 6 δ 0 = dot(r,r); 7 for( i=1; i ≤ imax; ++i ) { 8 df = K d; // dForce 9 q = M d -h 2 df; // mass 10 α = δ i-1 / dot(d,q);
14 β = δ i / δ i-1 ; 15 d = r + β d; 16 if ( δ i > ϵ 2 δ 0 ) break; 17 } 18 v = v + h a; // use solution 19 p = p + h v;
Blocked data layout: Necessary attention must be paid to the memory alignment to reach efficient performances on the GPU. A common approach to improve coalescing is to convert large arrays-of-structures into structuresof-arrays. A more efficient approach on a GPU was proposed using an array of structures where each member is replaced by a b-sized array, where b is the number of threads running on the same multiprocessor.

Matrix-less scheme: In a co-rotational FEM formulation, each tetrahedron contributes a dense 12×12 matrix, resulting from a product of 5 smaller matrices, some of which change at each step. Instead of computing and storing the final matrix, we store the few scalars required to express the small matrices, and each matrix-vector product is decomposed into five successive but more straightforward products.

Reducing launch overheads: Only a few milliseconds are available in an interactive application to compute one full-time step, involving hundreds of different computation steps. In such a context, the achieved performance can be significantly affected by the launch overheads of each kernel. Therefore, we proposed to merge as many steps as possible into a single kernel. The final CUDA-based implementation can simulate a deformable object with 45k tetrahedral elements at 212 FPS on an NVIDIA GeForce GTX 480, 18× faster than our most optimized sequential implementation on an Intel Core i7 975 3.33GHz CPU. We ran benchmarks to evaluate the speedup on several GPUs with varying mesh sizes (figure 2.10). With the largest meshes, the speedup can reach 23× on a GeForce GTX 480 and 13.5× on the previous generation GTX 280.

Asynchronous Preconditioners for Efficient Solving of Non-linear Deformations

In Courtecuisse et al. (2010a) we proposed an asynchronous preconditioning technique. The method relies on the assumption that A in equation (2.12) undergoes small perturbations between two consecutive time steps. Indeed, if

P t = A -1
t is available at a specific time t, it remains a "good" approximation for the following time steps. The preconditioner can then be updated at a low frequency on a dedicated CPU thread, and the last preconditioner available can be used to advance the simulation (Figure 2.11). Therefore, the overhead in computing the preconditioner is removed from the simulation loop, which allows using more advanced and computationally costly preconditioners such as a complete factorization of the system4 :

P = A = LDL T (2.35)
where D is a diagonal matrix and L is a sparse lower-triangular matrix.

Instead of Cholesky factorization, we choose LDL T factorization since it produces more numerically stable results than a Cholesky factorization. The factorization is performed by the cs sparse library [START_REF] Davis | [END_REF], using a single core on the CPU. Other libraries [START_REF] Toldeo | [END_REF]Schenk et al. (2008) propose parallel factorization, but we found that cs sparse provides sufficiently fast updates, and a sequential factorization enables to save CPU cores from computing the preconditioners of other objects in parallel. The preconditioner is updated asynchronously within a dedicated CPU thread. We use the last preconditioner available to advance the simulation so that the simulation never needs to wait for the computation of the current preconditioner to be complete.

An essential advantage of this factorization is that the resulting L matrix remains sparse, making the preconditioner application faster within the CG. This operation consists in solving two Sparse Triangular Systems (STS):

y = L T -1 b x = L -1 (D -1 y) (2.36)
where L is stored in Compressed Row Storage (CRS) (see S. et al. (1995) for details on the structure). Solving the STS is equivalent to performing a Gauß elimination, which is difficult to parallelize as it involves many dependencies. Therefore, the STS is solved on the CPU, taking advantage of caches and sparsity to solve the system efficiently. Combining the GPU-based CG with a CPU-based preconditioner implies transferring the solution vector b between the CPU and the GPU at each iteration of the CG. However, since the preconditioner is usually a good approximation of the actual system, only a few iterations are necessary, and the cost of such transfers remains limited. For large systems, the computation cost of the LDL T factorization can become prohibitive, and the resulting preconditioner can diverge from the actual system. However, we note that an essential part of the error is associated with the rotations Saupin et al. (2008b) which quickly vary between time steps. In order to limit the divergence of the preconditioner, we estimate the nodal rotations R t-h→t that were introduced since the last update of the preconditioner (i.e. between time t -h and t). Indeed, the nodal rotations approximate the co-rotational formulation where rotations are computed per element and accumulated in the global stiffness matrix. R t-h→t is a block diagonal matrix, easy to compute and easy to invert. This way, the most recent preconditioner P t-h is then rotated with the current rotation matrix R t-h→t as follows: 2.37) where the "rotated preconditioner" P r is less sensitive to geometrical nonlinearity. Finally, the method allows simulating the deformation of homogeneous as well as heterogeneous tissues in real-time.

P r = R T t-h→t L t-h D t-h L T t-h R t-h→t (
Convergence rate and computation time The method was compared with several preconditioners and update strategies. In the figure 2.13, all preconditioners use an Incomplete Cholesky factorization, but with different update strategies: CG and PCG-1 correspond to the classical CG algorithm with and without preconditioning. PCG-1 always provides an updated preconditioner resulting in the minimum iterations necessary to converge, but it also provides the worst performance since it requires computing an incomplete factorization of the matrix at each time step. PCG-0 corresponds to a strategy where the preconditioner is pre-computed and is never updated.

While it provides a reasonable convergence rate, it quickly diverges when large deformations are introduced. PCG-30 corresponds to a preconditioner updated every 30 simulation steps. The convergence rate is equivalent to PCG-1 at the time where the preconditioner is updated, but it introduces spikes in the computation time to compute the factorization. When the preconditioner is computed asynchronously (PCG-30 async), the factorization does not affect the computation time of the simulation. However, between updates, the quality of the preconditioner decreases. This drawback is alleviated in the final scheme (PCG-30 async+warp), which uses local rotations between updates. This scheme obtains a nearly optimal preconditioner at each step of the simulation, with both few iterations and fast computations. In figure 2.14 we show the importance of the method for heterogeneous structure. We produced a simulation where a heterogeneous beam composed of 6, 500 elements and 1, 470 nodes is deformed under gravity. Direct solvers must process a complete factorization of the matrix at each step to take into account the non-linearity of the model. This operation is expensive, but it can be parallelized using optimized libraries such as Pardiso Schenk et al. (2008). However, using 4 threads (Pardiso-4) is only 2.2× faster than the sequential version, and only 1.90× with 8 threads (Pardiso-8). Standard preconditioners (standard) require a large number of iterations due to the significant heterogeneities. Pre-computing the LDL T factorization (no update) and using it throughout the simulation enables removing the overhead of the factorization while keeping a limited number of iterations. However, in cases of large deformations, the actual stiffness of the material may be very different from the rest configuration, and a large number of iterations are necessary. Using our method, the factorization time of the preconditioner is still negligible, and the number of iterations remains very low throughout the simulation. For this preconditioner, the (warping method ) provides only limited improvements. Indeed, the preconditioner was updated on average every 8 time step, and the rotations between two consecutive updates remained limited. Inverse corresponds to the inversion of the diagonal matrix for the Jacobi preconditioner and the factorization of the system for LDL T preconditioners. Solving is the time taken to solve the system. Total is the total time of a single time step.

Real time Simulation of contact 2.4.2.1 GPU-Based preconditioner for contact problems

In Courtecuisse et al. (2011a) we proposed to rely on the asynchronous preconditioner of the previous section as an approximation of the compliance matrix W:

W = HA -1 H T ≈ HP -1 H T = H R LDL T R T -1 H T (2.38)
Indeed, since P represents a close approximation of the factorization of A, we propose to use it to compute W in equation (2.32). The above equation requires computing the product of the inverse of the preconditioner with the Jacobian of contacts H, which can be achieved by computing columns independently of the following matrix:

LDL T X i = H i ⇔ X i = LDL T -1 H i (2.39)
where X gives the result of the inverse of the preconditioner times the Jacobian of the contact. This operation is still expensive since multiple triangular systems must be solved for each right-hand side vector. However, the multiple right-hand side vectors stored in H can be computed independently. Therefore, we assign the computation of each column to an independent multiprocessor on the GPU. Each group is therefore fully processed by a single processor (Figure 2.15), which enables fast local synchronizations directly on the GPU. Then we use another level of parallelism where each STS is solved with several threads. Indeed, many data can potentially be treated in parallel during the solving process of each STS. This two-level strategy fits the GPU architectures where local synchronizations within a group of threads are fast, whereas global synchronizations over multiple groups are costly. In the figure 2.16, we measure the time needed to compute the compliance matrix according to the number of constraints. Solving the different STS on the GPU (LDL GPU) is much faster than solving each STS on the CPU (LDL CPU). Indeed, the computation time for the CPU version is linear according to the number of right-hand side vectors, whereas the GPU can process them in parallel. Therefore, for up to 78 constraints, the computation time remains almost constant with our GPU implementation. Indeed, below this limit, the GPU computing units are not fully utilized, and the different STS are processed in parallel. Beyond this number, some GPU processors will compute several STS successively, and the computational time curve takes a staircase appearance. Nevertheless, the GPU processors can overlay waiting times, due to synchronizations and access in memory, with computations for another STS. Thus, solving the system for 140 constraints is only 1.8 times slower than for 70 constraints.

Number of constraints

Virtual Cutting of Deformable Objects based on Efficient Topological Operations

The ability to simulate surgical cuts, dissection, soft tissue tearing or microfractures, is essential for augmenting the capabilities of existing or future simulation systems. A new remeshing algorithm with a fast FE approach was proposed in [START_REF] Paulus | [END_REF] to support such features. The method is generic enough to support a large variety of applications. The remeshing algorithm we present is an extension of the √ 3-subdivision scheme introduced by Burkhart et al. (2010) (Figure 2.17). This approach soundly and efficiently handles partial cuts, the emergence and propagation of cracks, and the existence of overlapping or crossing between separation surfaces. We show the benefits of our approach, evaluating the impact of cuts on the number of nodes and the numerical quality of the mesh. These points are crucial to ensure both accurate and stable real-time simulations. Beyond remeshing, the second issue is adequately updating the deformable model's mechanical properties and equation systems when the cut is introduced. Indeed, the preconditioner introduced above is updated with delay, and the modifications significantly affect its efficiency. Moreover, the delay of the updates can lead to instabilities and inaccuracies, particularly when treating contact with the instruments and self-collisions between different parts of the cut. In Courtecuisse et al. (2014a), we proposed to combine the method with the Sherman Morrison Formula (SMF) to compute the correction of the preconditioner due to the topological changes (Figure 2.18). An essential advantage of updating the preconditioner is that it helps to maintain the number of affected nodes by the perturbations at a minimum. Indeed, each new factorization implicitly contains all anterior modifications to the last update of the preconditioner. The SMF correction is, therefore, necessary only for the perturbation that appeared since the last update, which only involves a few affected nodes (Figure 2.19). When a modification is performed on the mesh, we first compute the correction of the current factorization. Then we compute the correction of the preconditioner, which was being calculated at the time of the cut. After two consecutive updates without topological modification, the preconditioner does not need any additional correction.

Haptic Rendering of Hyperelastic Models with Friction

In Courtecuisse et al. (2015) we introduced a method for the haptic rendering of hyperelastic materials. Such simulations are known to be difficult due to the non-linear behavior of hyperelastic bodies. Soft tissue dynamics are simulated in real-time (20 to 100 Hz) using the Multiplicative Jacobian Energy Decomposition (MJED) method. The compliance matrix W is asynchronously evaluated at low rates, while a local contact problem is computed at high rates to update haptic forces. To update haptic forces, our approach shares constraint equations and compliance between a haptic control loop (at high rates ≃ 1 kHz) and the simulation (at low rates ≃ 30 Hz). After performing the collision detection and the computation of deformable models, the NLCP is defined and solved in the simulation. The compliance matrix W and the Jacobian of the contacts H are shared with a separate haptic loop. In the haptic loop, the position of the device is refreshed to update the motion of the virtual instrument driven by the device motions (Figure 2.20). The displacement provides a new violation in the constraint space, then a new value of the contact response is computed and sent to the force feedback device.

Applications

To highlight the generic nature of the above contributions, they are used in various medical applications, i.e., a cataract surgery simulation, a hepatectomy in laparoscopic surgery, and cerebral tumor removal.

Application to cataract surgery: The cataract is an opacification of the eye's lens, which prevents the light from reaching the retina, resulting in partial or complete blindness. Millions of people are affected by this pathology, particularly in Third World countries. A surgical procedure known as Manual Small Incision Cataract Surgery (MSICS) [START_REF] Guo | [END_REF] requires only basic technology. This technique requires a small incision (around 5 mm) to extract the lens in a single piece. The eyeball and the lens must be simulated because they both undergo large deformations and high stresses. The heterogeneity of the lens must also be taken into account since the nucleus is stiffer than the periphery. The lens is modeled with 1113 nodes and 4862 tetrahedra to simulate the lens extraction, whereas the eye contains 1249 nodes and 3734 tetrahedra. The center of the lens is 5 times stiffer than the periphery. The meshes of the organs and the incision have been generated offline. The lens is removed with the help of deformation of the eyeball and friction with the surgical instrument (Figure 2.21).

To account for both the multiple contacts and the heterogeneity of the lens, we used our asynchronous preconditioner to ensure convergence. The preconditioned CG required an average of 11.6 iterations to converge to 10 -5 , despite strong deformations and heterogeneity. Using our method, we managed to maintain a computational speed from 18 to 25 FPS. Within a single time step, the distribution of the computation time was: 40.56% for the free motion and 44.69% for the corrective motion.

Application to liver resection

Nearly 100,000 European citizens die every year of cirrhosis of the liver or liver cancer. Surgical procedures remain the options that offer the maximum success rate against such pathologies. Although other simulators have been developed [START_REF] Bourquain | [END_REF]; Lamadé et al. (2002), our approach is based on patient-specific data. Meshes of the organs are obtained from a semi-automatic segmentation of a CT (see [START_REF] Soler | [END_REF] for details). We simulate 5 deformable bodies in interaction (liver, stomach, colon, intestines, and diaphragm). Each organ is composed of several hundred nodes and thousands of elements with complex shapes composed of several thousand triangles (Figure 2.22). Therefore, an essential issue concerns collision detection which is performed by the method introduced in Allard et al. (2010).

This simulation runs at 25 FPS, including during cutting phases. The distribution of the computing time was: 27.77% for the free motion, 11.01% for collision detection, 22.48% for the constraint motion. The preconditioned CG requires on average 5.82 iterations to converge. Finally, by taking into account the mechanical coupling between the contacts, we managed to produce consistent haptic feedback. For instance, users can feel the stiffness of ribs behind the liver by applying contacts to the soft organ. 2011). The brain is a relatively soft tissue, and most existing simulations rely on explicit integration. Nevertheless, during the surgery, the accuracy of the surgeon's interactions with the tissues is fundamental, as an error of a few millimeters may often have dramatic consequences for the patient. The primary needs of such simulation are precisely modeling and topological modifications, high heterogeneity in the deformation modeling (the tumor is often stiffer than the brain tissues), and real-time computation. We show that our method meets these two issues and may be beneficial to improve the realism of existing brain simulations. The brain is modeled as a heterogeneous deformable body, composed of 1, 734 nodes and 7, 680 linear tetrahedral elements. The tumor is 20× stiffer than the brain. The preconditioner is updated every 5.6 step during the simulation, and a new topological modification appears every 5.5 simulation step, affecting 24 nodes. A total of 553 modifications are performed, and the method remains stable with an average of 5.70 iterations to solve the linear system. The collisions and self-collisions are correctly solved while processing the modifications, and cut parts can instantaneously be separated upon contact with the instrument. Finally, we achieve between 20 and 40 FPS, and the method remains interactive.

Perspective

Perspectives related to numerical solvers and interactive simulations are vast and essential for my research project. Compared to Saupin et al. (2008b) in 2008, the GPU-based parallelization and the asynchronous preconditioning technique introduced in the previous section already allowed to increase the order of magnitude of the size of FE meshes (from approximately 250 to 2000 nodes) that can be simulated in real-time (i.e., 50-100 Hz) with an implicit integration time and interaction constraints. However, despite the important speedups already obtained, the bottleneck remains the computation of the compliance matrix W, which represents up to 60% of the computation time on representative simulations.

Real-time FE simulations for large-scale problems using precondition-based contact response and isolated DOFs constraints: We are developing a numerical method for interactive real-time simulations. The method relies on the asynchronous preconditioner to solve the linear system with multiple right-hand sides. Instead of assembling the compliance matrix in the constraint space as done in Courtecuisse et al. (2014a), we solve the problem in the isolated DOFs space corresponding to the unit DOFs concerned by the constraints. The isolated DOFs space is particularly well suited for GPU parallelization and to reuse information between time steps.

We propose a new strategy to improve the direct solver's performance while building the contact compliance matrix by exploiting the sparsity in the Jacobian matrix (right-hand side) and the factorized system matrix (left-hand side). An efficient GPU-based solver is under development to parallelize the problem. The method can be generally applied in contact resolution and is well suited for an asynchronous preconditioning technique that enables realtime computations for large-scale problems. In addition, based on the nature of the asynchronous scheme, which implies reusing the factorized system, we propose a strategy to share the resolution information between consecutive time steps and benefit from a further speedup. We evaluated the performance in different applications and compared it with existing solutions on the CPU and GPU (Figure 2.24). Preliminary results already show significant improvement compared to the methods introduced in the previous section. With this approach, the computation time of the compliance matrix is drastically reduced. Therefore, we envision the possibility to simulate FE meshes composed of 10000 to 20000 nodes at 50-100 Hz, in the same conditions. Large scale simulation with domain decomposition and Model Order Reduction: To decrease the computational expense, we will investigate Domain Decomposition Methods (DDM). The idea is to split the domain in a set of smaller independent equations that can be reduced and solved independently, allowing this way to significantly decrease the computation expense (Figure 2.25). On top of that, it should open new possibilities to update the domains at various frequencies, according to the desired accuracy of each part. DDM can also be extended to take advantage of the high potential of parallelism provided by recent processors.

Large scale simulation with Model Order Reduction: In collaboration with the Defrost team in Lille, we will investigate the possibility of using projection-based model order reduction. We will use a method of snapshot-Proper Orthogonal Decomposition (POD). The deformation of the organs will be expressed in a low-dimensional space spanned by a reduced basis. In this method, the reduced basis is constructed from a snapshot space containing the most informative movements. One key aspect is the inherent parallelism of the POD approach, and efficient parallelization strategies will be developed to reduce the computational expense of the simulation without scarifying the accuracy and guaranteeing stability. The method is under development, but we believe it will provide a compromise between the quality of FE elements after the cut and the respect of the cutting path performed by the user. More importantly, we believe that the mechanical operation can be significantly optimized since the topological operations, based on element duplication, enable predicting the filling pattern and the matrix modifications after the cut. This information should be very relevant to accelerating the simulation after the cut and improving the stability of the mechanical computations.

Conclusion

This chapter introduced a set of methods for real-time biomechanical simulation of soft anatomical structures, relying on an implicit time integration method. The proposed paradigm relies on an asynchronous preconditioner updated at low frequency, significantly reducing the number of iterations in the CG method. It also improves the contact response process by taking into account the mechanical coupling between contact points. We also extended the approach to handle topological modifications without compromising the interactivity of the simulation. The method is particularly beneficial for heterogeneous structures.

Chapter 3

Intraoperative registration of biomehcanical models

Introduction

This chapter is dedicated to the registration of biomechanical models with intraoperative data. Augmented Reality (AR) is becoming an increasingly helpful tool for guidance and navigation. It can be used to display the preoperative planning on top of medical images and add important information such as the visualization of tumors that are usually difficult to see in the intraoperative context.

Recently, biomechanical models have been used for their ability to regularize the ill-posed non-rigid registration problem (see [START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF] for a detailed discussion). Biomechanical models are built from preoperative images and used to extrapolate the displacements of organs (Figure 3.1), allowing this way to update the preoperative information to the intraoperative configuration, even where no intraoperative data are available. For that, patient-specific simulations must be employed, involving at least four additional challenges:

1. The geometry of the organ models must be acquired from medical images, and this process is still not automatic, requiring complex segmentation and mesh generation.

2. The material properties of living tissues need to be characterized. Such properties are patient-specific, and it is difficult to predict the mechanical behavior of in vivo tissues.

3. Boundary conditions are essential since they significantly influence the deformation. However, ligaments or contact with the environment is challenging to see in medical images.

4. Advanced models need to be employed to reproduce the mechanical response of organs with high accuracy, including, for instance, heterogeneous stiffness or multiple tissue materials.

Up to now, such generic tools providing surgical assistance with robustness and accuracy do not exist. Indeed, the non-rigid registration of intraoperative images remains an open problem currently addressed by a whole scientific community. To limit the scope of this chapter, we mainly focus on liver cancer surgery and brain shift compensation.

Related Works

The use of biomechanical models in the operating room is an active research topic (see [START_REF] Martin | Augmented Reality in Medicine: Systematic and Bibliographic Review[END_REF] for a recent survey). In the following, we briefly survey relevant methods in the context of preoperative planning, and intraoperative guidance Carter et al. (2005); Maier-Hein et al. (2013).

Liver motion and non-rigid registration

Primary cancer of the liver and metastases is the sixth most common cancer and the third most frequent cause of cancer death, resulting in 600,000 deaths annually. Hepatic resection is one of the few potential curative treatments for large tumors (i.e., > 3 cm), but the associated survival statistics are often limited.

The liver is a parenchymal organ with potentially deep tumors that are difficult to reach, and during the respiration cycle, the liver undergoes significant motion in the abdominal cavity, mainly induced by the motion of the diaphragm Kang et al. (2012). In Rohlfing et al. (2004), the craniocaudal translation is reported as the dominant component of the liver motion during respiration (12-26 mm), compared to anterior-posterior (1-12 mm) and lateral (1-3 mm) motions. However, it is stated that purely rigid transformation is not sufficient to model liver motion during the respiration cycles, as some parts of the organ undergo significant displacements (up to 35 mm), which cannot be captured with rigid registration.

Attempts of AR have been realized in laparoscopy to improve the visualization of the screen. It appeared as an ideal interface between the numerical reconstruction and the surgeon [START_REF] Kang | Stereoscopic augmented reality for laparoscopic surgery[END_REF]; Buchs et al. (2013). A 3D visualization and surgical planning software tools are proposed in Soler et al. (2014) to overlay a 3D preoperative model on top of the stereoscopic camera of the Da Vinci robot. Initially, the model registration was manually performed, but Bernhardt et al. (2016) proposed an automatic solution estimating the location of the tip of the endoscope in the volume image data.

Since pneumoperitoneum represents the primary source of shape variation, some work aimed at precomputing the deformations to reduce the initialization problem to a rigid model-based alignment Sánchez-Margallo et al. ( 2011). The method proposed in [START_REF] Bano | [END_REF] relies on biomechanical models to simulate the pneumoperitoneum as external pressures are applied inside the abdominal cavity. It is mainly used to optimize the trocar's placement rather than solving a registration problem. Similarly, a deformed mesh is computed using a mass-spring-damper model Nimura et al. (2015), and the insufflation pressure is applied to the volume as boundary conditions. Recently, a database-based was proposed to estimate the deformation of the organ under pneumoperitoneum [START_REF] Johnsen | [END_REF]. The database is built from both intra-operative images and preoperative segmentation. It is intended to be used as an atlas of deformations. However, due to many unexpected behaviors, these methods fail at reducing the problem to a pure rigid alignment and can hardly estimate the intra-operative deformation.

Other approaches rely on intraoperative imaging to directly estimate the deformation. Clements et al. (2008) proposes a rigid alignment approach based on salient anatomical features, extracted both in preoperative images and endoscopic data. In Clements et al. (2017) authors provide evidence of benefits of deformation correction with an evaluation across 20 patients. Nevertheless, the method remains sensitive to unavoidable occlusions of the liver (for instance, by the surgeons' hand), and the method may provide miss correspondences for large deformations. Besides anatomical landmarks, organ silhouette was recently considered by exploiting the organ's rigidity [START_REF] Collins | [END_REF]. Finally, Stoyanov et al. (2010) proposed to reconstruct a 3D point cloud located on the surface of the organ from a laparoscopic stereoscopic camera. The practical value of the proposed method is demonstrated by reconstructions of various in vivo images of robotic-assisted procedures.

Biomechanical models allow regularizing the ill-posed non-rigid registration problem and provide extrapolation between image-based constraints in areas where no data are available [START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF]. [START_REF] Suwelack | [END_REF] proposed an electrostatic-elastic problem formulation for the registration of a physically-based model derived from preoperative mesh and intra-operative surface computed from a laparoscopic stereo camera. An elastic body is electrically charged to slide toward an oppositely charged rigid shape. However, the method requires that at least 50% of the organ surface be captured by the camera, which is generally impossible. In Haouchine et al. (2013), a biomechanical model of the liver based on preoperative CT images is used to estimate the actual position of a tumor based on tracking of surface features extracted from a laparoscopic camera during the intervention. Elastic springs are used to couple tracked features and the biomechanical model. Plantefève et al. (2016) exploits salient anatomical features, identifiable in both preoperative and intraoperative images of the liver. The method was used to display in AR the internal structures of the preoperative scan on top of the intraoperative view obtained from a laparoscopic camera. In addition, the organ silhouette brings essential information to perform the registration process [START_REF] Haouchine | [END_REF]. However, the organ's contour can be incomplete and complex to extract due to inter-patient variability and may not be sufficient to solve the problem. To solve this problem, [START_REF][END_REF] proposed to combine the contour information with shading cues to estimate the deformation.

Additional imaging modalities can be used in the operating room to facilitate the initial alignment and limit uncertainties related to endoscopic images. Intraoperative ultrasound has been considered in Dagon et al. (2008) to register the liver vessel tree on a three-dimensional liver model. The registration of the vessels drives the deformation of the whole organ and thus approximates the deformation. Oktay et al. (2013) suggested exploiting intra-operative CT scans after insufflation as an additional constraint to drive the simulation.

Recently Bernhardt et al. (2016) proposed a method to scan the tip of the endoscopic camera to find its rigid pose w.r.t the intraoperative scanned organs. These methods provide accurate registration, but intra-operative scans are, for the moment, not available in standard operating rooms. 2D slices to 3D volume registration problems have also been considered Ferrante and Paragios (2017). In this case, it has been shown that the mechanical parameters, as well as boundary conditions, can be at least as crucial as the constitutive law itself Bosman et al. (2014). [START_REF] Peterlik | Modelbased identification of anatomical boundary conditions in living tissues[END_REF] showed that the mechanical action of boundary conditions could be directly extracted from the biomechanical registration process. Similarly, Coevoet et al. (2014) formulated the non-rigid registration problem as an inverse mechanical problem. The method provides a set of boundary conditions bringing the FE model to the desired geometrical deformations.

Brain-shift compensation methods

During brain tumor ablation surgery, accurate localization of the brain tumor is essential to ensure its total resection and reduce the morbidity of surrounding healthy tissues. Images of the patient's brain are acquired before surgery and used by image-guided neurosurgical navigation systems to assist the surgeon. However, due to the intraoperative deformation of soft tissues, called "brain-shift", these images no longer correspond to the brain morphology of the patient throughout the whole procedure.

The causes impacting the amount and direction of brain-shift are numerous: patient positioning during surgery, craniotomy size and dura opening, loss of Cerebrospinal Fluid (CSF), actions of the surgeon, swelling due to drugs, anesthetics or edema, etc. Even if the brain shift mainly occurs in the direction of gravity, it is challenging to predict, and the surgeon, following his experience, thus estimates brain deformations. Magnitude and direction of brain-shift have been the subject of several studies reporting cortical displacements up to 20 mm and subsurface movements up to 7 mm [START_REF] Gerard | Brain shift in neuronavigation of brain tumors: A review[END_REF]. The main methods to compensate for brain shift occurring during tumor ablation surgery are pure image-based and biomechanical simulation-based techniques.

Image-based registration methods: [START_REF] Nimsky | [END_REF] presented a rigid registration between pre-and intraoperative MR images, using fiducials placed around the craniotomy. A non-linear registration method based on mutual information is introduced by Hastreiter et al. (2004). Regardless of the accuracy, the use of such intraoperative MR devices is cumbersome. The acquisition procedure is complex and considerably increases the operating time: the patient must be transferred to the scanner, and specific surgical tools are required due to the magnetic field generated. Moreover, these devices are expensive and require large dedicated operating rooms. For these reasons, this intraoperative modality is rarely used in clinical routines. [START_REF] Pereira | [END_REF] proposed to use pre-and intraoperative Cone-Beam Computed Tomography (CBCT) enhanced with intravenous injection of contrast products. The volumetric deformation vector field is computed by combining a rigid registration of the skull and then an elastic transformation based on vessels, ventricles, and image intensities. However, CBCT images suffer from poor soft tissue contrast and are thus rarely used for brain tumor surgery.

Several methods rely on the acquisition of intraoperative US images. Ultrasound scanners are light systems, portable, already available in most operating theaters, and cost less than 10% of a classic intraoperative MR device. In addition, intraoperative US acquisitions do not necessitate essential changes in the operating procedure and are compatible with other surgical equipment. On the one hand, deformations of soft brain tissues can be visualized using B-mode US imaging. These images can be registered with preoperative MRI using mutual information [START_REF][END_REF], cross correlation (Rivaz and Collins, 2015) or linear correlation of linear combination [START_REF] Fuerst | [END_REF].

On the other hand, Doppler US imaging provides flow visualization. Vessels around the craniotomy region can thus be viewed using this US modality. [START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF] proposed a modified Iterative Closest Point (ICP) algorithm to rigidly register the blood vessels extracted from preoperative MR Angiography (MRA) to the ones extracted from intraoperative Doppler US acquisitions. A non-linear transformation using a thin-plate spline is also presented. Finally, this rigid registration approach was validated in a clinical study including seven patients [START_REF] Reinertsen | Intra-operative correction of brain-shift[END_REF]. After defining matching features between pre-and intraoperative MRI, displacements are imposed on a brain biomechanical model to perform the registration. However, as stated before, this MR modality is rarely used in clinical practice during surgery.

Several groups proposed to track the exposed cortical surface intraoperatively. On the one hand, Sun et al. (2005a), and De [START_REF] Lorenzo | Volumetric Intraoperative Brain Deformation Compensation: Model Development and Phantom Validation[END_REF] introduced models where the displacements of the brain surface are directly driven by the motions acquired with stereo-cameras. For a better integration to the surgical process, Sun et al. (2005b) suggested attaching these cameras to the stereoscopic operating microscope. Their system was evaluated on ten patients. On the other hand, Sun et al. (2014) precomputed an atlas of brain deformation, calculated using a distribution of boundaries and deformation-inducing force conditions. The brain shift is then determined using an inverse problem approach which linearly combines solutions from the atlas. This method has been recently validated by Miga et al. (2016) using a clinical study of over 16 patients. A limitation of these methods is the assumption that the whole 3D deformations can be extrapolated by the biomechanical model from the exposed brain surface only. However, according to Wittek et al. (2007), the registration accuracy improves when data are collected from the exposed and unexposed brain's surface. Moreover, no subsurface features (e.g., the ones visible with MR or US images) can be considered. Mohammadi et al. (2015) thus proposed to combine the tracking using stereo cameras with Doppler US imaging. The exposed cortical surface, as well as the vascular tree close to the tumor, are registered. For that, a constraintbased biomechanical simulation is used. However, this requires bringing two imaging systems (i.e. the stereo cameras and US scanner) in the operating room and therefore complicates the clinical process. Finally, Bucki et al. (2012) presented a similar model-based method. Only the blood vessels are registered, and a single Doppler US acquisition is necessary. In addition, results on one surgical case are provided.

Methodology

The contribution of this chapter is a unified constraint-based formalism and solving process combining intraoperative image data with biomechanical models.

Problem definition: The non-rigid registration of a preoperative 3D model S with intra-operative data is an ill-posed problem known to be very difficult to solve. It can be formalized as finding the 3D non-rigid transformation T mapping S into the image frame. However, since intraoperative data are usually incomplete, errors can only be measured with the information that can practically be extracted from the operating scene (image or sensors). Let m be the set of data extracted per-operatively, the non-rigid registration problem can be synthesized with: This process usually admits a unique solution for rigid transformation scenarios (without considering noise or other pairing issues). However, the rigid assumption is not valid for our purpose, and the deformations generate a problem having an infinite number of solutions.

To decrease the number of solutions, we assume an elastic behavior of the organ with preservation of the volume after the deformation. Therefore, we propose to rely on a non-linear biomechanical Finite Element model coupled with a set of constraints H to solve the registration problem of equation (3.1). This leads to a non-linear problem whose solution is given by the positions of the physical model providing the equilibrium between the internal forces and the external forces applied by the constraints:

F(p) + H(m, p) = 0 H(m, p) = δ (3.2) (3.3)
where F provides the internal forces of the FE model, and δ is the violation of the constraints, as described in the previous chapter. Lagrangian Multipliers λ are to impose displacements and enforce constraints of equation (3.3). The number and the direction of Lagrangian Multipliers will be described below; it depends on the position of the models p and the data extracted from intraoperative images m.

Note that equation (3.2) and (3.3) corresponds to the static equilibrium of the dynamic equation (2.22). However, in this chapter, the simulation time t is not considered anymore. Instead, the simulation steps correspond to the iterations of the iterative registration procedure.

At this point, it is important to underline that the solution of the equation (3.2) is still not unique for mainly two reasons:

1. The internal force F(p) depends on the initial configuration of the model. Indeed, from a numerical point of view, the rest configuration should correspond to the organ's shape where the FE model is not subject to any external forces. However, since FE meshes are generated from segmentation, it depends on the organ's configuration when the preoperative image has been acquired. An important issue is to estimate an initial transformation D(S) to cancel the initial stress of the FE model before running a simulation.

2. H(m, p) is a nonlinear function used to impose displacements to the model's position p in order to fit per-operative data m. However, the correspondences between m and the material points on the model are usually not known. An iterative process is employed to deform the model toward the shape-fitting image data and minimize the model's internal energy. However, this process is usually sensitive to the initial positions, and the method may converge towards the closest local minimum, which does not necessarily correspond to the global solution of equation ( 3.1).

Finally, the evaluation and validation of the solution is also a scientific challenge that needs to be addressed Thompson et al. (2015). Indeed, the problem is complicated since the ground truth is usually unavailable due to imaging and technical limitations.

Segmentation and meshing and initialization

Preoperatively, an anatomic patient-specific model is built from the images acquired before surgery. As preoperative images are often acquired one day before surgery, several hours are available to build the model. However, for clinical usage, the generation of the patient-specific model should be as automatic as possible.

The model is usually constructed using either contrast-enhanced CT data or volume MRI sequence. The image is first segmented (using, for instance, the semi-automatic methods available in ITKSnap2 ), then a tetrahedral mesh is obtained from the segmented maps (for instance, using CGAL3 ). In our applications, the target resolution of the FE mesh is around 2000 nodes to enforce fast computation times.

The FE model should correspond to a state of equilibrium between external and internal forces, but the segmented organ is inevitably subject to loads when the image is acquired. A typical example concerns gravity. Most biomechanical models ignore gravity and directly simulate the initial mesh with null stress. However, this will yield different results, especially for living tissues that are usually soft and undergo large deformations. When the gravity force is applied to a model, it will shrink unless other conditions are accounted for (boundary conditions, internal forces, muscle contraction. . . ) or other external forces or pressures. Therefore, a rest shape must be determined so that the model converges to the targeted equilibrium state when loads are applied.

Formally speaking, a transformation D(S) mapping the initial segmentation S to a shape minimizing the initial internal stresses is needed. In addition, the reference frame of the preoperative and the intraoperative images may be significantly different as they usually come from different devices and imaging modalities. A rigid initialization R(D(S)) is therefore needed. This initial transformation can be estimated automatically (using, for instance, stereotaxic markers) or provided by the user at the beginning of the operation. However, since the rigid transformation is estimated intraoperatively, it must be available within several seconds for clinical usage.

Geometrical binding

Vision-based algorithms are traditionally used to extract geometrical primitives (2D/3D positions, partial triangular mesh. . . ) from intraoperative images. These data provide the intraoperative shape (called target) on which the preoperative model (called source) must be registered.

A significant difficulty is related to the fact that it is complicated to relate the position of the resulting data to the FE model. Indeed, the correspondences between image data and material points on the model are usually not known (Figure 3.3). For instance, a 3D point cloud on the organ's surface may be obtained from a laparoscopic stereoscopic camera. However, identifying the corresponding visible part of the model is a difficult task since the visible part is usually limited, noisy, and the organ is already deformed during the image acquisition.

? This problem is addressed with an iterative process taking as input the target data m and the source positions p of the collision mesh. Indeed, in practice, the search is usually performed with positions p of a dedicated collision model, either to save computation time or to have a specific rep-resentation of the structures on which constraints will be applied (vessels, contours. . . ). Nevertheless, as explained in the previous chapter, a linear relation J can be precomputed using barycentric coordinates (see equations (2.16) and (2.17)) to transfer the constraint equations to the mechanical DOFs p. For the sake of simplicity, we now omit this relation, considering that the constraints are directly applied to the mechanical model.

As in the previous chapter, the constraint equations are linearized at each iteration H(p, m) ≃ ∂H ∂p = H, providing the directions along which Lagrangian Multipliers are computed. The linearization (corresponding to the collision detection in the previous chapter) corresponds to the pairing of the control points m with material points on the model. The iterative closest point (ICP) method Rusinkiewicz and Levoy (2001) is employed to associate the target points with their respective closest4 points on the source mesh. The linearization is performed at the beginning of the simulation step (corresponding to the iteration t) using the current position of the model and assumed constant during the solving process.

Outliers and geometrical filtering:

The quality of the registration strongly depends on the quality of the pairings performed during the binding process. Indeed, incoherent or antagonistic pairings (called outliers) can appear due to the noise (image artifacts, missing data. . . ) and must be ignored. Depending on the application, the following filters may be used:

Unique pairing: First, we enforce that control points m are associated with a unique element on the model; otherwise, it would lead to overconstrained problems. When several vertices of m are projected onto the same geometrical element of the mesh, only the nearest one is kept.

Filtering with distance: A pairing is considered as an outlier, and then ignored, if it does not satisfy the following condition:

d j ∈ [ d -d t ; d + d t ] (3.4)
where d is the median distance computed over all pairings and d t , a threshold. In addition, a filtering criterion d max can be added in order to reject pairing above a given distance d max .

Filtering with normals: When image data are sufficiently dense, the normals of the surface can be estimated using a simple least square plane fitting as implemented in the Point Cloud Library5 . In this case, the method can be improved by ensuring that the normals associated with image data have the same orientation as the FE surface (Figure 3.4). 3.4b the distance criterion is not sufficient since all the control points are associated with the same side of the object and the constraints cannot be satisfied. 3.4c, a set of normals aiming inside the surface are computed. 3.4d the control points are associated with the closer triangle whose normal is oriented in the same direction.

The method is improved by using a cubic Bézier interpolation of the source surface as described in Vlachos Jörg Peters et al. (2001). It provides a smooth description of the triangulation allowing for a continuous sliding of the constraints between edges and triangles, which helps stabilize the registration. The barycentric coordinates of the closest point on the cubic interpolation of the triangles are determined with the Newton-Raphson algorithm. At each time step, the control points m are associated to their respective closest points on the Bézier path, and the corresponding normal is evaluated on the Bézier interpolation.

Constraints definition

This section details the various constraints employed for our applications. At this stage, the control points m of the target shape are associated with material positions p on the biomechanical model. Lagrange multipliers are used to impose constraint forces on the FE model to cancel the violation of constraints.

Bilateral constraints Ω

Bilateral Constraints are holonomic constraints Ω(p) = 0, used to impose the relative positions of constrained points on the orthogonal plane defined by the constraint normal. For a given control point i, the violation of constraints δ i is computed with:

δ i = n i • (m i -p free i ) (3.5)
with • being the Euclidean dot product and n i is the normal of the constraint.

Depending on the desired behavior, n i may be computed in several ways. In can be either the normal of the Bezier source surface at position p i , or the direction of the pairing n j = m i -p i ||m i -p i || , or simply the unit vectors ⃗ x, ⃗ y or ⃗ z. At the end of each simulation step, the corresponding constraint force λ i , applied in the direction of n i , is computed so that no violation remains ( δi = 0) in that direction. It means that the vertex i, with position p i , will be located on the tangential plane given by n i (Figure 3.5). This formulation allows the control points to "slide" on the surface of the FE mesh to stabilize around the configuration, minimizing the energy and satisfying the constraints. The combination of two bilateral constraints forces the points to follow a line at the intersection of the two planes (Figure 3.5b). As introduced in the next section, this is particularly relevant for imposing the target points to "slide" along vessels. Similarly, three orthogonal bilateral constraints force the point p i to be located at the exact 3D location m i at the end of the step.

δ f =0 δ i >0 p i 𝑛 m p f δ f =0 δ i <0 p i 𝑛 m p f (

Unilateral constraints Ψ

Unilateral constraints are used to simulate collision on the source mesh. It can be either collision with other FE models, the surrounding environment, or even collision with medical instruments segmented in the images. Unilateral constraint satisfies the Signorini conditions λ ⊥ δ as explained in the previous chapter (see section 2.3.4.1).

The violation is computed the same way as in equation (3.5). At the end of the resolution, the constraint node i is located on the positive side of the constraint. Therefore, the additional difficulty compared to bilateral constraints is related to the fact that the constraint's normal n i must be oriented in the positive direction and remain consistent during the iterations of the registration procedure.

Projective Constraints φ

Projective constraints can be used to register the contours of the source with the target when using protective cameras. However, the first difficulty is to find the primitives of the model that correspond to the visible contour of the model in the image (Figure 3.6). Indeed, the 3D source model is a close surface (without outline), and the model's contour must be computed according to the projection matrix P c of the camera. Contour extraction: The algorithm6 is performed in two steps. We first select all the front-face triangles with respect to the camera direction dot(n, z) > 0, with n being the normal of the triangle and z the direction of the camera. However, this test is not sufficient since the source model is not convex. Indeed, many triangles remain selected (red triangles in fig. 3.6a) despite not being visible from the camera since other triangles are closer (green triangles). An additional intersection test is performed to keep only the visible triangles. Let a be the camera position and p (i) the set of points that belong to the remaining front-face triangles (red and green). For each position p (i) , if the ray --→ p (i) a intersects any other front-face triangles, all the triangles connected to p (i) are discarded.

Finally, the visible contour of the model (yellow in 3.6a) is obtained by selecting all the edges that are connected to only one visible and one non-visible triangle.

Contour pairing: Once the model's contour is obtained, it is projected in the image using the projection matrix P c . The Iterative Closest Point method is then used to bind each point of the contour m, with the closest point on the projected contour of the model.

Finally, the constraints are projected back in 3D and applied to the model. The constraints must only impose displacements in the 2D image's plane, leaving the model-free in the direction of the camera's depth. This way, the motion of the model will be driven by internal elastic forces to reach equilibrium with constraint forces, i.e., if the contour of the model is smaller than the image's contour, stretching forces will tend to bring the model closer to the camera, whereas compression will push away the model.

Constraint Solving

The linearization of equations (3.2) and (3.3), provides the following KKT system in a quasi-static scenario:

Ax + H T λ = 0 Hx = δ (3.6) (3.7)
Note that although we only consider static equations for the registration, the dynamic terms of equation ( 2.12) are usually conserved during the registration process, i.e. A = (1 + hα)M + h(h + β)K as in equation (2.12). Indeed, dynamic terms increase the stability of the registration procedure, and since we only consider the converged state (i.e., when the velocity update dv = 0), the solution of the dynamic and static equations are identical.

A new pairing procedure is called at each step to associate the control points m i with their respective closest point p i on the model. For each control point, a unique bilateral constraint is defined to perform the registration. Although for a single constraint, this would move the point p i to the exact location m i , it is not necessarily the case for multiple constraints. Indeed, when several constraints are applied simultaneously on the model, the constraints are coupled through the compliance matrix W. Therefore, at the end of the simulation step, p i and m i are located on the tangential plane defined by the normal n i , but the cumulative effect of other constraints may move the solution away in lateral directions (Figure 3.7). Therefore, the model tends to "slide" toward the shape, minimizing the energy necessary to enforce the constraints.

During the Gauss-Seidel iterations, the constraint equations are either activated with a non-zero force or deactivated if the violation is canceled by other constraints (Figure 3.8). As a result, only the constraints necessary to cancel the violation are activated, and λ minimizes the energy to register the model.

The following section will show that this observation can be used to identify the boundary conditions acting between the preoperative and intraoperative configurations. Finally, the iterative registration process is stopped when the deformable body is stabilized, i.e. when the variation of velocity between consecutive iterations t and t + 1 is lower than a threshold ||x (t+1) -x (t) || < ϵ.

Outliers and mechanical filtering:

The geometrical filters introduced previously may not be sufficient. Additional mechanical filters are employed to stabilize the process and facilitate the convergence.

Force clipping: During the solving phase, the constraint forces are computed as follows:

λ ∈ [-f max ; f max ] (3.8)
with f max a parameter defining the maximal force applied on the constraints. Since the constraint forces are bounded, all the constraints might not be satisfied at the end of the resolution ( δ ̸ = 0). However, the cumulative effect of all the constraints will create the general motion/deformation of the model. Strong forces are only necessary if the registration constraints require significant energy to act against the internal elastic forces of the deformable body, which is most of the time attributed to outliers.

Image compliance: In the above formulation, the data extracted from images are never impacted by the biomechanical model, but instead, the model deformed to fit the data. However, as stated in the introduction, the data may be noisy, and their positions are subject to inaccuracy.

Figure 3.9: The white contour is directly given by image-based tracking where the top right control point has been intentionally displaced above the liver. The orange shape is the wireframe surface of the FE mesh and the green shape is the resulting polyhedron using soft constraints.

A compliance factor is associated with the constraints to improve the robustness of the method:

(W + W soft ) λ = δ (3.9)
where W soft is a diagonal matrix whose coefficients w soft are chosen according to the mechanical parameters of the model and the confidence of the tracking method. If w soft = 0, image data m are not influenced by the constraints, whereas the higher w soft , is chosen, the more the image data will be projected onto the model rather than deforming it. Therefore, both the biomechanical model and the data m can be deformed during the constraints resolution (Figure 3.9).

Contributions

In this section, our methodology is applied and tested with several image modalities and clinical applications. The goal is to bring our methodology as close as possible to the Operating Room (OR). 2015) we introduced a method to compute the rest shape for highly deformable models of the brain. We proposed an iterative method that estimates the rest shape of the model, i.e., canceling gravity and the cerebrospinal fluid pressure (CSF) on the surrounding brain. The method is based on Sellier (2011), with the addition of collisions with a surrounding surface.

Rest shape (Figure 3.10a) was found with an average error equal to 0.746mm. The figure 3.10b shows that the highest errors are located on the surface of the model, close to contact areas. Nevertheless, inside the volume, the errors are low, and the shape is well interpolated. The figure 3.10c shows the von Mises stress repartition on a slice of the mesh at equilibrium. The blue area corresponds to the fixed cerebellum part: stress is quasi-null here. Above this region, there is a red area with high stress due to the gravity acting on the mesh. Finally, the heterogeneous region of the tumor also undergoes high stress because of its higher stiffness.

3.4.1.2 Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation Planning and guidance are based on preoperative images during brain tumor ablation, which does not account for the brain shift. However, this deformation is a significant source of error in image-guided neurosurgery and affects the accuracy of the procedure. In Morin et al. (2016Morin et al. ( , 2017d)), we introduced a constraint-based biomechanical simulation method to compensate for brain shift that integrates the deformations of the blood vessels and cortical surface, using a single intraoperative ultrasound acquisition. An overview of the method is shown in Figure 3.11. The method combines a biomechanical model built from preoperative images and image-based constraints extracted from intraoperative US images. Before surgery, the patientspecific FE model of the brain is constructed from preoperative MRI images. The model involves the morphology of soft tissues and the geometry of the vessels located around the tumor. During surgery, localized Doppler and Bmode ultrasound images are acquired directly in contact with the brain. The vascular tree and the footprint of the ultrasound probe are then extracted from these intraoperative images. A biomechanical simulation is then performed to compensate for the brain-shift deformation. We introduce several types of constraints allowing to 1) model contacts between the brain and the dura, 2) register the pre-and intraoperative vessels, 3) constrain the cortical surface under the footprint of the probe. Finally, the preoperative MRI images are updated using the displacement field calculated from the biomechanical model of the organ.

The method has been evaluated on five retrospective cases. Data were collected by the SINTEF Medical Technology Institute at St. Olav University Hospital (Trondheim, Norway) [START_REF] Reinertsen | Intra-operative correction of brain-shift[END_REF]. For each clinical case, T2-FLAIR and Angiographic MRI images were acquired before surgery. During the procedure, ultrasound images were acquired. These data were collected through a clinical study approved by the local ethics committees, and patients' consent was obtained before the procedure.

Our method has been evaluated and compared to the rigid registration of [START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF][START_REF] Reinertsen | Intra-operative correction of brain-shift[END_REF], available in the open-source platform for image-guided therapy CustusX Askeland et al. (2016). The first difficulty in obtaining quantitative results was to define reliable landmarks visible in both image modalities. For each case, 5 to 9 landmarks were first identified in blood vessels bifurcations, by two operators, on both the preoperative MRA and intraoperative Power Doppler US images. These landmarks were used to define a target registration error. In addition, anatomical structures such as sulcus were delineated by a clinician in the MRI and B-mode US acquisitions, allowing a measurement that is independent of the data used by the method to enforce constraints on the model. After registration, MRI is wrapped using the displacement field of FE models and displayed to the surgeon (Figure 3.13). A pointer first shows the borders of the exposed cortical surface (top row) then a deep sulcus bifurcation point (bottom row).

Qualitative observations, especially deep around the tumor, show consistent matching between the wrapped MRI and US images, whereas significant errors can be observed in the initial preoperative MRI, enhancing the importance of the deformation.

Although validated in a post-operative study only, our method seems usable in a clinical context. Indeed, most steps of our approach, particularly the per-operative ones, are performed without the need for interactions of an operator. The generation of preoperative models was performed in approximately 20 mn (excluding the segmentation of the brain). As shown by the figure 3.14, the construction of the soft tissue model includes the mesh generation (FE and collision meshes), the definition of boundary conditions, an inversion of the stiffness matrix (denoted K 0 -1 ). Except for the segmentation of the structures, all these steps are performed automatically. The most expensive part is segmenting the vessels in the ARM image and their skeletonization that takes about sixteen minutes. Manual interactions are still needed, but as said before, this is acceptable in the clinical workflow as these operations can be performed on the day before the surgery.

Preoperatively, the additional time of our method is shown in the figure 3 for simulation, and 10 seconds to warp the image. In addition, the mode B and Doppler ultrasound images being acquired simultaneously, the extraction of the footprint of the probe and the vessels are performed in parallel.

Given that brain surgery generally lasts for several hours, this overhead seems acceptable for a per-operative usage of our approach. In Chabanas et al. (2018) the method was extended to take into account the resection-induced deformations occurring during the tumor removal procedure. A qualitative evaluation of the compensation is also presented using initial and updated images. An analysis of three cases of surface tumors shows that the biomechanical registration can compensate for up to 63% of the brain shift, with an error in the range of 2 mm (Figure 3.16). While more cases must be considered, these first results have shown the ability of the method to compensate for resection-induced brain shift, without additional treatments, in the case of a peripheric tumor.

Resection-induced brain-shift compensation using vesselbased methods

Physics-based registration with 2D MRI slices

Modern MRI allows for dynamic scanning without any artifacts due to movement. For instance, the MRI MAGNETOM ® Aera SIEMENS 1.5 T allows getting 2.5 acquisitions per second of a given slice. This frequency is sufficient to capture the dynamic behavior of tissues during breathing cycles, but it is for now restricted to only one plan of acquisition, and given the absence of any volumetric information, some anatomical structures might be excluded. The outline of the organ is manually segmented in the first image of the MRI sequence. The tracking method takes as input the initial 2D segmentation and a continuous flow of dynamic MRI slices. An optical flow algorithm provided by OpenCV is then used to track the displacement of the liver's contour in the dynamic sequence (Figure 3.17d). Based on the position and orientation of the MRI slices, a set of sparse 3D control points located on the liver's surface can be reconstructed. However, due to off-plane motion, the control points provided by the MRI slices cannot be associated statically with the same material point in the liver. Therefore, we proposed a 2D/3D registration method combining 3D FE models with the 2D images data. The method provides a 3D extrapolation of the dynamic of the organs observed in 2D MRI slices. Given an average porcine respiratory rate of 15-20 breaths per minute, we acquired sequentially MRI dynamic data along nine different orientations chosen randomly (Figure 3.18). Assuming the breathing motion is periodic, slices have been synchronized manually, and the contour of the organ has been segmented manually in all the images. It provides a sparse set of 3D points that can be used to evaluate the registration errors.

The error was defined as the Hausdorff distance between the surface of the mesh and the control points. We tested various constraint configurations (i.e., combinations between active registration slices and validation slices).

We showed that 3 orthogonal MRI slices are sufficient to perform an entire 3D registration of the preoperative model. Yet, As the dynamic motion of the liver can only be acquired along a single plan at time, we showed that combining static slices with a single dynamic acquisition already provides acceptable registration error for the breathing motion (Figure 3.19). 2020) we extended the method for 3D registration of female pelvic organs using 2D dynamic magnetic-resonance images (MRI). The aim is to provide better knowledge and understanding of pathologies such as prolapsus or abnormal mobility of tissues. 2D dynamic MRI sequences are commonly used in clinical routines to evaluate the dynamic of organs, but due to the limited view, subjectivity related to human perception cannot be avoided in the diagnoses. The method was applied to the four main structures of the female pelvic floor (bladder, vagina, uterus, and rectum). The original contribution was to combine the 3D FE models, image data, and a priori knowledge of boundary conditions. Boundary conditions include contacts between structures as well as an approximative definition of ligaments and attach structures (Figure 3.20). Although the exact patient-specific modeling of these conditions is complex (since boundary conditions are usually not visible in medical images), we showed that even approximative information improves the solution.

The methodology was evaluated with two patient-specific data sets of volunteers presenting no pelvic pathology, and a sensitivity study is performed using synthetic data. The resulting simulations provide a 3D extrapolation of the dynamic of the organs observed in a single 2D MRI slice, facilitating diagnosis compared to 2D sequences. In addition, the methodology follows a protocol compatible with current clinical constraints presenting potential short-term medical applications in this way. we introduced a method for semi-automatic registration of 3D deformable models using 2D shape outlines (silhouettes) extracted from a monocular camera view (Figure 3.21). The method is composed of the following steps: 1) A direct simulation D(S) is applied to transform the reconstructed model obtained from the segmentation (red) in shape close to the 3D position observed in the image (green). 2) A Rigid transformation (blue) is provided by the user to roughly align the model with the contour of the organ segmented in the image (yellow). 3) Projective constraints are applied to the biomechanical model to fit the organ's contour and finally provide the 3D shape in the camera's frame.

Results show that the final elastic registration can be obtained in just a few seconds, thus remaining compatible with clinical constraints. We evaluated the sensitivity of our approach with both the initial alignment and the silhouette length and shape. Finally, we applied our method to surgical data (Figure 3.22). We first applied our method to the liver registration in la- paroscopy, where other organs partially exclude the contour. Our method rapidly converges to a realistic configuration using only contours located on the top and bottom of the organ. We also applied the same approach for kidney surgery, where the entire contour of the organ was visible in the image. The method was combined with infrared markers in Adagolodjo et al. (2018) to produce an Augmented Reality (AR) system for open liver surgery. The workflow of our method is introduced in the figure 3.23. Preoperatively, a FE model is built based on the CT scan of the patient. Markers are manually placed on the organ's surface after opening the abdominal cavity and tracked in real-time by a set of infrared cameras. Initial registration is then performed, combining contours and markers data to perform the initial alignment. Markers attached to the liver and tracked in real-time by infrared cameras, and used to deform the model. Finally, the video stream of a monocular camera is displayed on monitors in the operating room, overlaid by the deformed model.

The method was tested with both synthetic and ex-vivo samples. In addition, our method was applied in the operating room during a liver resection of a human patient. A biomechanical model composed of 1900 tetrahedral elements has been generated from the preoperative CT, including the segmentation of the tumor and the vessel tree. The cameras have been fixed on a metallic bar, mounted transversely on top of the patient with a direct view above the patient (Figure 3.24). After the patient's opening, seven markers of 0.5mm diameter have been sutured on the patient's liver around the tumor area in approximately 10 minutes. Their position has been chosen randomly to cover the maximum area of the organ visible in the optitrack system. The initial registration (including segmentation of the contour and initial rigid alignment) has been performed in approximately 5 minutes. Then the tracking method was applied for 10 minutes providing convincing registration results even for large deformations and occlusions by the hands of surgeons or when the organ is pushed inside the abdominal cavity. The simulation was interactive with a frame rate above 100 FPS thanks to GPU parallelization. This preliminary study provides promising results to improve the location of tumors and help surgeons plan the ideal resection intraoperatively.

Identification of Anatomical Boundary Conditions

Together with the material properties, the boundary conditions have a significant influence on the mechanical response of the organs. However, their modeling remains an open question, as the connective structures are challenging to identify in standard imaging modalities. 2014) we presented a novel method dealing with the identification of boundary conditions of a deformable organ. As an input, the method requires a set of scans acquired in different body positions. Using the constraint-based FE simulation, the method registers the two data sets by solving an optimization problem minimizing the energy of the deformable body while satisfying the constraints located on the surface of the registered organ. Once the simulation equilibrium is reached, the surface forces needed to satisfy the constraints provide a reliable estimation of location, direction, and magnitude of boundary conditions applied to the object in the deformed position. The method was tested employing two scenarios with a beam object, where the deformations were computed via simulations to have both the deformed shape and surface loads in the target deformation. The data was used as ground truth and compared to the von Mises stress and surface loads obtained in the registration process. The method was also tested with two abdominal CT scans of a pig acquired in flank and supine positions (Figure 3.25). We showed that while computing a physically admissible registration of the liver, the resulting constraint forces applied to the surface of the liver strongly correlate with the location of the anatomical boundary conditions (such as contact with bones and other organs) that are visually identified in the CT images.

Atlas-based Transfer of Boundary Conditions Biomechanical Simulation

In [START_REF] Plantefève | [END_REF], we introduced a method for automatic modeling of boundary conditions in deformable anatomical structures, which is an essential step in patient-specific biomechanical simulations. The method is based on a statistical atlas that gathers data defining the connective structures attached to the organ of interest. The atlas is registered with the patient's data using a physics-based technique, and the resulting boundary conditions are defined according to the mean position and variance available in the atlas (Figure 3.26). Preliminary parameter sensitivity analysis shows that the results are not affected by small changes of the boundary conditions elasticity. Two simulations are performed: (i) ground-truth simulation where the liver model is constrained using the boundary conditions obtained from manual segmentation and (ii) atlas-based simulation employing the boundary conditions estimated by our method. In both cases, a perfect match between the resulting shapes is obtained after the equilibrium is achieved. The method is then evaluated using abdominal scans of ten patients. The results show that the atlas provides sufficient information about the boundary conditions which can be reliably transferred to a specific patient.

Perspective

As stated previously, the non-rigid registration of biomechanical models still involves many open scientific questions. The methods introduced in this chapter bring solutions for specific medical applications, but they also have limits, and several aspects could be improved.

Numerical models and parameters: Concerning the modeling aspects, the parameters of the biomechanical are not patient-specific. Even if these parameters had limited importance in the context of imposed displacements (Wittek et al., 2009), it would be interesting to study their impact on the final registration procedure. In particular, the effect of heterogeneous elasticity for the tumor has to be evaluated. These parameters could be estimated intraoperatively using US elastography. While computing true static Young's modulus using shear waves is complex, the ratio between the tumor and surrounding soft tissue elasticity could at least be estimated.

2D/3D registration with fiducial markers: We are particularly interested in the robotization of percutaneous procedures in the liver, which is the purpose of the next chapter. However, none of the methods described above are adapted for this application. Indeed the liver is a large organ, and US images are too limited for the vessel-based registration. In addition, the marker-based registration method cannot be directly adapted since the abdominal cavity is not open during percutaneous interventions. Finally, the solution based on dynamic MRI images significantly raises the difficulty of developing a non-ferromagnetic robot. The imaging plane is represented along with the position of the radiating source (in blue), the internal marker's projections (black dots), and the resulting projection lines (red lines). To enforce registration, the 3D markers (green dots) are constrained on the projection lines using 2D projective constraints (blue and green arrows). The mechanical solver then finds the minimal energy position of the liver, enforcing the constraints.

We are working on a 2D/3D registration method using fluoroscopy and fiducials. We consider that a set of gold fiducials are manually placed percutaneously in the liver before the registration. A similar method has already been employed in Ohta et al. (2016) for radiotherapy of malignant liver tumors to guide the Cyberknife. Thanks to recent advances of Artis Zeego imagers, an initial volume image acquisition will be performed, providing 3D positions of fiducial at the initial step. During the automatic insertion, 2D projective positions of fiducials will be tracked in fluoroscopic images (Figure 3.27). During the breathing motion, fiducial markers are tracked with fluoroscopy, and the biomechanical model is registered in real-time with 2D projective constraints. Thanks to our collaboration with the IHU7 Strasbourg, we aim at testing our methods during an animal experiment.

Intra-operative confidence assessment of non-rigid registration based on biomechanical models: Another critical challenge concerns the validation of the registration procedure. Like other intraoperative compensation methods, we are sensitive to localization errors. However, validation is a significant obstacle to the actual usage of these numerical tools in the OR. Although the validation in a general context seems extremely difficult, we are developing a method to assess the uncertainty of the solution to provide clues of the reliability of the solution intraoperatively.

We are currently working on a method for online measurement of a confidence map associated with the non-rigid registration of a biomechanical model Baksic et al. (2019). The general workflow is illustrated in the figure 3.28. An augmented view combining the deformed model with medical images is displayed at high frequency, as done in this chapter. Asynchronously, at low frequency, image data and positions of the registered model are used as input of a mechanical study allowing evaluating the model's mobility. The confidence map can be used intraoperatively to only show the reliable parts of the numerical model in the augmented view. The mobility is defined as the maximal displacement of the model subject to both an estimated force identified to prevent injury of tissues and imagebased constraints used to perform the registration. The rationale behind this is to consider that image data provide ground truth positions positions of the organ, whereas other parts of the organ are subject to uncertainty mainly related unknown boundary conditions and mechanical parameters.

We acquired clinical ex-vivo data of a human's liver thanks to our collaboration with the Paul Brousse Hepatobiliary Hospital. An explanted liver was prepared with internal and external radiopaque markers and scanned in 4 deformed positions. The four shapes were then segmented from the CT scans along with the internal and external markers. Then, the marker-based registration method described above was applied to perform the registrations between shapes using the external markers. We computed the confidence map of the registered solution, and we found a strong correlation between registration errors and the mobility of the model. We also found that the confidence map is dependent on the marker's location, but interestingly, it also depends on the current shape of the organ (Figure 3.29). It can be explained by the mechanical coupling between constraints that are more or less important according to the stiffening of the deformation.

Conclusion

A new constraint-based formalism was introduced for non-rigid registration between a preoperative model and the intraoperative configuration. We applied the methodology on experimental data and fit as much as possible the clinical constraints to easily integrate the solutions in the operating room.

The methodology was first applied to compensate for the brain-shift phenomenon observed during tumor ablation procedures. A constraint-based simulation was processed to register the pre-and intraoperative vascular trees and the cortical surface with the probe footprint. Based on the displacement field of the model, preoperative data are then updated to provide the surgeon with navigable images corresponding to the current brain shape.

A semi-automatic registration of a pig liver has then been proposed. The method relies on the registration of 3D FE models to 2D dynamic MRI slices, providing a biomechanical extrapolation for the off-plane dynamic motion of organs. We showed that a priori knowledge of boundary conditions could be combined with intraoperative image-based data to provide a 3D extrapolation of the dynamic motion observed in a single 2D MRI slice.

Finally, a biomechanical registration method was proposed for liver surgery. We proposed to retrieve the 3D shape and position of the organ only given its 2D contour in a medical image and a set of 3D markers placed on the surface of the liver. The method is fast, robust, and compatible with intraoperative constraints but remains challenging to apply in a clinical environment.

Introduction

Robotic systems can assist percutaneous needle insertion in overcoming limitations due to human factors and increasing tool positioning accuracy. In the last decade, numerous solutions have been proposed for automatic needle steering. However, the deformation of structures (organs, needles) remains a significant obstacle for more task automation.

Traditionally the problem of the deformation is addressed by extracting a set of features from live images (also called visual servoing) and adjusting locally the pose/motion of the robot to compensate for the deformations Hutchinson et al. (1996). Nevertheless, visual servoing raises several limitations, in particular for needle insertion applications:

1. Per-operative images usually offer poor visibility of structures (such as tumors or vessels), and it is challenging to extract essential data in real-time. In addition, some information may be missing or not visible intraoperatively, such as the so-called disappearing liver metastases Robinson (2009). Due to chemotherapy effects, the appearance of liver tumors may change or become invisible in intraoperative images, even if the lesions still contain active tumors.

2. When large deformations occur, the control law of the robot can be significantly modified, which is extremely difficult to relate with imagebased displacements. For instance, when the needle is deeply inserted inside the tissue, the needle shaft becomes fully constrained, preventing any lateral motions of the needle.

3. Traditional controllers do not have access to any biomechanical or anatomical models capable of predicting the deformation of organs in real-time. However, the needle's path at the beginning of the insertion significantly impacts the possibility of reaching the target later. The needle trajectory must therefore be adapted as soon as possible.

Our goal is to insert a flexible needle in a deformable environment. A trajectory is manually defined based on a tomographic reconstruction (CT) in an undeformed configuration (Figure 4.1). A flexible needle is attached to the end effector of an articulated robot, considered as being infinitely rigid.

The soft structure is attached to a support, positioned within the working space of the robot. This work aims to provide Cartesian displacements of the robot's base such that the needle's tip remains on the predefined trajectory for any occurring deformation. It is important to note that an inverse problem must be solved to perform this insertion. Indeed, although the tra-jectory is completely defined at the initial step, it will deform as the insertion proceeds. Indeed, once the needle is inserted inside the volume, the needle becomes coupled with the tissue. Given that surgical needles are much stiffer (while still being deformable) than organs, any displacement of the robot's base will deform the structure and modify the planned trajectory. The input displacement commands must then be adapted accordingly. We underline that our goal is neither to provide an optimal trajectory nor to enforce the feasibility of the trajectory. Instead, we assume the trajectory being generated by a planning system taking into account mechanical effects such as it was proposed in Hamzé et al. (2016), or given by an expert (surgeon). If the input trajectory is not reachable, we expect the robotic system to stop falling at a local minimum.

Related works

Automatic needle insertion has received considerable interest over the past decades. Yet, recent surveys [START_REF] Abolhassani | [END_REF]; Cowan et al. (2011); Kaye et al. (2014); [START_REF] Elgezua | Survey on current state-of-the-art in needle insertion robots: Open challenges for application in real surgery[END_REF]; Siepel et al. (2021) still highlights the need for control models when dealing with deformable tissues. In this section, we review the main methods that have been proposed for accurate robotic needle steering.

Robotic needle steering

The advantages of surgical robots and manipulators are well recognized for daily assistance to telesurgery. Medical robots are developed for biopsies Seifabadi et al. (2012), prostate cancer Abdelaziz et al. (2011), breast therapy [START_REF] Zhang | [END_REF] and even neurosurgical applications Masamune et al. (1995); Ye et al. (2020). Typically, clinicians select a couple of images at the beginning of the intervention and manually define a trajectory and some landmarks to fix an entry point, a needle orientation, and an approximative path. The robotic system will then automatically or semi-automatically proceed towards the target following the set trajectory.

Robotic systems for percutaneous interventions are usually classified into two categories Cleary et al. (2006): table-mounted and patient-mounted systems.

While some patient-mounted robots partially compensate for physiological motions [START_REF] Piccin | [END_REF], the motions of internal organs are by far too complex and different from skin motions to be compensated for using purely passive solutions. Indeed, according to De [START_REF] Jong | [END_REF], the required tip positioning accuracy for liver tumor ablation is ≈ 3 mm. This level of accuracy is rarely obtained after the first needle insertion, and the practitioner usually needs to correct the trajectory several times, possibly stopping the patient breathing for a few seconds to limit disturbing motions [START_REF] Widmann | [END_REF]. Such iterations involving successive needle insertion and removal increase the risk of tissue damage and the probability of spreading tumor tissues along the needle path while removing the needle. Nevertheless, the complete procedure being generally too long for a single apnea, the deformation of both the organs and surgical needles remains an open problem that limits the development of automatic tasks performed by the robots in the operating room Kulkarni et al. (2019).

Needle and target tracking

Image-guided techniques extract information (for instance, needle and target positions) from vision sensors. This information is used in a closed control loop to guide the needle tip towards the target (called visual servoing).

A real-time system tracking is then necessary to take into account the deformations caused by external forces or natural motions (such as breathing) Neumann et al. (2013). Imaging modalities used for percutaneous procedures usually depend on the tumor visibility, the practitioner preference, and the local availability [START_REF] Lencioni | [END_REF]; Crocetti et al. (2010).

The most spread imaging modality is ultrasound (US) due to its low cost, harmlessness, and real-time capability. [START_REF] Kobayashi | [END_REF], track the needle shape using ultrasound (US) images. In Neubach and Shoham (2010) the authors estimate flexible needle's tip position from intraoperative US images. The method is combined with tissue stiffness estimation (from localized tissue displacements) and applied to the automatic needle insertion inside soft tissues. Okazawa et al. (2006) presented two methods to detect the needle in 2D ultrasound that specifically addresses needle curvature. Vasconcelos et al. ( 2016) introduced a calibration method for both 2D and 3D ultrasound probes that involves scanning an arbitrary region of a tracked needle in different poses. Although ultrasound is fast, portable, widely available, and easy to combine with robotic systems, the needle visibility in US images remains complicated. In addition, the US is limited to superficial tissues and subject to artifacts due to air-filled structures such as lungs or bowel Pua and Sofocleous (2010).

Other image modalities have been investigated for needle tracking (see [START_REF] Abolhassani | [END_REF] for a survey). Real-time imaging such as CT fluoroscopy or real-time MRI is sometimes chosen for complicated cases Hawkins et al. (2016); Jiao et al. (2018). For instance, [START_REF] Navab | [END_REF] uses X-ray fluoroscopy to 3D align a needle, held by a medical robot, inside a porcine kidney. Seifabadi et al. (2013) presents an MRI-compatible robot for teleoperated bevel-tip needle steering under real-time MRI guidance. The image quality of MRI is often much better, but the acquisition time is usually slower, and it raises several difficulties to align the images' plane with the structures. In addition, due to the difficult access to the patient in MRI bores and overexposure to X-rays in CT fluoroscopy, the most common alternative for percutaneous procedures into the liver remains conventional CT Puijk et al. (2018). In this case, insertions are done outside of the imaging device, and 3D images are only taken at critical moments of the operation but not in real-time. As a result, it raises significant difficulties for needle guidance during the whole procedure.

Alternatively to medical imaging systems [START_REF] Abayazid | [END_REF] used an optical fiber embedded into the needle's shaft for direct measurement of the deflection and even for three-dimensional reconstruction of the needle shape. [START_REF] Kim | [END_REF][START_REF] Kim | [END_REF] use a set of Fiber Bragg Gratings (FBG) sensors and general elastic rod theory to reconstruct the shape of the needle in 3D even for significant deflection of the needle.

The solutions mentioned above share two main limitations. First, they all rely directly on images to track the needle and define the target point, either using tomographic reconstruction (CT, MRI) or live images (echography, fluoroscopy), entailing many challenges for image features extraction and noise correction. In addition, image quality is most of the time antagonistic with the acquisition frequency. It limits applications either to off-line insertions or raises significant difficulties on image-processing algorithms and images localization. Augmented Reality is a promising technique for mini-invasive surgery assistance, and robotic-assisted procedures Fischer et al. (2007) but most augmented reality systems used for percutaneous procedures are still limited to rigid registration. The second significant limitation is related to the fact that the correction of the needle trajectory is always performed a posteriori (i.e., only when errors are detected in the images or sensors). None of these methods can predict the behavior of the tissues to generate a priori deformations to reach the target, which is essential for needle procedures since the path taken by the needle inside tissues depends on all the history of the insertion and may prevent it from reaching the target.

Needle steering and control

The most common strategy to control a needle with a robotic system is to manipulate the needle base position to modify the needle tip path inside tissues Kaye et al. (2014). In this case, the needle base's six degrees of freedom (DOF) are controlled by the robot end-effector where the needle is attached. Relating the 6 DOF to the needle tip motion is a challenging inverse problem, especially if this problem has to be solved in real-time.

The method named duty cycling consists of spinning the needle along its insertion axis to add three additional Degrees of Freedom (DoFs) at the tip of the needle Rucker et al. (2013), hence allowing the steering of the needle inside the volume during the insertion. [START_REF] Bernardes | [END_REF] developed a robot-assisted approach for the automatic steering of flexible beveled needles in percutaneous procedures. The method uses a duty-cycled rotation of the needle to perform insertion with arcs of adjustable curvature. The method used closed-loop imaging feedback with intraoperative motion replanning strategies to compensate for system uncertainties and disturbances. Krupa (2014) presented a duty-cycling robotized system for steering beveled needles allowing this way for the creation of complex non-straight trajectories to reach a target and avoid obstacles. Secoli et al. (2016) proposed a bio-inspired multi-part needle and validated the control strategy by fitting experimental models. [START_REF] Reed | [END_REF] described a robot-assisted needle steering system using three integrated controllers and proposed a stochastic motion planner with paths replanning for duty-cycling techniques.

Instead of steering the needle to the tumor, another approach proposed by [START_REF] Mallapragada | [END_REF] is to move the tumor toward the needle trajectory. The method is used for breast biopsies and takes as input real-time fluoroscopic images in which the tumor is located.

Kinematic and mechanical model-based control

One of the first models for needle-deformable tissue interactions was proposed by [START_REF] Dimaio | [END_REF]. Based on a linear model, the method allows simulating stick-slip behaviors. A significant amount of work was also performed by Misra et al. (2010); Abayazid et al. (2016), including experimental validations with robotic systems to steer bevel-tip needles.

DiMaio and Salcudean Dimaio and Salcudean (2005) were among the pioneers to investigate robotic needles' steering through soft tissue. They compute the Jacobian matrix numerically from mechanical models, from which the needle base velocity is derived and used as input in an open robotic control loop. The method is combined with trajectory planning strategies: attractive fields drive the needle towards the desired target, whereas repulsive fields avoid obstacles, but the method is not real-time. [START_REF] Glozman | [END_REF] proposed a real-time steering system that integrates planning and control in a closed control loop for dynamic systems. A mechanical model (springs with different stiffness coefficients along the needle shaft) simulated the interaction between the needle and soft tissues.

Modeling the interactions between needle and tissue is essential to predict the behavior of the needle inside tissues. Misra et al. (2008) studied interaction force at the tip of asymmetry bevel-tip needles to estimate the rupture toughness of the tissues. Rucker et al. (2013), proposed a closed-loop control system for asymmetric-tipped needles. [START_REF] Abayazid | [END_REF] proposed an image-guided control system to steer flexible needles. Both kinematics and mechanical models are used to predict the deflection of the needle model.

Advanced mechanical models have been used to predict the behavior of the tissue. Barbé et al. (2007a) uses a linear Kelvin-Voigt (KV) model to estimate online the forces involved in percutaneous interventions, but the method is limited to the forces applied in the direction of the needle. Khadem et al. (2015) presented a mechanics-based model for the simulation of needle insertion in soft tissues. The proposed model is based on beams' theory, allowing the prediction of the deflection of the needle. Robotic experiments are conducted to identify the model's parameters, then used to steer the actual needle. [START_REF] Kobayashi | [END_REF] uses a nonlinear viscoelastic model calculated intra-operatively to manipulate a needle while considering organ deformations. The force upon the needle is measured using a force sensor and applied to a mechanical model to estimate tissue deformations. Nevertheless, no interaction models between the needle and the tissue are proposed restricting the method to superficial insertions. The method is then extended in [START_REF] Kobayashi | [END_REF] using an ultrasound-guided manipulator combined with a physics-based model of the liver. The biomechanical registration provides internal stresses during the insertion, but this method is limited to 2D insertions. Recently [START_REF] Chiroiu | [END_REF] showed that considering deformations of the needle and the liver results in a better needle navigation mechanism. They use the Cosserat elasticity to describe the interaction be-tween the needle and the human liver.

Numerical simulation and trajectory planning

Large-scale simulations of needle insertion (including advanced FE models of needle, tissue, and interactions) have also received considerable interest for training purposes. Chentanez et al. (2009) presented a FE approach based on the beam's theory to predict needle deflection, but the method relies on expensive re-meshing steps along the needle shaft. Interactive models were proposed in [START_REF] Duriez | [END_REF]. The method does not require any re-meshing process when the needle goes through the tissue. In addition, it allows for the simulation of complex phenomena such as tissue deformations, needle-tissue friction, and puncture force.

Beyond modeling, several researchers have developed motion planners for flexible needles in 3D tissues Park et al. (2008); [START_REF] Jijie | Motion planning for steerable needles in 3D environments with obstacles using Rapidly-exploring random trees and backchaining[END_REF]. Duindam et al. (2010) derived an inverse kinematics solution to reach the desired position and orientation in 3D. However, all these methods assume that the tissue is rigid. Alterovitz et al. (2008) presented trajectory planning algorithms, including probabilistic methods considering uncertainty. The method was improved in [START_REF] Wen | [END_REF] to explicitly consider motion and uncertainties while guiding the needle to a target in 3D anatomy. Hamzé et al. (2016) relied on the method proposed in [START_REF] Duriez | [END_REF] for optimal trajectory planning for liver surgery. The method considers deformations (breathing, needle deflection, friction) to avoid obstacles identified at the planning step. The main limitation of off-line trajectory planning strategies is that significant changes might occur between preoperative and intraoperative configurations and invalidate the chosen trajectory in an actual application. Li et al. (2017) proposed a path planning (and online replanning) approach for steerable needles based on discrete potential fields in 3D anatomical structures, but the method is not real-time.

Advanced FE simulation and needle/tissue interaction models can predict the behavior of needles inside tissues, but such complex models have not been used to control robotic systems. One of the main reasons is that once the needle is inserted, the overall system (robot, needle, and tissue) can be seen as a deformable robot with infinite degrees of freedom, which cannot be controlled with standard robotic approaches. The control of soft robots is a recent research topic. Largilliere et al. (2015) proposed an advanced control strategy based on FE models to deform a soft robot. An inverse problem based on a QP (quadratic programming) algorithm is used to solve the equations of motion and control a soft robot in real-time. However, interactions between deformable structures with complex behaviors, such as nonlinear friction, are not considered. Similarly, Coevoet et al. (2019) solves equations of motion and control soft robots in real-time. However, the optimization problem is defined using the actuators' space of the considered robot, which cannot be directly applied to needle insertion where the deformation of the tissue is entirely independent of the robotic system.

Telemanipulation and shared control

Telemanipulation systems can filter the operator's tremors, scale clinician's movements to enhance accuracy, and reduce X-ray exposure for the medical staff. Robots dedicated to percutaneous interventions can be classified into three different types, depending on their level of autonomy: manual, automatic, and with shared control:

Manual: Robot for the assistance of interventional radiology are generally operated manually Kulkarni et al. (2019) A recent example of a teleoperated system starting clinical trials is the Zer-oBot developed by Hiraki et al. (2020). Another example is the DEMCON needle placement system Arnolli et al. (2018). The rotation of the needle guide is done automatically using CT imaging following a predefined straight path. Such manual systems may solve the problem of initial registration and planning Heerink et al. (2019), but do not take into account the needle-tissue interaction.

Automatic: Automatic needle steering methods can compensate for such complex interactions (see Li et al. (2018) for a recent review) and therefore received significant interest in the last two decades. Methods have been proposed for both symmetric [START_REF] Alterovitz | [END_REF]; DiMaio and Salcudean (2003); [START_REF] Glozman | [END_REF]; Kojcev et al. (2016) and beveled-tip needles Alterovitz et al. (2008); Chevrie et al. (2016); Lapouge et al. (2021). Chevrie et al. (2016) rely on a needle/tissue interaction model combined with base manipulation and tip deflection to steer a beveled-tip needle. Though models in these methods may increase the tip position control's predictive nature, they require associating observed features with the underlying models. It makes these methods highly sensitive to registration errors, which are extremely difficult to compensate for per-operatively.

Shared Control: Shared control methods enable the user to be proactive during the task and thus more responsive to errors Abbink et al. (2012). However, due to the difficulty of merging user and automatic inputs, only a few shared-control systems have been proposed for the needle steering inside soft tissue. Several authors Abayazid et al. (2016); Wartenberg et al. (2016); Chevrie et al. (2019) propose methods to steer beveled-tip needles. Chevrie et al. (2019) gives the user complete control over the needle tip velocity through a haptic interface. The user feeds the desired tip velocity to the automatic needle steering method proposed in Chevrie et al. (2016) that is used to control the needle. Wartenberg et al. (2016) only gives the possibility to the user to act on the insertion velocity, while the bevel orientation is obtained automatically.

Methodology

Imaging System Simulation Robot In this section, we introduce the closed-control robotic loop for a complete automatic insertion (Figure 4.2). As stated in the introduction, we propose to rely on inverse FE simulations to predict the behavior of deformable structures. A critical difficulty concerns the computation time of inverse steps. Indeed, the behavior of the overall system (robot, needle, and tissue) being highly non-linear, the validity domain of the simulation is limited to small deformations and remains valid only for a small amount of time. To overcome these limitations, we introduce a numerical method that allows solving inverse Finite Element (iFE) simulations at high frequency. Our approach relies on a forward FE simulation of a needle insertion (involving complex non-linear phenomena such as friction, puncture, and needle constraints). Control commands are then derived from two critical steps:

Corrective Step: As for visual servoing, we extract a set of features from live images to enforce the consistency of the models with actual data. However, instead of directly steering the needle toward these features, we first register FE models with the observations. The advantage of relying on FE models is that it provides a regularization technique to extrapolate the displacement field of extracted data.

Predictive

Step: Input commands of the robot are obtained from an optimization process based on iFE. It allows anticipating the behavior of mechanical structures to adapt commands much faster than waiting for a correction from the images. Inverse steps are performed to numerically derive the socalled Jacobian of the Simulation, which relates Cartesian displacements of the base of the needle with displacements of the tip inside the volume, allowing this way to compensate, or even induce, necessary deformations to reach a target.

Needle insertion model

We first introduce the numerical model of the needle and the interaction model between the needle and volume elements.

Model of the needle

The model of the needle is based on the beam theory Petyt (1969). The needle comprises beam elements, i.e., segments of 2 nodes with 6 DOF per node (translations and rotations). The local stiffness matrix Kn e for the first beam element of the needle is computed as follows:

Kn e =                     a 0 b zy 0 0 b yz 0 0 0 f sym 0 0 -c yz 0 d zy 0 c zy 0 0 0 d yz -a 0 0 0 0 0 a 0 -b zy 0 0 0 -c zy 0 b zy 0 0 -b yz 0 c zy 0 0 0 b yz 0 0 0 -f 0 0 0 0 0 f 0 0 -c yz 0 e zy 0 0 0 cyz 0 d zy 0 c zy 0 0 0 e yz 0 -c zy 0 0 0 d yz                     (4.1)
with a = EA l , b zy = 12EIz l 3 (l+Φy) , b yz = 12EIy l 3 (l+Φz) , c zy = 6EIz l 2 (l+Φy) , c yz = 6EIy l 2 (l+Φz) , d zy = (4+Φz)EIy l(l+Φz) , d yz = (4+Φy)EIz l(l+Φy) , e zy = (2-Φz)EIy l(l+Φz) , e yz = (2-Φy)EIz l(l+Φy) and f = GJ l . The above coefficients are defined with the following geometrical and mechanical parameters: A is the cross-sectional area of the beam, and l is the length. G = E 2 * (1+ν) where E is the Young modulus of the beam, and ν is the Poisson ratio. J = Iy + Iz is the Polar moment of inertia, where Iz and Iy are the cross-section moment of inertia about the respectively z and y axis of the beam. Assuming mass ratio = 1, Ix = Iz = π * (r 4 -r 4 )

4

, where r and r are respectively the external and internal radius of the beam. Finally, Φ y = 12EIz

GSyl 2 and Φ z = 12EIy GSzl 2 with S y and S z being the shear area in the y and z directions.

The above matrix is computed in a local reference frame oriented along with the frame of the first node of the beam. Therefore transformation matrices are necessary to change the reference for the remaining beam elements composing the needle. The local 12×12 stiffness matrix K n e is given by:

K n e = Λ Kn e Λ T (4.2)
where Λ is a matrix obtained from the direction cosines of angles between the local and global coordinate systems. In addition, the model is combined with the corotational formulation allowing this way for the simulation of large rotations of the beam while being restricted to small elongations. Deformable models composed of beams produce a particular numerical system. Indeed, if the points are organized in ascending order (Figure 4.3) along the beam, the matrix generated by the system is expressed in a Block Tri-Diagonal (BTD) structure. It creates three non-zero data blocks on each row of the matrix since the displacement of a point only influences its two neighboring points. Although the values change at each time step, the inversion of this system is very fast using the algorithm of Thomas Kumar et al. (1993), which has linear complexity in the number of nodes.

Needle mechanical characterization

The needle is parameterized with geometrical and mechanical parameters that must be identified. We used a 21 Gauge, 12 cm Sterican ® needle (B. BRAN Melsungen AG). The Poisson ratio ν = 0.3 is chosen equally to steel Dehghan et al. (2006).

While the geometrical parameters can be measured directly, the Young modulus and the discretization of the beam still have to be evaluated. We evaluate these parameters using image correlation (see below for the experimental setup), allowing optimizing the models by fitting the parameters with the recorded measurements (Figure 4.4). The needle is positioned horizontally (thanks to the robot) with various loads attached to its tip. Several simulations were performed varying the value of E and the number of beam elements discretizing the model. We measured the back projection (in pixels) between the simulated model and the actual bent needle. We found experimentally E = 200GP a, which is consistent with the literature Khadem et al. (2015), and a converged solution with 28 elements discretizing the shaft.

Needle-tissue interaction model

When the needle is inserted or in contact with the volume, the two models become coupled, and the mechanical actions of one model need to be transferred to the other one. The computation of these forces and displacements requires common DoFs. Lagrange multipliers λ are used to impose constraints, as proposed by [START_REF] Duriez | [END_REF]. According to the positions p n of the needle and p v of the volume, the following constraints are used (Figure 4.5):

Penetration constraint H ϕ (p n , p v ) is applied before penetrating the tissue, between the needle's tip and its closest point on the surface. The surface is defined as the external subset of the triangles of the tetrahedral or hexahedral volume mesh. H ϕ is a 3Dofs constraint: a unilateral contact force is applied along the normal of the triangular surface to avoid the penetration, whereas Coulomb friction is added in the tangential direction. Depending on the situation the constraint has one of the following behaviors (Figure 4.6):

1. When the needle is not in contact with the volume, the penetration constraint is applied between the needle's tip and the respective closest point on the surface triangle. H ϕ satisfies the Signorini conditions λ ⊥ δ, i.e., while the needle tip is located on the positive side of the constraint, no contact force is applied. However, considering the constraint in the mechanical system at this step enables applying a reaction force as soon as the needle touches the surface. This is particularly important for both haptics that runs asynchronously and may reach the contact before the simulation and our robotic application to predict the contact and the behavior of the needle/tissue according to the robot's motion.

2. When the needle is in contact, a reaction force is applied to cancel the violation (green arrow in the figure 4.5a). In addition, friction (parameterized with µ s ) is added in the tangential plane (red arrow in the figure 4.5a). Objects are considered in contact until λ < p f , leading to a deformation of both the needle and the volume.

3. When the force exceeds the threshold p f , the contact stress is transformed into a trajectory constraint allowing the tip of the needle to cross the membrane of the tissue.

Needle constraints H ψ (p n , p v ) Once the needle has crossed the surface, a set of trajectory constraints is used to enforce the shaft of the needle to follow the path created by the needle's tip. The constraints are created dynamically as the needle moves through the volume. The distance d n controls the discretization between constraints, i.e., a new constraint is created, at the needle's tip location, as soon as the tip is at least at a distance d n of all the other trajectory constraints. Each constraint is defined by its barycentric coordinates into the tetrahedral mesh, allowing their definition at an arbitrary location within the volume without needing expensive re-meshing operations. Moreover, when the volume is deformed, the trajectory is recomputed accordingly using the barycentric coordinates.

Needle constraints are created along the path taken by the needle inside the volume. Each trajectory constraint H ψ is composed of a set of 3DoFs constraints applied in the direction of the following vectors: the first vector n i 0 (red arrow in the figure 4.5b) aims in the direction of the following trajectory constraint (or to the needle tip for the last point) and is used to apply resistance against the penetration along the needle's shaft. Two additional vectors n t 1 , n t 2 (blue arrows in the figure 4.5b) are defined in the orthogonal plane in order to prevent any relative displacements of the needle in the tangential plane. These vectors define the constraint directions as explained in the previous chapters. When the needle is inserted into the tissues, Okamura (2004); Misra et al. (2008) showed that resistance force acts on the tip of the needle. Therefore, the interaction model has an additional parameter, 0 ≤ µ n ≤ 1, (0 no friction, 1 sticking), the penetration resistance coefficient along the shaft. If µ n ̸ = 1, the needle can penetrate and move along the trajectory. In this case, the trajectory constraint H ψ is applied as follows:

1. The constraints H ψ are associated to the respective closest point on the needle, defining this way the material point on which constraint will be applied on the needle's model.

2. The constraint are solved as bilateral constraints (both in the directions n i 0 and n t 1 , n t 2 ), resulting in respective constraint forces λ i 0 , λ t 1 , λ t 2 preventing any relative displacement in the 3D space.

3. Finally, the constraints force λi 0 = λ i 0 × µ n × 1 dn is scaled with both the friction coefficient µ n and the inverse of the constraint distance in order to get independent reaction force from the discretization of the constraints.

The proposed model provided a constraint-based formulation of needle insertion simulations. It is particularly relevant for the robotic control strategy proposed in the next section. The constraint-based formulation is also essential for haptic applications since the user can feel the resistance along the shaft and the mechanical coupling between the needle and the volume.

Registration and experimental setup

As stated above, to enforce the system's consistency, the models need to be registered with image data. The registration is only performed for the tissue and the robot. Indeed, the needle is considered rigidly attached to the end effector of the robot. Therefore, our method's essential advantage is that the needle's positions are derived from the mechanical model and interaction constraints with the volume without explicitly tracking the needle in images.

The non-rigid registration of the volume procedure is performed with the marker-based method introduced in the previous chapter [START_REF] Adagolodjo | [END_REF]. This choice is motivated by the fact that the goal of this chapter is only to evaluate the control strategy based on iFE simulations, and for this purpose, it is necessary to reduce registration errors at a minimum. Therefore, although the method is not directly applicable to percutaneous procedures (since it is based on markers), it is currently the most stable and accurate solution for fast non-rigid registration with large and dynamic deformations.

Experimental Setup

To test our approach, we developed the following experimental setup (Figure 4.8). It includes: i) optical tracking system; ii) monocular camera iii) robotic arm with its needle holder; iv) foam and its support. (i) -Tracking system: We used an OptiTrack motion-capture system1 which includes both hardware and software components for the calibration and localization of cameras. The system is composed of 6 Flex13 cameras arranged around the working space of the robot. The tracking system provides the 3D positions (at high-speed sampling, i.e., at least 120 frames per second) of a set of markers placed at the surface of the foam. After the calibration, the system reports back-projection errors of 0.025 mm on average. (iii) -Foam: We chose a polyurethane foam for its high resistance against fracture in case of lateral motion of the needle and its low friction resistance against the penetration of the needle. The foam is attached and fixed at its four corners within the working space of the robot. The central part is free and can be deformed both tangentially and laterally to create complex 3D deformations. The foam dimensions are voluntarily long and thin [12 × 6 × 1] cm to create potentially large deformations and minimize errors during the non-rigid registration step. Nevertheless, it also raises significant difficulties to maintain needle tip within the volume of 1cm thick, including during large in-depth deformations. (iv) -Robot: The 6DOF Mitsubishi RV1A robot is a 6-arm anthropomorphic robot. The kinematic model R(q)=X provides the position of the end effector X according to each joint q. The Denavit-Hartenberg formulation is given in the figure 4.10. The robot includes a 3D printed needle support mounted on its end-effector. The CAD model was designed to include the location of 6 markers allowing the definition of a rigid body attached to the support (Figure 4.9). These markers allow solving the Hand and Eye Problem, transforming the robot's frame coordinates to the tracking system coordinates. After the registration with the cameras, the system's accuracy is evaluated, solving a kinematic problem to position the tip of the needle (assumed rigid) at several 3D positions given by markers. The needle base is known from the CAD model, whereas the tip position is evaluated with the beam model. For the rigid positioning of the tip, we reported a back-projection error in the monocular camera's views of less than 1.0 pixels between the needle's tip and the markers, showing consistent registration of the overall system.

Registration and Constraints

The registration of the system is performed using Lagrange multipliers that are associated with the following constraints:

Bilateral constraints H χ (p n , X ) are used to fix the position of the needle's base to the terminal part of the robot (Figure 4.11a). H χ is a holonomic 6DoFs constraint (position and orientation) whose violation is defined as the relative displacement between the needle's base and X the position of the end effector located on the terminal part of the robot. Since the robot is infinitely rigid, only the needle is affected by H χ .

Observation constraint H Ω (p v , m) H Ω is a set of 3DoFs bilateral constraints, used to register the model of the foam with respect to the observations (see Fig.4.11b). These constraints are applied between observation points m given by an external tracking system and a set of 3D points on the model's surface m, segmented during the meshing step. An Iterative Closest Point (ICP) algorithm is used to bind each point of m with its respective closest point in m. The figure 4.12 shows the accuracy of our approach using 25 markers uniformly distributed at the surface of the foam. We measured the Hausdorff distance between the projected contour of the model and the outline of the foam manually segmented in images. A mean error varying between 1 and 3 pixels is reported with maximum values of 6 pixels.

Contrary to the needle, volume model parameters may not be accurately known during the insertion (boundary conditions, attach points, mechanical parameters. . . ). The Young Modulus of the foam was evaluated at E = 1.3kP a Zieliński and Witkiewicz (2006). However, constraints embedded in H Ω enforce a small displacement error of volume's position p v in the neighborhood of the observations m; which significantly decreases the sensibility of the method with the mechanical parameters of the volume. A sensitivity study was performed varying E ± 20% (which corresponds to the standard uncertainty of liver parameters in healthy subjects Roulot et al. (2008)), without a significant impact on the registration accuracy.

Objective function

The numerical simulation of a needle insertion described above can be written as a forward non-linear problem S (X , p n , p v , m), where X is the position of v of the needle and the volume that fulfills constraints H χ , H ϕ , H ψ and H Ω , for any position of the robot X (t) and observation m (t) .

The trajectory is defined as a set of connected points s whose positions are given by a linear relation s = J T p v , where J is derived from barycentric coordinates of the trajectory at the initial step. Although the trajectory is virtual and cannot be tracked by any imaging system, it is possible to estimate its position for any occurring deformation using the position p v of the registered mechanical model of the foam. Let n (t) be the position of the needle's tip after integration, and t (t) be the desired point on the trajectory. The target point on the trajectory t (t) is defined by a parameter c ∈ [0..1] allowing moving the target from the first to the last point of the trajectory. The needle must satisfy the following conditions to advance the target point c on the trajectory (and perform the insertion): First the needle's tip must be located on the target point (i.e. ||n (t+1) -t (t+1) || = 0). In addition, to anticipate the displacement of the needle inside the volume, the needle must be oriented in the same direction as the target trajectory (Figure 4.13). Formally, the objective function I(X , p n , p v , m) is given by the following equation:

I(X , p n , p v , m) = n -t η θ = 0 (4.3)
where θ = arccos(dot(⃗ o, ⃗ t)) is the angle between the needle's tip direction and the tangent of the target trajectory and η the gain to weight the constraint function.

An inverse problem must be solved to compute the next position of the robot X (t+1) , minimizing the objective function I.

Contributions

The critical contribution of this section is to develop a numerical framework to solve inverse Finite Elements (iFE) simulations involving complex needletissue interaction constraints to derive robotic commands. The non-rigid registration process allows maintaining low deformation errors between FE models and actual structures. Then, we extend the method with a shared control strategy to increase the safety, stability, and accuracy as well as the acceptance of the developed solution (Figure 4.14). The primary motivation is to leave potentially dangerous decisions and actions to the practitioner, whereas complex non-intuitive manipulations of the needle are performed automatically, in particular, to compensate for breathing motions. 2016) we proposed to compute the Jacobian of the simulation J to steer the needle and reach the objective function I. J is a 6 × 4 matrix which is numerically derived applying successive Cartesian perturbations δX i to the end effector X , and estimating the resulting objective function in the simulation:

J[: i]= I(X (t) , p (t)
n , p

v , m (t) )-I(X (t) +δX i , p

n , p

(t) v , m (t) ) ∥ δX i ∥ (4.4)
where J[: i] is the column i of the Jacobian and I(X (t) , p

n , p

v , m (t) ) is the value of the objective function at the beginning of the simulation step.

The displacement of the robot minimizing the objective function I is obtained with a Newton method. For this purpose, a pseudo-inverse J * (using Singular Value Decomposition) is used to solve the inverse problem. However, the above equation requires solving 6 independent forward FE simulations S for each input robotic command. Therefore, the computation of the Jacobian is the most critical step in terms of computation time and must be carefully optimized to maintain adequate insertion time. Increment Target: c+ = δc 8: end if 9: for i = 0 to 6 do 10:

Algorithm 1 Inverse Simulation Loop

Compute Violation: δ i = H(p, X + δX i , m (t) ) 11:

Solve Constraints: Wλ i = δ i 12:

Corrective Motion: p i = x free -A -1 H T λ i 13:

Compute Jacobian: J i = I( X ,p)-I( X +δX i ,p i ) ∥δX i ∥

14:

Reload data: p = p 15: end for 16: dX = J -1 • e 17: Move End Effector: X (t+1) = X + dX 18: Compute Violation: δ = H(p, X , m (t) ) 19: Solve Constraints: Wλ = δ 20: Corrective Motion: p = x free -A -1 H T λ

The essential contribution was to show that the Delasus operator W can be computed only once per simulation step to compute the columns of the Jacobian. Instead of solving the complete FE simulation for each column of J, they are obtained by only solving a local constraint problem for each perturbation δX . The Jacobian is therefore evaluated in the constraint's space which is much smaller than the motion's space. This way, its computation only involves the lines 9 to 15 of the algorithm 1:

1. Compute a new violation δ i for perturbations δX i (line 9).

2. Solve constraints providing Lagrangian's multipliers λ i (line 10) . We evaluated the accuracy of the control approach with synthetic data. We compared two strategies where the robot is controlled using: i) our inverse simulation ii) a naive approach where the needle and the trajectory are considered constant and rigid during all the insertion. Although the naive approach could lead to a successful insertion over the first centimeters (when the trajectory stays straight), significant errors are progressively introduced due to a miss evaluation of the deflection of the needle in the control loop. Instead, our solution automatically followed complex trajectories, including during large deformations of the models.

4.4.1.2 Robotic insertion of a flexible needle in deformable structures using inverse Finite Element simulation In Adagolodjo et al. (2018), the method was extended to an experimental evaluation. In order to maintain reasonable insertion time, we proposed an asynchronous control system (Figure 4.16). The external imaging system provides observations m located on the surface of the foam at high frequency. A dedicated thread computes at high frequency a joint-based interpolated motion between the current position of the robot and the desired position of the end effector X (t+1) provided by

S (t) = S(X (t) , p (t) n , p (t) v , m t ); I (t) = I(X (t) , p (t) n , p (t) v , m t ).
the simulation. The inverse simulation steps are executed asynchronously to compute the Jacobian J (t) . Based on the Jacobian J (t) and the error I, a new desired position of the end effector is computed and sent to the robot's thread asynchronously (i.e., before X (t+1) is reached).

Evaluation of the accuracy during needle insertion: We developed the experimental setup introduced in the previous section and evaluated the method in a real scenario. To define the trajectory and enforce its feasibility, we manually inserted the needle where a metallic thread was beforehand slipped within the shaft of the needle. During the manual insertion, deformations were applied on both the needle and the foam, which created a curved path. After the insertion, the needle was removed, letting the metallic wire within the volume. A CT scan including the foam, the metallic thread, and markers was performed and segmented to create FE meshes. A robotic insertion was then performed to follow the desired path automatically. During the insertion, vertical and lateral deformations were applied (see the figure 4.17), leading to significant modification of the undeformed trajectory and necessary bending of the needle. Despite deformations, the automatic method maintained the tip of the needle within the thickness of 1cm the volume and followed the desired path without any human intervention. Before the insertion, another metallic thread was placed within the shaft of the needle, allowing for error measurement between manual and robotic paths (see Fig. 4.18). We reported an average error along the trajectory of 1.62 mm with a maximum error of 3.73 mm which is acceptable for a medical application. Finally, the total insertion time has been performed in 6 minutes, which is similar to what is obtained by other methods Schulz et al. (2013), allowing this way for clinical applications. 4.4.1.3 Robotic needle insertion in moving soft tissues using constraint-based inverse Finite Element simulation Although the above method provides a significant speedup to compute iFE steps, the method was limited to quasi-static insertions where deformations only result from the needle interactions with the volume and displacements of the robot. In Baksic et al. (2020), the problem of needle insertion during respiratory motions is explicitly addressed.

The above solution requires evaluating the objective functions in the motion's space (i.e., line 12 of algorithm 1). Although this operation is usually not the bottleneck for traditional FE simulations, it raises computation time issues when performed several times in a close robotic loop. Therefore, we proposed to formulate the objective functions as constraints objectives functions. These additional constraints do not have a physical meaning, but instead, they allow to compute the variation of the objective functions directly in the constraint's space according to the needle's base perturbations.

The W matrix is augmented as follows:

         w 1,1 . . .
                  λ 1 . . . λ k 0 . . . 0          =          δ 1 . . . δ k ē1 . . . ēp          (4.5)
Lines 1 to k correspond to the original Compliance matrix used to impose constraints in the inverse simulations, whereas lines k + 1 to k + p correspond to the objective constraint functions. The mechanical influence of the constraint objective functions on all the other constraints is removed by setting columns k + 1 to k + p to 0. It enforces that mechanical constraints are not impacted by the constraint objective functions. Instead, the mechanical influence of all the other constraints is conserved (lines k + 1 to k + p), allowing retrieving the value of the objective functions without any need for re-projection in motion's space.

The time saved with this approach enabled combining the approach with stabilization strategies which are particularly relevant for needle insertion applications in a dynamic environment (i.e., breathing). We proposed to compute the centered Jacobian of the simulation, requiring the computation of 6 additional inverse simulation steps (12 steps). In addition, we also proposed to use a Tikhonov regularization while inverting the Jacobian.

The method was evaluated using a synthetic simulation of a pig liver during respiratory motion. The models have been generated from in-vivo data, but the motion of the liver is generated by a direct FE simulation involving a complex interaction between the needle and the surrounding environment. Several methods were compared to automatically steer the needle along the trajectory (Figure 4 J). In addition, the current strategy based on constraint objectives brings an additional speedup of more than 2×. In a context where the liver moves independently from the robot and the needle, the proposed solution was the only one to enforce the system's stability. Indeed, even if the solutions are numerically identical, the delay introduced by other methods adds additional errors, and the methods cannot compensate for the deformations induced by the breathing motion.

4.4.2 Shared control 4.4.2.1 Shared control strategy for needle insertion into deformable tissue using inverse Finite Element simulation For responsibility and acceptance reasons, it is necessary to let the practitioner control the needle insertion, at least in combination with the robotic system. In [START_REF] Baksic | [END_REF] we proposed a shared control strategy (Figure 4.20). The method is composed of two parts: the needle steering and the decision-making. The steering part is done by the automatic method designed to follow inputs provided by the user. Contrary to usual teleoperated systems, the user controls the needle tip position relative to the tissue without considering complex motion such as breathing. In addition, haptic guidance is added to help the user follow a planned path without any need for compensation for the global motion of the system. A direct simulation involving multiple deformable organs and contacts was implemented. This simulation represents the physical system shown in figure 4.20 in which the user has to place the needle. Based on this simulation, virtual X-ray images are generated and displayed to the users. The objective is to compare both the needle tip positioning accuracy and the insertion time in different cases. The method was evaluated by several users performing a simulated needle insertion to reach a tumor located at a depth of 60 mm within the liver (100 mm deep under the patient's skin). The shared-control method introduced in this article (SC) is compared with a fully teleoperated insertion (Man) on the one hand, and with a fully automatic one (Auto, described above) on the other hand. A total of six untrained users and the system developer, denoted as the highly trained user, did the experiments. The users were asked to place the needle tip tumor's center without any limit of time. Five successive insertions of each type were performed, where users had as much time as needed to familiarize themselves with the interface. The results show that the untrained users divided by five the absolute placement error using any SC method compared to Man.

Advanced modeling

Among the difficulty to accurately place needles, additional constraints lie in the fact that needles may lacerate and damage surrounding tissues due to the physiological motions of organs during the insertion. Although the lacerations are limited when the needles are held by a practitioner (avoiding applying too much stress on the needle inside tissues), they are mainly observed when needles are held by a stiff, robotic system Bayle et al. (2014), which is why some robotic devices were conceived to avoid this behavior [START_REF] Piccin | [END_REF]. In this case, lacerations significantly impact the needle tip motion inside the tissue. Cutting models are therefore needed to produce relevant training simulators. The proposed method relies on an intermediate plan defined independently from the needle and tissue mesh resolutions. The rupture and extension of the cut are triggered by the interaction forces between the needle and the tissue.

The cut plane is used to simulate the additional degrees of freedom of the needle inside the tissue, resulting in an interactive model of the cut, without any need for expensive remeshing operations (Figure 4.22). An essential advantage of the method is that the cut surface is generated from mechanical information and not just geometric data. The model was employed during an automatic needle insertion inside the liver during the respiratory motion (see Fig. 4.23). A robotic motion has been precomputed to compensate for breathing motion. The needle insertion is performed with both a flexible (E f lexible = 50GP a) and a rigid (E rigid = 200GP a) 12cm long needle following the same precomputed motion of the robot. As expected, the flexible needle deforms and does not reach the threshold to cut the tissue. Instead, the rigid needle introduces significantly higher stress in the volume, resulting in the cut of the tissue. The cutting path is therefore extended, allowing the needle to slide along the cut direction.

Perspective

The future directions of this work will be: i) test and apply the method to an experimental environment ii) better numerical models and parameterization iii) test new control strategies based on extreme learning machines.

Shared Control and teleoperated controller The current developments are focused on the implementation of the experimental setup to test and validate the shared control approach (Figure 4.24). We already managed to test the automatic control strategy in an actual experimental setup, but it was limited to quasi-static scenarios. We are currently developing a new setup composed of a KUKA arms LBR IIWA14. A specific needle gripper (pneumatic actuator) has been developed to grasp and release the needle in case of complications immediately. The registration of the robot, needle, and silicone gel is still performed with the optitrack system. We use a haptic interface to control the position and orientation of the robot's endeffector (Omega 6 from force dimensions) and share the control between the automatic commands and user inputs.

In terms of research, the user study showed encouraging results showing that the method allows for a tip placement precision equivalent to an automated needle insertion for untrained users. Trained users can aim with more precision thanks to the visual feedback, reducing the in-plane precision. However, further research is needed to validate the approach, including a more populated user study and an actual robotic system. In addition, the KUKA arms LBR IIWA14 are redundant robots. In the future, we will investigate the possibility of taking the kinematic redundancy in the inverse FE loop to improve the insertion path and stability.

Integrated Needle/tissue interaction model and parameterization:

The elastic formulation described before offers the possibility to model needles with various diameters, stiffness, and Poisson coefficients, but it suffers from the fact that the length of needles is challenging to preserve. It represents an important issue since the variation of the length of the needle may significantly influence the Jacobian of the simulation, especially for rotations that are currently very limited. Soon, we will investigate the possibility of extending the Cosserat model and adapt the formulation of our needle/insertion interaction simulation. For this purpose, we will propose specific constraints to couple the Cosserat model with the volume FE mesh of the liver. Specific considerations will be necessary to model the nonlinear behavior of the needles, such as the friction along the shaft and the penetration force when the needle punctures the tissues. The model will also be used in the control loop and compared with the current strategy.

Among the numerous parameters used in FE simulations, the behavior of the needle and the haptic feedback is strongly influenced by the number and the distance between constraints d max . If too large discretization is chosen, it may result in "free" areas where the relative motion between the needle and the tissue is not constrained, which is unrealistic. On the other hand, if too many constraints are introduced, it may significantly slow down the performances, and more importantly, it may create an over-constrained problem leading to unsolvable situations. We will propose new solutions to integrate constraints over all the volume intersected by the needle. Similar to Allard et al. (2010) proposed for collision detection and response at the arbitrary resolution, we will propose to group constraints along the shaft of the needle and report their influence on the models numerically. This solution should entirely prescribe the displacements of both models while allowing to control the number of constraints and find the best trade-off between accuracy and computation time. We expect difficulties to provide a consistent definition of the orientation of the constraint groups. Indeed, the misalignment of constraints is known to raise convergence issues during the constraint solving process. Therefore, we will propose numerical strategies to recompute the constraint groups' directions during the solving process to enforce stability. Beyond modeling, we are also interested in the possibility of parameterizing the interaction model with accurate data. However, some parameters are specific such as the friction or the cutting threshold. For this purpose, we are currently working on the possibility to detect "events" in the force measurement to control the simulation state (Figure 4.25). For instance, transitions and ruptures can be detected in the insertion force profile as proposed in Barbé et al. (2007b). In this case, the penetration's detection can replace the penetration force λ max , which is unnecessary. We wish to extend the method in case of friction and to detect cutting events. However, we insist that our goal is not to estimate the ground truth parameters but to fit the experimental data according to the chosen models as much as possible.

Learning robotic needle steering inverse finite element simulations: Tissue motion compensation during automatic needle steering is a challenging research topic. While simulation-based control strategies capture the deformable nonlinear coupling between needle and tissue, they significantly increase the computational cost of the control. In the future, we will investigate the possibility of relying on machine learning methods to enable autonomous robotic needle steering with fast computation times. We propose to use an Extreme Learning Machine (ELM) to learn an inverse model which accounts for needle-tissue interaction. The ELM is trained with synthetic data generated from multiple needle insertions controlled by the iFE simulations method introduced in this document.

The use of ELM to perform needle steering into soft tissue simulations presents promising perspectives. Preliminary results indicate the method can achieve compatible clinical precision, and it is robust to previously unseen trajectory shapes and variable tissue elasticity parameters while using only a third of the computational time demanded simulation-based methods.

In future works, the ELM will be tested in scenarios with external perturbations to address its robustness. Finally, an experimental validation in an anatomical phantom will be performed to address the method's viability in physical scenarios.

Conclusion

This chapter proposed a new approach for automatic control of robotic needle insertion in a deformable environment. Input commands are derived from the inverse FE simulations allowing for the prediction of the behavior of deformable structures. Errors of FE models to real material are controlled thanks to a non-rigid registration performed at high frequency. We proposed a constraint-based formulation allowing for the computation of inverse steps in constraint space providing this way total insertion time compatible with clinical applications.

A shared-control method was also proposed, allowing for percutaneous needle placement inside of soft tissue. The method allows the user to have full control over the tip position inside the tissue through a haptic device, while complex tasks are performed automatically: i.e. the difficult needle steering required to follow the tip reference is done automatically, whereas the user is guided with virtual fixtures to keep the needle on a predefined path. The user study showed encouraging results showing that the method allows for a tip placement precision equivalent to an automated needle insertion for untrained users. Trained users are capable of aiming with more precision thanks to the visual feedback, reducing the in-plane errors.

Research projects and Related activities This chapter presents the funded research projects and collaborations that have supported the contributions described in the previous chapters. Then we discuss technological developments and transfer activities. All developments are carried out on the open-source SOFA platform. Finally, we detail some related activities and other responsibilities.

Main research projects

Between 2013 and 2021, I participated in 10 funded projects, either as a principal investigator or as a collaborator, representing an amount of nearly 1.2 Me in total. These projects were obtained from very selective organizations. The most critical grant is the ANR JCJC Sperry project, which perfectly fits the objectives described in this manuscript and allowed me to hire and be responsible for a group of students composed of a postdoctoral fellow, two doctoral students, and a research engineer.

CONECT: Coupling rObotic and simulatioN for nEedle proCedures auTomatized, IDEX grant, (180 Ke): This work uses finite element simulations to control a medical robot during needle insertion into a deformable gel. The robot, the needle, and the gel can be considered as a deformable robot with an infinite number of degrees of freedom. This project uses finite element simulations to anticipate deformations and propose new control strategies for the robotic system. LOSAR: Liver Open Surgery with Augmented Reality, Inria grant, (140 Ke): Our goal is repeatedly applying our method for one or several essential publications in medical conferences. However, current steps are still insufficiently automated, and the algorithm needs to be improved for better reliability. These essential elements are outside the mission's research and require development and engineering effort. Indeed, an effort for automation and ergonomics will have to be made to make the software sufficiently simple to be used in the operating room. The accuracy of the registration model should also be verified and validated through experimentations.

Sperry: SuPervisEd Robotic suRgerY -Application to needle insertion, ANR JCJC grant, (300 Ke): We aim at developing an autonomous robotic system for RFA of deep liver tumors. This project follows the Conect project, which represents the first steps towards our objective. However, many simplifications had to be done, which prevented using the method for a percutaneous insertion in living tissues. Mainly, the current solution was limited to a quasi-static scenario, and we did not consider the possibility that deformations may be generated by external factors (contact with other organs, breathing. . . ). It represents essential scientific obstacles that are addressed in this project. The Sperry project aims to improve the numerical models, computation time, and total insertion time essential for experimental validation. Finally, we have a solid motivation to bring our research results to an experimental prototype allowing for qualitative and quantitative evaluation.

Technology development

The code developed in my research activity is integrated into SOFA either as private plugins or directly in the open-source framework. SOFA has already been the leading development framework for more than 20 successful Ph.D. students. In addition, through our previous projects, we have learned how to benefit from medical inputs and expectations. To have a maximum impact, SOFA allows the design of medical applications quickly and includes medical expertise and the fundamental research and development phases. My activity related to robotics also led me to use and contribute to the robotic IRIS platform at ICube.

SOFA framework

Software Open Framework Architecture (SOFA) Faure et al. ( 2012) is a globally leading open-source platform for interactive medical simulation. This toolkit is highly modular, with separate models for constitutive mechanics, visual rendering, and haptic collision detection. Its architecture emphasizes separate models for constitutive mechanics, visual rendering, and collision detection (Figure 5.1). SOFA exploits mapping software classes, depicted with arrows, to coordinate these three models. Moreover, it provides a scene graph designed explicitly for interactive medical simulation and offers several competing components enabling users to interchange them flexibly to optimize the result. SOFA also supports additional functionality by a highly flexible plugin mechanism

The Inria Foundation administratively hosts the SOFA Consortium as a nonprofit organization. The first objective of the SOFA Consortium is to maintain the code of SOFA and manage the coordination of all community developments. The consortium also fosters the growth of the SOFA ecosystem by providing training and support, thus increasing the size of the open-source community. Finally, the SOFA consortium promotes all SOFA-based technologies and eases the associated technology transfers. Therefore, in the context of a research project, the SOFA Consortium guarantees code stability, industrialization, and the best accelerator of technology transfer.

IRIS platform

The IRIS platform (Imaging, Robotics, and Health Innovation) of ICube laboratory is an experimental platform dedicated to developing and implementing medical and surgical aid solutions. The platform offers services for designing and implementing robotic systems and software integration to control robots and imagers. IRIS is organized around five areas: Human Imaging, Preclinical Imaging, Robotics, Biomechanics, and Information / Image Processing, grouping experimental means and human resources for the accompaniment of research in GMCAO (Gestes Médico-Chirurgicaux Assistés par Ordinateur). The robotics department offers state-of-the-art robotics equipment, including a collaborative robotic cell (two KUKA arms LBR IIWA14 and accessories acquired thanks to the IDEX University of Strasbourg) and a robotic micrometric positioning cell that allows for a robotic micrometric positioning cell the development of systems for needle insertion.

The platform is located in the compound of the civil hospital of Strasbourg, near IHU, IRCAD, surgeons, and radiologists of the NHC, which facilitates interactions with the actors of the medical domain. The IRIS platform is part of the national equipment and robotics network Equipex ROBOTEX (thématiques Robotique Médicale et Robotique dans l'Usine du Futur) and is also exploited within the framework of the Labex CAMI.

The Sperry project is developed within the IRIS platform. I co-financed, and I supervise a research engineer who works on the project demonstrator, keeping in mind the possibility of sharing and reusing (via the IRIS platform) the results for other research projects. It includes robotic developments such as the needle gripper, camera calibration, registration, and communication protocols between the KUKA robot and SOFA.

Technology Transfer and Valorization

Outcomes of the research described in this manuscript have high transfer potential. Indeed, this project's technical and industrial results are particularly relevant for digital information and computer-assisted surgery companies, and medical imaging companies. Therefore, it is essential to facilitate exchanges and interoperability of software without sacrificing the possibility to commercialize our results.

I authored 4 APP repositories. Two of these modules (related to numerical solvers) were transferred to InSimo in 2013. In August 2018, another plugin for needle insertion simulation was licensed to Marion Surgical (Canada) through the SATT Conectus.

In addition, we are working with Infiny Tech 3D to integrate my research plugins in demonstrators and increase their visibility to get new contracts with other industries and encourage new collaborations that can lead to technological transfers. Supervision of interns: Finally, since I got my position in 2013, I supervised eight interns. On average, each intern worked in the team for six months. All the students passed their defense with excellent ranking, and they all found a job very quickly related to the topic of numerical simulations. My expertise and review activity concerns several domains, including numerical simulations, non-rigid deformation, and robotics. On average, I'm reviewing between 10 and 20 papers per year.

Other research activities

Teaching activity

I took the responsibility of creating a teaching module focusing on the Real-time simulation of biomechanical models. This module introduces the main challenges of biomechanical simulations, augmented reality, and robotic control, emphasizing numerical aspects and optimizations, allowing reaching computation time and performances compatible with interactivity and real-time. This module (around 30 hours) is given to 4 groups of students (TI Santé, Master Imagerie, Robotique et Ingénierie pour le Vivant (IRIV), Formation d'ingénieurs en partenariat (FIP), and MECC Master informatique. The complete material of this class is available online at https://hadrien.courtecuisse.cnrs.fr/home/teaching/.

Each group of students has its specificity. The TI Santé students are biomedical engineers, where we try to implement in C++ a complete simulation loop including numerical solvers and deformation models. The master IRIV gathers biomedical engineers and medical doctors. My lecture is, therefore, more oriented toward medical applications and modeling. The FIP specialty includes students that were half employee and half students. Finally, my lecture focuses on computer science and programming issues with the master's students Math Info.

In the coming years, I will be involved in ITI HealthTech, characterized by a great diversity between doctors and surgeons, researchers in computer science and engineering sciences, and creativity and innovation management.

Conclusion

The main objective of this research project is to advance towards the use of simulation as a tool for training and assistance in medical procedures. We identified three levels of applications of our work, lying at the intersection of several scientific domains.

Real-time simulation of complex FE models and interactions: Chapter 2 focused on the use of numerical methods to speed up the computation times of simulations. We were particularly interested in the contacts and interactions between instruments and deformable parts since these steps are the most expensive ones in the simulation loop. We proposed a method of asynchronous preconditioning which makes it possible to deal with contact problems (including nonlinear models) in real-time.

The ability to simulate surgical cuts, dissection, soft tissue tears, or microfractures is essential for most surgical procedures. For this purpose, it is necessary, on the one hand, to update the topological information of the mesh and, on the other hand, to update mechanical information in real-time.

We proposed a new real-time remeshing algorithm. The proposed method is generic and makes it possible to limit the number of nodes and enforce the quality of the elements created. We also introduced a numerical method, based on the asynchronous preconditioner, to consider the cut from the numerical point of view.

Augmented reality and biomechanical registration: In Chapter 3, we introduced a generic constraint-based framework for the non-rigid registration of a preoperative model with intraoperative data. To regularize the problem, we extract a sparse (and potentially noisy) set of data from live images later to impose displacement to the intra-operative biomechanical models. The information collected from images is then used to minimize errors' model during the registration step and find a standard coordinate system for the patient, the imaging system, FE models, and potentially the robot.

Regarding intra-operative assistance, we proposed a complete registration for brain surgery. In collaboration with the Sintef in Norway, we developed a neuro-navigation system for updating the preoperative MRI images and the localization of a tumor during the resection of the tumor. A linear elastic model combined with image-based constraints was used to register vessels extracted from the preoperative MRA and intra-operative Doppler ultrasound images. Results showed that the biomechanical registration was able to predict displacements of internal structures (landmarks, tumors) with the accuracy of less than 2 mm errors even if few or no intraoperative data were available. The originality of this work lies in the fact that we only use data compatible with current clinical routines, allowing this way to consider short-term medical applications.

Per-operative human liver registration is another important topic. We developed a method to display in AR the internal structures of the preoperative scan on top of the intraoperative images. We proposed a markers-based method allowing for the tracking and the deformation of a preoperative model in real-time during open liver surgery. Our framework is composed of both a non-rigid initial registration method to estimate the location of the markers in the preoperative model and a real-time tracking algorithm to deform the model during the surgery. The method was validated on both synthetic and ex-vivo samples. In addition, we demonstrated its applicability in the operating room during a liver resection surgery on a human patient.

Robotic Control Using Inverse Finite Element Models Chapter 4 concerns the development of numerical methods for the control of a robotic system interacting with deformable structures. The most original aspect is to rely on inverse Finite Element simulations to control an articulated robot interacting with deformable structures. We showed that it was possible to derive the Jacobian of the simulation in real-time, relating the Cartesian displacements of the base of the needle with the displacements of the tip inside the volume. Then, the Jacobian was used to compensate, or even induce, necessary deformations to reach a target. The ability to predict the behavior of deformable structures is particularly relevant for needle insertion applications since the path taken by the needle is highly impacted by the whole history of the insertion and must be anticipated as soon as possible.

Our method is based on the computation of several forward simulation steps where each DOF of the articulated robot is successively perturbed to measure the corresponding displacement of the tip inside the tissue in the simulation. Specific numerical strategies were proposed since several steps of a Finite Element simulation involving complex nonlinear constraints are necessary. Nevertheless, since the Jacobian is used in the control loop of the robot, it is necessary to optimize this process in order to limit delays that may create instabilities to control the system. We showed that our method could provide new updates of the Jacobian at frequencies between 40 -100 Hz, which is sufficient to enforce the system's stability.

We applied this method to an actual experimental setup. The needle was inserted inside a polyurethane foam using a Mitsubishi RV1A anthropomorphic robot arm. During the insertion, vertical and lateral deformations were generated, leading to a significant modification of the undeformed trajectory, necessary needle bending, and even an off-plane shift between the base of the needle and the insertion point. Despite these substantial modifications, the method maintained the tip of the needle within the thickness of 1 cm of the foam and followed the desired curved path with an accuracy lower than 1 mm without any human intervention.

Perspectives Combining the research topics addressed in this manuscript allows us to envision large-scale applications with actual medical outcomes. The perspectives of this work are therefore numerous.

Regarding numerical aspects, we are currently investigating the possibility to solve constraints in an intermediate space, called isolated dofs, much more suitable for GPU parallelization and to reuse the computed information between time steps. Other strategies such as domain decomposition or model order reduction will be investigated to significantly increase the dimension of FE meshes that can be simulated at interactive rates. Finally, we wish to develop innovative solutions for the simulation of cuts and collisions after cutting. Although it is essential in most surgical practice, it still represents an open scientific challenge.

Augmented reality and non-rigid registration methods for visual assistance are of significant importance for the surgical community. In the coming years, we plan to apply our solution to actual fluoroscopic medical data.

Numerous engineering issues will be addressed owing to collaborations with IHU Strasbourg. An important question that still needs to be answered is the availability of providing metrics and quality assessments in order to make medical decisions based on the data provided by the models.

The control of medical robots for percutaneous procedures is also an important topic that will be continued. In the short future, we aim to share control strategies between the user and the autonomous system to allow animal experiments. For this purpose, many scientific and engineering steps still need to be overcome. We will also study the possibility of steering a needle using a learning-based approach. Machine learning control strategies can shift the workload to an offline training phase. It is then possible to apply the results of the patient-specific inverse simulations and derive control commands without any simulation in the operating room. It presents critical scientific obstacles that we wish to address.

  . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Interactive medical simulations . . . . . . . . . . 1.1.2 Computer Assisted Interventions . . . . . . . . . 1.1.3 Percutaneous procedures . . . . . . . . . . . . . . 1.2 Contributions and medical applications . . . . . . . . . 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.1: RFA at Strasbourg Hospital. Radiologists use CT volume reconstruction to position the needle at the target location.

Figure 1 . 2 :

 12 Figure 1.2: Percutaneous procedures under fluoroscopy. The surgeons have to place the needle's tip at a pre-identified location using 2D projective fluoroscopic images, providing no information about the depth of structures.

Figure 1 . 3 :

 13 Figure 1.3: Clinical workflow:The pre-operative steps are ideally achieved several days before the surgery. Based on a tomographic image, we acquire the geometry of the organ and plan the trajectory for the needle insertion. Intraoperatively, the robotic system automatically performs the needle insertion. FE models, including needle/tissue interactions and collisions with organs, are used to predict the behavior of internal structures. To limit errors, FE models are permanently registered with intraoperative data extracted from medical images. Input commands of the robot are computed automatically in real-time with an inverse simulation using the registered FE models. In addition, for safety reasons, the user shares the control between the automatic method and remote controllers.

Figure 2 .

 2 Figure 2.1: Application for biomechanical models. Predictive and accurate models are necessary for a clinical usage, but strong computation time constraints must be enforced for any type of application.

Figure 2 . 2 :

 22 Figure 2.2: Corotational FEM. Left: a tetrahedral element in the reference configuration with its local frame. Middle: the deformed tetrahedron and its local frame. Right: deformation (blue arrows) as measured after aligning the frames.

  model  FE model

Figure 2 . 3 :

 23 Figure 2.3: Contact mapping. The collision detection is performed on the collision meshes (orange). The proximity information are defined at arbitrary location on the collision surface. They are transferred to equivalent nodal forces, first on the collision surface and then on the FE mesh using mapping functions (red).

Figure 2 . 4 :

 24 Figure 2.4: Constraint linearization.

Figure 2 . 5 :

 25 Figure 2.5: Illustration of the physical phenomena simulated for the contact's response.

Figure 2 . 7 :

 27 Figure 2.7: Contact force distribution in different scenarios and using different approximations of the mechanical coupling. 2.7a is homogeneous, 2.7b is heterogeneous and contacts are solved without coupling, 2.7c is heterogeneous and contacts are solved with coupling.

Figure 2 . 8 :

 28 Figure 2.8: CPU/GPU task graphs.

Figure 2 . 9 :

 29 Figure 2.9: Euler implicit + CG

Figure 2 . 10 :

 210 Figure 2.10: Performance comparison between CPU and GPU in interactive FEM simulations with varying object complexity (number of tetrahedra).

  Figure 2.11: The preconditioner is updated asynchronously within a dedicated CPU thread. We use the last preconditioner available to advance the simulation so that the simulation never needs to wait for the computation of the current preconditioner to be complete.

Figure 2 . 12 :

 212 Figure 2.12: An estimated rotation matrix is computed for each node, between updates of the preconditioner.

Figure 2 . 13 :

 213 Figure 2.13: Performances of different preconditioning schemes during the simulation of a co-rotational FEM beam bending under gravity, measured as (a) the number of iterations required by the CG to converge, and (b) the computation time per simulation time step.

Figure 2 . 14 :

 214 Figure 2.14: Simulation of a heterogeneous deformable beam falling under gravity and average computational time for 200 simulation steps and convergence rate for different preconditioners. Red parts are 50× stiffer than blue parts. Inverse corresponds to the inversion of the diagonal matrix for the Jacobi preconditioner and the factorization of the system for LDL T preconditioners. Solving is the time taken to solve the system. Total is the total time of a single time step.

Figure 2 . 15 :

 215 Figure 2.15: First level of parallelism achieved for solving a Sparse Triangular System with multiple right-hand side vector on the GPU.

Figure 2 . 16 :

 216 Figure 2.16: Computation time for solving a STS with multiple right-hand side vectors.

Figure 2 .

 2 Figure 2.17: A set of tetrahedra before (left) and after(right) the cut. The re-meshing is based on the 2-3 flip operation, which flips edges of intersecting tetrahedra.

Figure 2 . 18 :

 218 Figure 2.18: Incremental update of the mesh structure for the cut.

Figure 2 . 19 :

 219 Figure 2.19: Correction of the preconditioner during topological modifications. When a modification is performed on the mesh, we first compute the correction of the current factorization. Then we compute the correction of the preconditioner, which was being calculated at the time of the cut. After two consecutive updates without topological modification, the preconditioner does not need any additional correction.

Figure 2 . 20 :

 220 Figure 2.20: Multithreaded approach for the haptic rendering of hyperelastic models. A preconditioner representing an approximation of the compliance matrix is updated at low frequency. A local contact problem is updated at a haptic rate. The contact directions and the compliance are assumed constant, and only the violation δ is updated according to the device's motion.

Figure 2 . 21 :

 221 Figure 2.21: Simulation of the lens extraction using MSICS technique.

Figure 2 . 22 :

 222 Figure 2.22: Simulation of a hepatectomy with haptic feedback, and associated datasets used. (left) FE mesh used, (middle) boundary conditions from the collision detection, (right) simulation of cutting.

2. 4 . 3 . 2

 432 Application to brain tumor resection Surgical resection of a brain tumor is an active research topic, for instance for the prediction of the brain shift Joldes et al. (2010) or for the Deep Brain Stimulation Bilger et al. (

Figure 2 .

 2 Figure 2.23: Real-time simulation of a brain tumor resection.

Figure 2 . 24 :

 224 Figure2.24: Evolution of contact space as well as isodof in a collision simulation between a deformable liver mesh and a rigid plane. For each 6 time steps, the blue elements represent the isolated DOFs appeared in the previous steps, and the red elements are the new isolated DOFs; The isodof elements are reflected in the matrix above the liver.

Figure 2 .

 2 Figure 2.25: DDM: Each color can be computed in parallel. Borders (black) must be processed sequentially.

  Simulation of cutting: We are currently developing a new method to turn the problem of remeshing into a non-rigid registration problem (see next chapter). We first generate the cutting surface from a point cloud representing the cutting path of the instrument. The point cloud is meshed using the Point Cloud Library providing this way a 3D surface composed of triangles describing the cut. Elements crossing the cutting path are duplicated. Then, we rely on node snapping technique to move the duplicated points onto the cut surface. For this purpose, an elastic registration problem is solved to move the neighborhood nodes. Finally, the topology structure is changed, and the mesh separates into two parts.

  Figure 2.26: Remeshing with constraints.

Figure 3 .

 3 Figure 3.1: Augmented reality using biomechanical models. A preoperative biomechanical model is built from preoperative data and merged interpretively with data extracted from intraoperative images in order to add information on the screen.

  Biomechanical simulation-based registration methods: Clatz et al. (2005) and Vigneron et al. (2012) performed MR acquisitions during surgery.

  p are the positions of the model S, and T (S) corresponds to the transformed model observable in the same condition as the per-operative data. For instance, in the context of laparoscopic images, T (S) corresponds to the projection 1 of the deformed model on the camera frame.

Figure 3 . 2 :

 32 Figure 3.2: Multiple solutions: both configurations have the same projection in the camera view and minimize the energy between internal forces of the FE model and projective constraints. Both of them are a local minimum of equation (3.1).

Figure 3 . 3 :

 33 Figure 3.3: Initial Problem: A set of data can be extracted on the intraoperative shape, but it's difficult to relate their position with the pre-operative model.

Figure 3 . 4 :

 34 Figure 3.4: Binding process of the control points and the FE structure. 3.4a the control points (green) are associated (red lines) with the closer surface of the object (orange). 3.4b the distance criterion is not sufficient since all the control points are associated with the same side of the object and the constraints cannot be satisfied. 3.4c, a set of normals aiming inside the surface are computed. 3.4d the control points are associated with the closer triangle whose normal is oriented in the same direction.

Figure 3 . 5 :

 35 Figure 3.5: The broken blue lines represent the Bezier surface defined over the triangles of the source mesh. m (in red) are paired with positions p i (in blue) at the beginning of the step. After the resolution, p f (in green) must be located on the tangential plane defined by the constraints. Note that m does not move, but instead, the model is deformed (assuming linear translations of the solution plane) to enforce the constraint.

Figure 3 . 6 :

 36 Figure 3.6: Projective constraints definition. The visible outline of the model is extracted with respect to the camera position 3.6a. An ICP method is performed to bind the segmented contour in the image with the projected contour of the 3D model 3.6b. Finally constraints are projected back on the 3D model without any constraints along the depth of the camera 3.6c.

Figure 3 . 7 :

 37 Figure 3.7: At the end of each time step, the constraints are satisfied i.e. m i and p i are located on the tangential plane (blue) given by the normal (red) of the nearest triangle.However, the cumulative effect of other constraints may move the solution away from the closest projection. 3.7a, 3.7b, 3.7c and 3.7d show an example of consecutive time steps, respectively before and after the constraint resolution.

Figure 3 . 8 :

 38 Figure 3.8: The binding process (left) and constraint forces evaluation (right) to register a deformable object (blue) with the control surface (brown). The Gauss-Seidel algorithm iteratively activates (red) or deactivates (gray) the constraints according to the respective violations.

  shape computation for deformable model of brain In Morin et al. (

Figure 3 .

 3 Figure 3.10: (a) simulated models (gray) with rest positions (blue) and target positions (red) (b) error repartition on FE mesh from yellow (close to 0 mm) to dark blue (2.69 mm) (c) von Mises stress (in Pa) on a slice of the mesh at equilibrium.

Figure 3 .

 3 Figure 3.11: Brain-shift compensation

Figure 3 . 12 :

 312 Figure 3.12: Brain shift compensation using constraint-based simulation. Pre-and intraoperative blood vessels, respectively in blue and orange, are registered. Other constraints maintain the brain in sliding contact with the dura mater and under the US probe. The Hausdorff distance between landmarks and delineated anatomical structures manually segmented are shown in graphs 3.12b and 3.12c.

Figure 3 .

 3 Figure 3.12a shows an example of simulation. Figures 3.12b and 3.12c show the corresponding errors for 4 patients. On average, 72% of the deformation

Figure 3 . 13 :

 313 Figure 3.13: Example of preoperative MRI (left), B-mode intraoperative US (center) and MRI updated with our method (right).

Figure 3 . 14 :

 314 Figure 3.14: Execution time of the preoperative steps (without taking into account the time required for segmentation).

. 15 .Figure 3 . 15 :

 15315 Figure 3.15: Execution time of the intraoperative stages, after acquisition of the ultrasound images and 3D reconstruction of the image volume.

Figure 3 . 16 :

 316 Figure 3.16: Compensation during resection.

3. 4

 4 .2.1 3D Physics-Based Registration of 2D Dynamic MRI Data In Courtecuisse et al. (2014b); Trivisonne et al. (2016), we introduced a method for the registration of 3D models with 2D slice images.

  Figure 3.17: Outline tracking of the liver with openCV.

Figure 3 . 18 :

 318 Figure 3.18: Constraint slices.

Figure 3 .

 3 Figure 3.19: 3D registration of the liver's model with the 2D dynamic MRI slices. Large deformations and rotations are applied to the mesh, and the method remains stable and real-time during the breathing cycle.

3. 4

 4 .2.2 3D Physics-based Registration of Pelvic System Using 2D Dynamic MRI Slices In Courtecuisse et al. (

Figure 3 . 20 :

 320 Figure 3.20: Constraints set used in the simulation: Outline constraints (left), Contact Constraints (center), anatomical constraints (right).

3. 4 . 3

 43 Non rigid registration for liver surgery 3.4.3.1 Silhouette-based Pose Estimation for Deformable Organs -Application to Surgical Augmented Reality

Figure 3 . 21 :

 321 Figure 3.21: Overview of the methodIn[START_REF] Adagolodjo | [END_REF] we introduced a method for semi-automatic registration of 3D deformable models using 2D shape outlines (silhouettes) extracted from a monocular camera view(Figure 3.21). The method is composed of the following steps: 1) A direct simulation D(S) is applied to transform the reconstructed model obtained from the segmentation (red) in shape close to the 3D position observed in the image (green). 2) A Rigid transformation (blue) is provided by the user to roughly align the model with the contour of the organ segmented in the image (yellow). 3) Projective constraints are applied to the biomechanical model to fit the organ's contour and finally provide the 3D shape in the camera's frame.

Figure 3 . 22 :

 322 Figure 3.22: Application to liver surgery (left) and kidney (right).

3. 4 Figure 3 . 23 :

 4323 Figure 3.23: Workflow of the method.

Figure 3 .

 3 Figure 3.24: (Top) Augmented view projected on the monitor for different steps of the surgery.

3. 4

 4 .4.1 Model-Based Identification of Anatomical Boundary Conditions in Living Tissues In Peterlik et al. (

Figure 3 . 25 :

 325 Figure 3.25: Evaluation of the method on porcine liver deformation induced by repositioning the pig from supine to flank positions. Predicted surface loads (b,c).

Figure 3 . 26 :

 326 Figure 3.26: Main steps of the atlas creation process and its application to the transfer of boundary conditions.

Figure 3 .

 3 Figure 3.27: (Left) Clinical workflow. (Right) Registration of a liver with virtual data.The imaging plane is represented along with the position of the radiating source (in blue), the internal marker's projections (black dots), and the resulting projection lines (red lines). To enforce registration, the 3D markers (green dots) are constrained on the projection lines using 2D projective constraints (blue and green arrows). The mechanical solver then finds the minimal energy position of the liver, enforcing the constraints.

Figure 3 . 28 :

 328 Figure 3.28: General workflow of the proposed solution. The method computes asynchronously a confidence map of the registered model applying forces on the registered configuration in various directions ⃗ x, ⃗ y, ⃗ z in order to estimate the mobility of the model. The mobility map is then combined with the augmented view of the organ only to display the reliable parts of the model fitting the desired accuracy for the surgical procedure.

Figure 3 . 29 :

 329 Figure 3.29: Confidence map associated with the desired accuracy of 1 cm.

Figure 4 . 1 :

 41 Figure 4.1: Needle insertion inside a deformable environment. (a) A curved trajectory is defined (dashed lines) to avoid an obstacle (gray circle). (b) Since the needle is stiffer than the foam, it is necessary to deform the foam (with tangential motion of the needle's base) in order to deform the trajectory to a straight line aligned with the needle's shaft.

  , being either teleoperated masterslave systems Piccin et al. (2009); Hiraki et al. (2020) or even passive guides Fischer et al. (2008); Arnolli et al. (2018); Han et al. (2019); He et al. (2020).

Figure 4 . 2 :

 42 Figure 4.2: Control loop: The robot is controlled by an inverse FE simulation; itself registered with intraoperative images.

Figure 4 . 3 :

 43 Figure 4.3: Assembly of the matrix of the system described by a deformable object composed of beams connected in series.

Figure 4 . 4 :

 44 Figure 4.4: (Top) the real needle and the model under different loads (from 6.5 to 24 grams). (Bottom) The back-projection errors (pixels) in images between the real and the virtual needle deformation for various parameters.

Figure 4 . 5 :

 45 Figure 4.5: Constraints applied during needle insertion simulations. Blue arrows are bilateral constraints. Red arrows are friction constraints. The green arrow is a unilateral constraint. Fixed parts of the gel (screws) are shown in black.

Figure 4 . 6 :

 46 Figure 4.6: The three stages of penetration: stage 1, δ > 0 no contact force is applied (λ = 0). Step 2, there is contact, a positive force (λ > 0) is applied to cancel the penetration (δ = 0). Until λ < p f the objects are considered in contact which leads to deformation. Step 3, if λ ≥ p f the needle penetrates.

Figure 4 . 7 :

 47 Figure 4.7: Trajectory constraints are applied along with the closest respective points on the needle. This pairing is recomputed at each time step to let the needle advance along the trajectory.

Figure 4 . 8 :

 48 Figure 4.8: (left) Experimental setup. (Right) view from the top and front cameras.

(

  ii-a, ii-b) -Monocular camera: We are using two monocular cameras (Logitech Webcam C920). The positions of the optical cameras are estimated, solving the Perspective-n-Point (PnP) problem based on the location of 3D positions provided by optitrack and the 2D coordinates manually segmented in the images. In the figure 4.8, red rectangles indicate the region of interest of each camera where the calibration was performed. After the registration step, the 2D back-projection error of optical markers is under 1 pixels in the whole zone of interest.

Figure 4 . 9 :

 49 Figure 4.9: Needle support.

Figure 4 .

 4 Figure 4.10: Denavit-Hartenberg model of the Mitsubishi RV1A robot.

Figure 4 . 11 :

 411 Figure 4.11: Constraints were applied during the needle insertion simulation. Each arrow represents one bilateral constraint.

Figure 4 .

 4 Figure 4.12: (Top) top camera view (Down) the front camera view. Three large deformations of the foam were manually created. After registration, the outline of the model is projected on the image and shown in blue. the end effector, p n is the positions of the needle, p v is the position of the volume, and m is the 3D position of the markers. Solving S provides at any time t, the positions p (t) n and p (t)

Figure 4 . 13 :

 413 Figure 4.13: Objective function.

Figure 4 . 14 :

 414 Figure 4.14: Overview and implementation of the proposed solution.

4. 4 . 1

 41 Robotic control using iFE simulation 4.4.1.1 Inverse real-time Finite Element simulation for robotic control of flexible needle insertion in deformable tissues In Adagolodjo et al. (

1: 2 : 5 :

 25 Free Motion:x free = A -1 b Store data(p n , p v , X ): p = p (t) , X = X (t) 3: Constraint Definition: H 4: Compute Compliance: W = HA -1 H T Compute error: e = I( X , p) 6: if dot(e, e) < ε then 7:

3 .

 3 Project back λ i in the motion space (line 11).

Figure 4 .

 4 Figure 4.15: (bottom) Simple trajectory (red) with 1 obstacle (green). (top) Complex trajectory (red) with 3 obstacles (green).

Figure 4 . 16 :

 416 Figure 4.16: Control loop diagram; HF:high frequency; LF:low frequency; S (t) = S(X (t) , p

Figure 4 . 17 :

 417 Figure 4.17: Augmented reality view at different steps of the insertion. The desired trajectory is shown in green. The needle model (black) overlays the real needle (gray) in the projective view. The color map on the foam shown Von Mises stress indicating the deformation of the model compared to the initial configuration.

Figure 4 .

 4 Figure 4.18: CT scan after the robotic insertion. The desired trajectory is shown in green and the path taken by the needle is shown in red.

Figure 4 . 19 :

 419 Figure 4.19: Distance in mm (in the forward simulation) between the needle tip and the desired position on the trajectory (corresponding to the norm of the objective function ∥ ēp ∥) according to the insertion time (s), averaged over 22 simulations.

Figure 4 . 20 :

 420 Figure 4.20: Overview of the method. The user controls the needle tip placement during robotic percutaneous procedures through a haptic interface. Virtual fixtures are provided to guide him on a predefined path. An automatic needle-steering algorithm provides a reference to the robot to follow the user's target, relying only on fluoroscopic images. Augmented fluoroscopic images representing the path are provided to the user as visual feedback, along with the registered 3D model of the organ.

Figure 4 . 21 :

 421 Figure 4.21: Trials setup.

4. 4

 4 .3.1 Interactive Finite Element model of needle insertion and lacerationFor percutaneous applications, the extents of the cuts are usually not sufficient to separate the tissue. Indeed, the compression of surrounding organs and the limited extension of the cut surface prevent the organs from mechanically separate into two parts in the abdominal cavity. Based on this observation, we introduced a new model for needle insertion and lacerations inPerrusi et al. (2021).

Figure 4 . 22 :

 422 Figure 4.22: Needle insertion and cutting: first the needle penetrates the tissue. Then needle/tissue constraints prevent the relative displacements of the models. Finally, in case of lateral motion, the cutting plane is generated from mechanical forces.

Figure 4 .

 4 Figure 4.23: (left) External view of the needle insertion with robotic assistance simulations. (right) Flexible (top) and rigid (bottom) needle.

Figure 4 . 24 :

 424 Figure 4.24: Experimental setup under development in Strasbourg.

Figure 4 . 25 :

 425 Figure 4.25: Needle insertion simulations with cutting from real data measurements.

Figure 4 . 26 :

 426 Figure 4.26: Data inputs used to train the ELM neural network: for a given time step i, the figure shows the robot position P (i), the needle tip position n tip (i), trajectory point n traj (i) sliding through the pink trajectory. Lateral noise is added to the trajectory position data to avoid overfitting over the limited trajectory shapes of the training set.

Figure 5

 5 Figure 5.1: SOFA-based interactive surgery simulation: (a) screenshot of hepatic surgery simulator; (b) biomechanics, visual rendering, and collision models.

I

  am an active contributor to the open-source SOFA framework. I mainly contributed to the GPU-related parallelization and, more generally, to numerical solvers and matrix implementations. The code developed during my research is integrated as SOFA plugins and promoted to academics and companies through the Sofa Consortium. It also highly simplifies and improves the maintenance of the code developed by Ph.D. students, postdoctoral fellows, or interns, which tend to be outdated quickly. A binary version of the code developed during my research activity can be found at this address: https://hadrien.courtecuisse.cnrs.fr/home/software/

5. 4

 4 List of publications I started publishing in 2008. In 2012, at 26 years old, I became a research scientist at CNRS in section 07/06 (computer science applied to medical systems). In 2020, I published 50 articles that received more than 1000 citations (h-index = 15). Since 2008, I produce between 3 and 4 articles per year in the most prestigious reviews and journals of the medical simulation and robotics community (Media, Pbmb) and international conferences (Siggraph, Miccai, Iros. . . ). These reviews and conferences are very selective with a rate of acceptance often lower than 30%. During these years, I demonstrated my ability to contribute to the multidisciplinary areas required for my research project. During my Ph.D., I focused on surgical simulation, with a particular emphasis on real-time modeling of soft tissues and interactions with surgical instruments. I have proposed contributions in several communities such as high-performance computing Courtecuisse and Allard (2009), computer graphics Courtecuisse et al. (2010a), medical simulations Courtecuisse et al. (2011a). In a very short period of one year as a postdoctoral researcher at Cardiff University, I managed to propose new methods for the simulation of cutting in real time Courtecuisse et al. (2014a), in the mechanical engineering community. Since I arrived in Strasbourg, I have been working on medical-oriented applications and augmented reality Morin et al. (2017d); Courtecuisse et al. (2020) with potential applications in the operating room. In addition, according to my research project, I also focus on the development of new systems of robotic assistance to the surgical gesture, by proposing new solutions of inverse simulations for robotic needle insertion Adagolodjo et al. (2019); Baksic et al. (2021).Book Chapters:[1] F. Faure, C. Duriez, H. Delingette, J. Allard, B. Gilles, S. Marchesseau, H. Talbot, H. Courtecuisse, G. Bousquet, I. Peterlik and S. Cotin. SOFA: A Multi-Model Framework for Interactive Physical Simulation. Studies in Mechanobiology, Tissue Engineering and Biomaterials, vol. 11, pages 283-321, 2012.

Table of contents 3

 of3 Perspective . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . 85 • I. Peterlik, H. Courtecuisse, C. Duriez and S. Cotin. Model-based identification of anatomical boundary conditions in living tissues. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), vol. 8498 LNCS, pages 196-205, 2014. • H. Courtecuisse, I. Peterlík, R. Trivisonne, C. Duriez and S. Cotin. Constraint-Based Simulation for Non-Rigid Real-Time Registration. Medicine Meets Virtual Reality, vol. 196, pages 76-82, 2014b. • R. Plantefève, I. Peterlik, H. Courtecuisse, R. Trivisonne, J. P. Radoux and S. Cotin. Atlas-based transfer of boundary conditions for biomechanical simulation. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), vol. 8674 LNCS, no. PART 2, pages 33-40, 2014. • F. Morin, H. Courtecuisse, M. Chabanas and Y. Payan. Rest shape computation for highly deformable model of brain. Computer Methods in Biomechanics and Biomedical Engineering, vol. 18, no. sup1, pages 2006-2007, 2015. • R. Trivisonne, I. Peterlik, S. Cotin, H. Courtecuisse, I. Peterlík, H. Courte-• M. Chabanas, F. Morin, H. Courtecuisse, F. Le Lann, O. Palombi, I. Reinertsen and Y. Payan. Resection-induced brain-shift compensation using vesselbased methods. In Proc. SPIE Medical Imaging 2018: Image-Guided Procedures, Robotic Interventions, and Modeling, vol. 10576, edited by R. J. Webster and B. Fei, page 17, SPIE, 2018. • Y. Adagolodjo, N. Golse, E. Vibert, M. De Mathelin, S. S. Cotin and H. Courtecuisse. Marker-Based Registration for Large Deformations -Application to Open

	71
	3.4.2 Physics-based registration with 2D MRI slices . . 75
	3.4.3 Non rigid registration for liver surgery . . . . . . 78
	3.4.4 Identification of Anatomical Boundary Conditions 80
	3.5 Liver Surgery. Proceedings -IEEE International Conference on Robotics and
	Automation, vol. 6, pages 4007-4012, 2018.
	• H. Courtecuisse, Z. Jiang, O. Mayeur, J. F. Witz, P. Lecomte-Grosbras, M. Cos-
	son, M. Brieu and S. Cotin. Three-dimensional physics-based registration of
	pelvic system using 2D dynamic magnetic resonance imaging slices. Strain,
	vol. 56, no. 3, page e12,339, 2020.

.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . 54 3.2 Related Works . . . . . . . . . . . . . . . . . . . . . . . 55 3.2.1 Liver motion and non-rigid registration . . . . . . 55 3.2.2 Brain-shift compensation methods . . . . . . . . . 58 3.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . 60 3.3.1 Segmentation and meshing and initialization . . . 62 3.3.2 Geometrical binding . . . . . . . . . . . . . . . . 63 3.3.3 Constraints definition . . . . . . . . . . . . . . . . 65 3.3.4 Constraint Solving . . . . . . . . . . . . . . . . . 68 3.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . 70 3.4.1 Brain-shift compensation with the intraoperative US . . . . . . . . . . . . . . . . . . . . . . . . . . cuisse, S. Cotin, H. Courtecuisse, I. Peterlik, S. Cotin and H. Courtecuisse. 3D Physics-Based Registration of 2D Dynamic MRI Data. In Medicine Meets Virtual Reality, vol. 220, page 432, IOS Press, 2016. • F. Morin, I. Reinertsen, H. Courtecuisse, O. Palombi, B. Munkvold, H. K. Bø, Y. Payan and M. Chabanas. Vessel-based brain-shift compensation using elastic registration driven by a patient-specific finite element model. In International Conference on Information Processing in Computer-Assisted Interventions (IP-CAI), 2016. • F. Morin, H. Courtecuisse, I. Reinertsen, F. Le Lann, O. Palombi, Y. Payan and M. Chabanas. Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation. Medical Image Analysis, vol. 40, pages 133-153, 2017d. • Y. Adagolodjo, R. Trivisonne, N. Haouchine, S. S. Cotin and H. Courtecuisse. Silhouette-based pose estimation for deformable organs application to surgical augmented reality. IEEE International Conference on Intelligent Robots and Systems, vol. 2017-Septe, pages 539-544, 2017.

  • Y. Adagolodjo, L. Goffin, M. De Mathelin, H. Courtecuisse and Anonymous. Inverse real-time Finite Element simulation for robotic control of flexible needle insertion in deformable tissues. In IROS, vol. 2016-Novem, pages 2717-2722, IEEE, 2016. • Y. Adagolodjo, N. Golse, E. Vibert, M. De Mathelin, S. S. Cotin and H. Courtecuisse. Marker-Based Registration for Large Deformations -Application to Open Liver Surgery. Proceedings -IEEE International Conference on Robotics and Automation, vol. 6, pages 4007-4012, 2018. • P. Baksic, H. Courtecuisse, C. Duriez and B. Bayle. Robotic needle insertion in moving soft tissues using constraint-based inverse Finite Element simulation. In IEEE International Conference on Robotics and Automation, page 1, IEEE, 2020. • P. Baksic, H. Courtecuisse and B. Bayle. Shared control strategy for needle insertion into deformable tissue using inverse Finite Element simulation. In IEEE International Conference on Robotics and Automation, pages 12,442-12,448, 2021. • P. Perrusi, P. Baksic and H. Courtecuisse. Interactive Finite Element model of needle insertion and laceration. EUROGRAPHICS 2021 (Short Paper), vol. l, pages 1-4, 2021.

  .[START_REF] Courtecuisse | Preconditioner-based contact response and application to cataract surgery[END_REF]): purple denotes a method where both the trajectory and the needle are assumed rigid in the inverse simulation. black stands for the complete resolution of the FE resolution for each inverse step. green is the method introduced in[START_REF] Adagolodjo | [END_REF]. yellow is the current method without stabilization. blue is the method introduced in Adagolodjo et al. (2019) with stabilization. red is the current method with stabilization.

	We found that the method introduced in Adagolodjo et al. (2019) (green)
	already provides a speedup of 3.5× compared to solving the complete FE
	simulation for each inverse step (i.e., computing the column of the Jacobian

Table of contents 5

 of5

.1 Main research projects . . . . . . . . . . . . . . . . . . 5.2 Technology development . . . . . . . . . . . . . . . . . 5.2.1 SOFA framework . . . . . . . . . . . . . . . . . . 5.2.2 IRIS platform . . . . . . . . . . . . . . . . . . . . 5.2.3 Technology Transfer and Valorization . . . . . . . 5.3 Other research activities . . . . . . . . . . . . . . . . . . 5.3.1 Supervision activities . . . . . . . . . . . . . . . . 5.3.2 Editorial Activities . . . . . . . . . . . . . . . . . 5.3.3 Teaching activity . . . . . . . . . . . . . . . . . . 5.4 List of publications . . . . . . . . . . . . . . . . . . . .

  This project aims to propose a biomechanical model of soft brain tissues to simulate their deformations and adapt the preoperative images of the patient's intraoperative morphology in real-time. The brain model is completed by integrating the cerebrovascular tree, reconstructed from a preoperative magnetic resonance angiography. At different times of the procedure, the surgeon can then perform a Doppler ultrasound scan of the region of interest to measure the deformations of the vascular tree. These images are used as boundary conditions for our biomechanical model, which then estimates the position and deformation of internal anatomical structures. In addition, to remain compatible with potential use in the operating room, the biomechanical registration must satisfy strong computation time constraints. Inria grant, (140 Ke): This project aims to significantly reduce the computation time of real-time numerical simulations such as those used for medicine learning or intraoperative assistance. These simulations are often based on hypotheses, simplifications, and approximations to optimize computations. These assumptions are generally acceptable in the context of generic simulations but are no longer valid when one seeks to precisely reproduce the behavior of an organ in a patient-specific context for preoperative assistance. To speed up computations, minimizing code rewrite, we propose to rely on both techniques and tools from the research work of the CAMUS team. In particular, we use the results obtained in terms of fine-grained automatic parallelization and also in terms of application development.

	currently exist in medicine.
	ASNAP: Accélération des Simulations Numériques pour
	l'Assistance Peropératoire,
	RESET: Retinal Surgery System for Training, ANR grant, (45 Ke):
	We seek to develop a complete simulation of the eye during all stages of retinal
	surgery. For this purpose, the physics-based simulation must allow making
	errors and generate surgical complications. This level of realism is equivalent
	to what is called "Level D" for flight simulators. These simulations do not

GUIDANCE: Non-linear simulations for intraoperative guidance in neurosurgery, CAMI Labex/French region grant, (110 Ke):

HiPerNav: High Performance Soft-tissue Navigation, European ITN grant, (150 Ke): The overall goal of HiPerNav is to train and educate future surgeons in the multidisciplinary field of image-guided interventions. The scientific and clinical aspects aim to develop a navigation platform to assist in treating cancer and liver metastases to improve eligibility and prognosis for liver surgery. This platform must allow complete management of surgical operations in (1) preoperative surgical planning, (2) navigation during intraoperative resection, (3) postoperative quality control.

  5.3.1 Supervision activitiesSince 2008, I have supervised 4 Ph.D. students, two postdoctoral students, six research engineers, and eight internships.Supervision of doctoral students: I participated in the supervision of two Ph.D. students. The first one was carried out by Fanny Morin on the intra-operative registration of the brain using ultrasound Doppler. Fanny completed a master's degree at ENSIMAG. Yoann Payan directed her thesis at the University of Grenoble. Her work initiated a strong collaboration with the TIMC and Sintef in Norway. Fanny defended on the 5th of October 2017 with five articles in leading international conferences on biomedical registra-tion, one publication in Medical Image Analysis, and two book chapters.Yinoussa Adagolodjo carried out the second thesis on robotic control during interactions with deformable structures for percutaneous applications. Yinoussa completed a bachelor's degree in mathematics, then a master's degree in computer science in Lille. His thesis was based in Strasbourg and directed by Michel de Mathelin. He defended on the 6th of September 2018, with six articles published in the top international robotics and medical simulation conferences. In addition, an article has been published at IEEE Transaction of Robotics which is one of the leading journals in the robotic community. I am currently supervising the Ph.D. of Paul Baksic and Ziqiu Zeng, who started respectively in October 2018 and June 2019. The first Ph.D. is focused on the robotic control for needle insertion in living tissue (including during breathing motion). The second Ph.D. is dedicated to numerical solving and parallelization strategies to accelerate the computation of the simulations. Supervision of post-doctorate students: I supervised two postdoctorate students. Phuoc Bui was funded by a project obtained by Stéphane Bordas, working at the University of Luxembourg. During this period, Phuoc published five articles on mesh fracture and error control. Maciej Bednarczyk recently started as a postdoctoral position to develop the demonstrator of the ANR Sperry project. Supervision of engineers: In 2015 I supervised Yinoussa Adagolodjo and Raffaella Trivisonne, who worked as engineers on the CONECT project. They both decided to continue as Ph.D. students in the Mimesis team. In 2015 I co-supervised Maxime Mogé in collaboration with the Camus team in Strasbourg. Maxime worked as a research engineer on the automatic parallelization of FE simulations. In January 2019, I supervised Mohamed Omar Boukhris, a research engineer working on an Augmented Reality software compatible with clinical constraints. They both were funded by an "Action de Développement Technologique" (ADT) from inria. The first one was led by Camus, whereas I was the principal investigator of the Losar project. Finally, I hired two engineers to work on my ANR JCJC Sperry project.

	Name		Years	Project	Research topic
	Yinoussa		11/2014		Idex:	Implementation of the kinematic model of the Mitsubishi
	Adagolodjo	02/2015	CONECT	RV1A robot
	Raffaella		01/2015		Idex:	Registration of biomechanical models and design of the
	Trivisonne		10/2015	CONECT	robotic system
	Maxime Mogé	01/2017 10/2018	ADT: ASNAP	Automatic parallelization of finite element simulations
	Mohamed Omar Boukhris	01/2019 12/2020	ADT: LOSAR	Augmented reality in the operating room
	Frédéric Roy	12/2018 03/2019	ANR: Sperry	Computer graphics and visualization for aug-mented reality
	Name Pedro Henrique Starting date 09/2020 Defense ANR: date Laboratory Software design, implementation and research activity Research topic Director Perrusi 09/2021 Sperry related to the Sperry project
	Fanny Morin	10/2014	10/2017	TIMC	Biomechanical modeling of the brain the brain shift for intraoperative compensation of	Y. Payan
	Yinoussa Adagolodjo	02/2015	09/2018	ICube	Automatic needle steering using inverse finite element simulation	M. de Mathelin
	Paul Baksic 10/2018	expected 10/2021	ICube	Shared control strategy between insertion commands automatic and teleoperated needle	B. Bayle
	Ziqiu Zeng	07/2019	expected 06/2022	ICube	Fast numerical solvers for stable constraint-base simulations	S. Cotin
	Name		Years	Project	Research topic
	Phuoc Bui		09/2014 10/2016	USIAS	Mesh fracture simulation with error control
	Maciej Bednarczyk	01/2021 12/2021	ANR: Sperry	Demonstrator of the Sperry project

  am a regular reviewer for several international conferences and scientific journals such as IEEE Haptics Symposium, Transactions on Haptics, The Visual Computer, Medical & Biological Engineering & Computing, Eurographics, Iros, and Icra... In addition, I have been program share for topranked international conferences such as Eurographics, Surgetica. Vriphys and Icra.

	Name	Years	Research topic
	Virginie Marec	2021	Needle insertion simulation and parameterization with a force sensor
	Pedro Henrique Perrusi	2020	Predictive filters for breathing motion
	Rayan Bounoua	2020	Validation and estimation of registration error using biomechanical models
	Julia Coste-Marin	2019	Real-time simulation of cutting
	Rami Assi	2018	Development of the augmented reality platform and its evaluation
	Alexandre Dolle	2016	Parallel collision on CPU using Layer Depth Images
	Nicolas Gautier	2016	Development of a demonstrator for retina surgery
	Asmaa Ait Hadouch	2015	Development of the communication protocol between Sofa, Mitsubishi MRV1A robot, and optitrack system
	5.3.2 Editorial Activities

I

Note that the matrix may be multiplied by mapping matrices, but this operation is ommited as stated above.

Complementarity is noted ⊥. It states that one of the two values δ n or λ n must be null whereas the other is positive.

Note that even if we compute an exact factorization of A t-h , the preconditioner remains an approximation since its computation is based on a previous configuration of the objects. We use it with delay in the simulation.

In this case, this process is known as a perspective-n-point problem which can be solved given a set of correspondences between 3D points and their

2D projectionsLepetit et al. (2009).

www.itksnap.org

www.cgal.org

Note that, depending on the type of constraints (2D, 3D or projective) the distance and the linearization should be performed accordingly. For instance, 2D norms must be employed for 2D constraints.

https://pointclouds.org/documentation/tutorials/normal_estimation.html

Note that the selection of visible triangles could also be implemented with an OpenGL rendering, selecting only the triangle whose depth corresponds to the depth buffer. However, the above algorithm was sufficiently fast for our applications.

https://www.ihu-strasbourg.eu/

http://optitrack.com/

Acknowledgements