Guillaume Mestdagh

An optimal control formulation for organ registration in augmented surgery

Keywords: Contrôle optimal -Simulation biomécanique -Chirurgie augmentée Optimal control -Biomechanical simulation -Augmented surgery ∂Ω D, ∂Ω N Dirichlet boundary, Neumann boundary

An optimal control formulation for organ registration in augmented surgery

Résumé

La réalité augmentée est utilisée en chirurgie minimalement invasive pour permettre au personnel médical de suivre en temps réel les mouvements du foie du patient. Pour mettre à jour la déformation d'un organe virtuel, une méthode de recalage élastique aligne un modèle biomécanique pré-opératoire du foie avec une surface partielle observée pendant l'opération. Tandis qu'une grande partie des méthodes de recalage élastique consistent à introduire des forces fictives dans le modèle direct, notre approche vise à reconstruire la vraie densité de forces surfaciques qui a créé la déformation observée. Nous exprimons le problème de recalage dans le formalisme du contrôle optimal, en utilisant comme variable d'optimisation la distribution de forces qui s'applique à la surface de l'organe. En permettant de définir à l'avance un ensemble de forces admissibles, cette approche favorise les champs de déplacement ayant un sens physique. Nous commençons par étudier l'existence de solutions pour le problème continu et nous calculons des conditions d'optimalité de premier ordre. Puis nous présentons la méthode d'adjoint que nous avons implémentée afin de traiter le problème numériquement. Finalement, nous validons notre méthode au moyen de cas-test liés à l'application en chirurgie augmentée. Lors de ces essais, nous mesurons l'erreur de recalage, et nous cherchons également, dans un cas particulier, à donner un sens à la distribution de forces obtenue.

Introduction en français

La naissance de ce projet résulte des discussions entre Yannick Privat, professeur à l'Institut de recherche mathématique avancée de Strasbourg, et Stéphane Cotin, chef de l'équipe Inria Mimesis. L'équipe Mimesis est spécialisée, entre autres, dans les systèmes de chirurgie assistée par ordinateur, et, plus généralement, dans la simulation biomécanique en temps réel.

Comme la plupart des organes abdominaux, le foie est un organe mou, et par conséquent il subit souvent des déformations importantes. Lors d'une opération chirurgicale, ces déformations peuvent venir du changement de position du patient, de sa respiration ou encore de la manipulation de l'organe par le médecin. Par ailleurs, au cours d'une opération faiblement invasive, le médecin manipule les organes indirectement. Pour cela, il se sert d'instruments insérés à travers la peau du patient (voir figure 1), tandis qu'un retour vidéo affiché sur un écran lui permet de voir ce qu'il fait. Le mouvement du milieu, permanent, et le champ de vision, limité, sont autant de facteurs qui rendent les opérations faiblement invasives difficiles.

Pour aider les soignants à visualiser la déformation du foie au cours du temps, des systèmes de réalité augmentée ont été mis au point. Ceux-ci produisent une vue tridimensionnelle de l'organe dans sa configuration actuelle, qui fait apparaître les structures internes, tumeurs et autres vaisseaux sanguins. Dans ce manuscrit, nous nous intéressons principalement à la méthode de recalage élastique qui met à jour la déformation du foie en combinant des données acquises avant et pendant l'opération. D'une part, un modèle biomécanique du foie est extrait d'images tomographiques réalisées quelques jours avant l'opération. D'autre part, une caméra laparoscopique fournit en temps réel un nuage de points montrant une partie de la surface du foie dans sa configuration actuelle. Le recalage vise à déterminer pour le modèle pré-opératoire une déformation qui correspond aux données observées.

Plusieurs difficultés sont liées au recalage élastique du foie. Tout d'abord, la déformation du foie doit être mise à jour en temps réel pendant le déroulement de l'opération. Notons cependant que le calcul en temps réel ne fait pas partie de nos prétentions dans ce manuscrit. Par ailleurs, les données per-opératoires, ne montrent qu'une partie de la surface du foie, et par conséquent le problème d'appariement de surfaces admet plusieurs solutions, qui n'ont pour la plupart rien à voir avec la physique du problème. Afin de calculer un champ de déplacement aussi proche que possible de la réalité, la méthode de recalage doit donc réduire au maximum l'ensemble des déplacements pouvant être associés à une observation donnée.

Chapitre 2. Revue de littérature et modélisation du problème

Les méthodes de recalage élastique reposent sur un modèle direct décrivant la mécanique du foie du patient. Ce modèle est souvent discrétisé par la méthode des éléments finis. En général, beaucoup d'attention est accordée au choix du modèle élastique décrivant les déformations du foie et ses interactions avec son environnement, si bien qu'une grande variété de modèles est utilisée dans la littérature. Habituellement, ces modèles sont de type hyperélastique. Cependant, lorsque le calcul en temps réel est requis, on optera plutôt pour un modèle dit co-rotationnel, qui permet d'engendrer des déformations nonlinéaires en ne résolvant que des systèmes linéaires. Les conditions aux frontières utilisées sont également d'une grande variété et vont de la simple condition de Dirichlet au modèle complexe, décrivant les ligaments qui maintiennent le foie en place comme des solides hyperélastiques.

Concernant le processus de recalage lui-même, il consiste à appliquer des forces sur le maillage de foie pour créer une déformation. Une grande classe de méthodes est directement inspirée des méthodes de recalage d'image. On y retrouve des forces artificielles qui attirent le maillage vers le nuage de points observé, et le modèle élastique n'y joue qu'un rôle de régularisation. Dans ce type de méthode, les forces appliquées ne reflètent pas l'origine réelle de la déformation. Elles n'ont même rien à voir, puisqu'elles sont créées par un nuage de points qui n'existe pas réellement. D'autres approches se soucient davantage de savoir si le champ de déplacement engendré est plausible au vu de la physique. Dans cette deuxième classe de méthodes, un ensemble de forces ou de déplacements admissibles est défini à l'avance en se basant sur des hypothèses physiques. Les déformations calculées sont ainsi plus réalistes.

Notre approche s'identifie plutôt à la seconde catégorie. Nous exprimons le problème de recalage comme un problème de contrôle optimal. La formulation qui en découle est Ici, le champ de déplacement u est obtenu en appliquant la distribution de forces g sur la surface du maillage. L'ensemble des forces admissibles est défini par

G M = g ∈ L ∞ (∂Ω N ) ∥g∥ L ∞ (∂Ω N ) ⩽ M .
Le problème de contrôle optimal fait intervenir la fonctionnelle d'attache aux données J(u) = 1 2 ˆΓ d 2 (y, S u ) dy, où S u est la partie de surface de foie déformée qui doit atteindre la surface observée Γ, et d(y, S u ) est la distance de y à S u . Les notations sont illustrées en figure 2 dans le cas particuler S u = ∂Ω u . Les coefficients M > 0 et α > 0 peuvent être ajustés pour régulariser le champ de déplacement malgré la présence de bruit dans les données.

Chapitre 3. Formes et conditions aux limites mixtes

Pour étudier le problème de contrôle optimal d'un point de vue mathématique, il peut être utile de connaitre certains outils et résultats. Le recalage étant une opération de nature géométrique, les outils que nous présentons dans ce chapitre tournent autour des formes, des distances et projections orthogonales, mais aussi des solutions continues du système de l'élasticité linéaire. On peut définir la régularité d'une forme. Celle-ci correspond à la régularité des transformations qui permettent d'aplatir sa frontière. Nous présentons également la dérivée de forme d'Hadamard, qui donne la possibilité de calculer les dérivées d'une fonction dépendant d'une forme.

Les distances et projections orthogonales sont des outils utiles pour repérer un point par rapport à une forme. La régularité de ces applications dépend en réalité de la régularité de la forme concernée.

Les résultats de continuité pour le système d'élasticité linéaire avec conditions aux bords mixtes s'appuie sur le formalisme des ensembles réguliers au sens de Gröger. Nous présentons cette notion ainsi que certains théorèmes dont nous nous servons par la suite pour obtenir des résultats théoriques.

Chapitre 4. Existence de solutions et conditions d'optimalité

Avant de passer au problème de contrôle optimal proprement dit, nous étudions la fonctionnelle J. Bien qu'elle ne soit pas toujours différentiable au sens de Gateaux, la fonction J admet toujours des dérivées directionnelles, et est continue pour la norme de la convergence uniforme. Sa dérivée dans la direction v s'écrit dJ(u)(v) = ˆΓ min x∈P y (u) [v(x) • (x + u(x) -y)] dy, où les points x de P y (u) vérifient ∥x + u(x) -y∥ = d(y, S u ). Dans le cas où presque tous les points de Γ possèdent une unique projection sur S u , dJ(u) est une forme linéaire et J est différentiable au sens de Gateaux.

Nous étudions l'existence de solutions pour le problème de contrôle optimal en utilisant le formalisme des domaines Gröger-réguliers. La présence de conditions aux bords mixtes représente une difficulté pour montrer l'existence de solutions, car la régularité des déplacements obtenus s'en retrouve affaiblie. Dans ce travail, nous exploitons la régularité W 1,p des solutions de l'élasticité linéaire. Nous proposons un résultat partiel limité à la dimension 2, le cas de la dimension 3 s'avérant nettement plus complexe. Toutefois, quand le système d'élasticité est remplacé par un système elliptique, nous montrons que des solutions existent quelle que soit la dimension.

Nous calculons les conditions d'optimalité de premier ordre pour plusieurs versions du problème de contrôle optimal. Elles font intervenir un état adjoint, qui représente la différentielle de J dans l'espace des contrôles, ainsi que des multiplicateurs de Lagrange pour gérer la contrainte L ∞ sur g. Du fait de cette contrainte, la précaution est de mise dans le calcul des conditions d'optimalité, puisqu'il n'est pas sûr que J soit Gateauxdifférentiable à l'optimum. Le théorème ci-dessous est une version simplifiée du résultat principal de cette section. Il énonce les conditions d'optimalité pour un point critique où J est Gateaux-différentiable. Le résultat présenté dans le manuscrit est plus général et ne fait pas l'hypothèse que J est Gateaux-différentiable. Par ailleurs, nous présentons des exemples où J est différentiable ou pas afin de fixer les idées.

Théorème. Soit g ∈ G M un minimum local du problème de contrôle optimal et u g l'état associé. Supposons en outre que J est différentiable au sens de Gateaux en u g et notons p l'état adjoint correspondant. Alors il existe un multiplicateur de Lagrange λ ∈ L 2 (∂Ω N , R), avec Notre implémentation de la méthode d'adjoint est constituée de trois procédures, chacune ayant un rôle particulier : évaluation de la fonctionnelle J, gestion du modèle d'éléments finis, résolution du problème d'optimisation (voir figure 3). Nous discutons successivement des choix techniques effectués pour chacune des procédures.

En premier lieu, la fonction d'attache aux données J dépend désormais de la position des noeuds du maillage et non plus d'un champ de déplacement continu. L'évaluation de J nécessite d'exécuter une procédure coûteuse pour trouver l'élément de surface du foie le plus proche de chaque point y ∈ Γ. Plusieurs procédures, dont une faisant appel à la distance signée par rapport à l'organe déformé, ont été testées. Notre choix s'est finalement porté sur une méthode dans laquelle les triangles candidats à l'appariement sont rangés dans une structure d'indexation spatiale, permettant ainsi d'effectuer des recherches de plus proche voisin très efficacement. Quant aux dérivées de J, elles peuvent Figure 3 : Structure de la méthode d'adjoint. La chaîne directe (en haut) vise à évaluer la fonctionnelle étant donné un contrôle b, tandis que la chaîne adjointe (en bas) transforme le gradient ∇J(u) en un gradient dans l'espace des contrôles. être évaluées en appliquant directement la formule obtenue dans le chapitre sur l'analyse du problème.

La méthode qu'on utilise pour résoudre le problème d'élasticité direct n'est pas la même selon le modèle élastique choisi. Dans le cas linéaire, la matrice de raideur est assemblée puis factorisée en début d'exécution, puis la forme factorisée est utilisée pour résoudre aussi bien le problème direct que le problème adjoint. Dans le cas non-linéaire, c'est une autre paire de manches. La méthode de Newton est un outil classique pour traiter ce type problème variationnel. Cependant, bien qu'elle affiche un taux de convergence quadratique dans ses dernières itérations, la méthode de Newton sous sa forme traditionnelle n'offre aucun garantie de convergence si le point de départ est situé loin d'une solution. Comme la convergence du solveur direct est une condition sine qua non pour utiliser la méthode d'adjoint, nous décrivons une méthode de Newton à région de confiance, issue du domaine de l'optimisation numérique, pour laquelle stabilité et convergence sont garanties.

Afin de résoudre le problème d'optimisation, nous avons initialement implémenté quelques algorithmes de premier ordre à la main, avant de nous tourner vers des méthodes toutes faites disponibles dans la bibliothèque Scipy, et plus particulièrement une méthode de quasi-Newton à mémoire limitée. Dans le cas où on applique la contrainte L ∞ , cette dernière est gérée de manière primale-duale. Nous illustrons la convergence des différents algorithmes à l'aide de problèmes simples.

• En guise d'application, nous proposons d'utiliser la méthode d'adjoint pour estimer une force ponctuelle appliquée par un instrument.

Les deux derniers points font l'objet de l'article de conférence suivant :

• Guillaume Un autre article, davantage orienté vers l'analyse mathématique du problème, sera très prochainement soumis à une revue internationale à comité de lecture.

Chapter 1

Overview of the manuscript

This project results from discussions between Yannick Privat, professor at the Institut de recherche mathématique avancée in Strasbourg, and Stéphane Cotin, leader of the Inria Mimesis team. The domains of expertise of the Mimesis team include computer-assisted surgery systems, and, more broadly, real-time biomechanical simulation.

Compared to open surgery, minimally-invasive surgery is known to improve the outcome of surgical interventions, reducing pain, bleeding and risk of infection. During a minimally-invasive intervention, the surgeon manipulates organs using instrument inserted through small incisions in the patient's abdomen, and receives visual feedback on a screen, thanks to a laparoscopic camera introduced through one of the incisions (see Figure 1.1). In this context of indirect interaction and feedback, operations such as hepatic tumor resections remain challenging, as the surgeon navigates in an opaque organ, avoiding blood vessels and following the motion due to breathing and heartbeats.

When an augmented reality system is used, an up-to-date view of the liver in its current configuration is displayed to help the medical staff keep track of the position of internal structures, such as tumors and blood vessels. In this manuscript, we are interested in the elastic registration procedure that computes an updated deformation field in the liver volume by combining pre-and intra-operative data. Namely, this procedure aligns a mesh representing the liver as it appears in pre-operative images with a part of the liver surface as it appears on images from the laparoscopic camera during the intervention.

The elastic registration procedure raises several challenges. First, to update the deformed view as the intervention goes on, real-time computation is required. Note that we are not expecting to meet the real-time requirement in this manuscript. In addition, intra-operative data only show a part of the current liver surface, and as a consequence, many solutions to the matching problem exist, and most of them are not consistent with the problem physics. A successful registration approach should reduce as much as possible the range of admissible displacements for a given observed surface, to return a displacement field as close as possible to the true liver displacement.

Chapter 2. State of the art and problem modelling

An elastic registration procedure relies on a mechanical model of the patient's liver, and usually involves the finite element method. In general, much attention is given to the elastic model that describes the liver deformations and interactions with its surroundings (Section 2.2). Existing methods in the literature involve hyperelastic models, and switch to the co-rotational model when real-time computing is required. A wide range of boundary conditions is also used, going from Dirichlet boundary conditions on certain zones of the liver surface to complex elastic model for the ligaments holding the liver.

The registration process itself is often driven by forces applied on the liver mesh to create a deformation (Section 2.3). In a whole class of methods, inspired from the image registration domain, artificial forces simply attract the liver mesh toward the observed point cloud, while the elastic model mostly plays the role of a regularizer. In these methods, we observe that forces driving the registration do not reflect the real causes of displacement. In particular, these forces are created by intra-operative data points, which do not really exist. Other methods are more concerned with reconstructing physically-relevant displacement fields. In this second class of methods, forces and displacements are chosen among an admissible set defined from physical hypotheses, resulting in more realistic displacement fields.

Our approach is closer to the second category. We formulate the registration problem using the optimal control formalism, resulting in a flexible formulation which can be studied and solved numerically using a wide variety of standard optimization tools (Section 2.4). Namely, the optimal control problem studied in this manuscript has the form inf

g∈G M J(u) + α 2 ∥g∥ 2 L 2 (∂Ω N ) subject to constraint      div(Aε(u)) = 0 in Ω 0 u = 0 on ∂Ω Dirichlet Aε(u) • n = g on ∂Ω Neumann .
In this formulation, u is the displacement field generated by the force distribution g. The set of admissible forces is defined by

G M = g ∈ L ∞ (∂Ω N ) ∥g∥ L ∞ (∂Ω N ) ⩽ M .
The optimal control formulation features the data attachment functional

J(u) = 1 2 ˆΓ d 2 (y, S u ) dy,
where S u is the part of the deformed liver boundary that should match the observed surface Γ and d(y, S u ) is the distance from y to S u . Figure 1.2 illustrates the notation in the case S u = ∂Ω u . The L ∞ bound M > 0 and the penalty coefficient α > 0 can be tuned to enforce displacement regularity despite the presence of noise.

Chapter 3. A few tools around shapes and mixed boundary conditions

To obtain mathematical insight about the optimal control problem, some tools and existing results might be of help. As the registration problem is rather geometric, the tools presented in this chapter revolve around shapes, distances, projections, and also making sure that elastic displacements are continuous. It is possible to define the regularity of a shape, based on transformations that should be applied to its boundary to obtain something flat (Section 3.1). We also describe the Hadamard boundary variation framework, which is useful to differentiate a function whose variable is a shape.

Distances and orthogonal projections are useful to locate a point with respect to a shape (Section 3.2). The regularity of those operations depend on the regularity of the concerned shape.

Continuity results for solutions to the linear elasticity system with mixed boundary conditions rely on the framework of Gröger-regular sets (Section 3.3). We describe this notion and cite a few theorems that are of use in the theoretical chapter.

Chapter 4. Existence of solutions and optimality conditions

Before we state some results about the optimal control problem, we study the differentiability of the functional J (Section 4.1). Though it is not always Gateaux differentiable, J always has directional derivatives, and is continuous for the uniform convergence norm. Its derivative in the direction v reads dJ(u)(v) = ˆΓ min x∈P y (u) [v(x) • (x + u(x) -y)] dy, where points x ∈ P y (u) satisfy ∥x + u(x) -y∥ = d(y, S u ). If almost every point in Γ has a single projection onto S u , then dJ(u) is a linear form and J is Gateaux differentiable at u.

We study the existence of solutions to the optimal control problem using the framework of Gröger-regular sets for partial differential equations with mixed boundary conditions (Section 4.2). The presence of mixed boundary conditions represents a difficulty when it comes to proving the existence of solutions, as they degrade the regularity of resulting displacement fields. Our results rely on the uniform W 1,p regularity of solutions to the linear elasticity system. We present a partial result limited to dimension 2, as the three-dimensional case is more complex. However, when the elastic system is replaced by an elliptic system, we obtain existence of solutions in any dimension.

First-order optimality conditions are derived for several versions of the optimal control problem (Section 4.3). They involve an adjoint state, which represents the derivatives of J in the space of controls, as well as Lagrange multipliers to manage the L ∞ constraint on g. Due to the L ∞ constraint, computing optimality conditions requires special attention, as J is not guaranteed to be differentiable at an optimum. The theorem below is a simplified version of the main result of the section. It expresses the optimality conditions for a minimizer where J is Gateaux differentiable. Note that the result presented in the manuscript (Theorem 4.4) is more generic and does not assume Gateaux differentiability. In addition, some examples where J is differentiable or not are given in Section 4.1.

Theorem. Let g ∈ G M a local minimizer of the optimal control problem and u g the associated state. Assume that J is Gateaux differentiable at u j and denote by p the adjoint state. Then there exists a Lagrange multiplier

λ ∈ L 2 (∂Ω N , R) with for a.e. x ∈ ∂Ω N λ(x) = 0 if ∥g(x)∥ < M λ(x) ⩾ 0 if ∥g(x)∥ = M
such that g satisfies the first-order optimality condition for a.e. x ∈ ∂Ω N p(x) + (α + λ(x))g(x) = 0.

Optimization solver Finite element package Functional

New iterate b Solve direct problem

F(u) = b Evaluate J(u) and ∇J(u) Objective gradient ∇Φ(b) = p + ∇R(b) Solve adjoint problem F ′ (u) T p = ∇J(u)

Chapter 5. An adjoint method to solve the registration problem

We discretize the mechanical problem using the finite element method. In the numerical formulation, the optimal control problem reads

min b∈B Φ(b) with Φ(b) = J(u b ) + R(b),
where b is the nodal force distribution at the mesh vertices, and u b is the displacement field defined by the discretized elastic system F(u b ) = b. In the discretized formulation, the only variable seen by the optimization solver is the control b. To compute derivatives of J(u b ) with respect to b, an adjoint method is used (Section 5.2). Namely, the objective gradient reads

∇Φ(b) = p + ∇R(b),
where the adjoint state p is solution to the adjoint problem F ′ (u b )p = ∇J(u b ).

Our implementation of the adjoint method is divided into three parts. The three parts are responsible for evaluating the functional J, managing the finite-element discretization, and solving the optimization problem, respectively (see Figure 1.3). We successively discuss technical options regarding each one of these parts.

First, the discrepancy functional J is now a function of the mesh vertices positions and not of the continuous displacement field (Section 5.3). Evaluating J involves a costly procedure to find the closest triangle in the liver boundary to each point y ∈ Γ. Several procedures were implemented, including one based on the signed distance to the deformed organ. We finally chose a procedure where triangles that are candidates to match Γ are stored in a spatial indexing structure, resulting in efficient nearest-neighbor searches. Concerning the derivatives of J, they can be evaluated using the derivative formula obtained in the theoretical chapter.

The algorithm used to solve the direct and adjoint elasticity problem depends on the chosen elastic model. When a linear model is used, the stiffness matrix is assembled and factorized as the procedure starts and the factorization is used for solving both direct and adjoint problems. However, solving the direct problem in the nonlinear case requires much more precaution. The Newton method is a standard tool for such nonlinear variational problems. Though it enjoys a fast convergence rate in the last iterations, the Newton method in its standard form is not guaranteed to converge when the starting point is too far from a solution. As enforcing convergence of the direct solver is critical for the adjoint method to succeed, we describe a trust-region Newton method, developed by the optimization community, where convergence and stability are guaranteed (Section 5.4).

To solve the optimization problem, we first implemented some artisanal methods, but we finally turned to off-the-shelf algorithms available in the Scipy package (Section 5.5). Namely, we use a limited-memory quasi-Newton method. In the L ∞ -constrained case, the constraint is managed by a primal-dual strategy. We show results on simple toy problems to illustrate the convergence of our procedure.

Chapter 6. Numerical results

The results section aims to illustrate the performance of the adjoint method in terms of registration error. First, we show a few examples with spherical and circular meshes (Section 6.1). In these examples, we illustrate the regularizing effect of the adjoint method, and we show that no guarantee in terms or registration accuracy can be expected from the procedure. We also evaluate the efficiency of regularization approaches involving the L 2 penalty or the L ∞ constraint.

To evaluate the relevance of our approach in augmented surgery, we produce registration results using the Sparse Data Challenge dataset (Section 6.2). The dataset, generated from phantom data, consists of a liver mesh and 112 point clouds representing as many deformed configurations. Registration error statistics are provided by the organizers after we submit the reconstructed deformation computed for each point cloud. We obtained the second-best submission in the challenge, with an average 3.3 mm displacement error, far below the 5 mm required by our clinical collaborators.

In the last result, we try to estimate the net force applied by a robotic instrument on the liver boundary (Section 6.3). As several medical robots are shipped without force feedback, measuring the applied force is useful to avoid causing damages to the patient's liver. We create synthetic deformations and associated point clouds to simulate the liver being manipulated by a robotic instrument. Then we reconstruct a force distribution restricted to a small zone on the liver boundary (see Figure 1.4), and add up nodal forces to obtain a net force estimate. Each one of our five test cases features a time-dependent force (50 frames). We obtain an average force estimation error of 10 %, and an average

Contributions of the thesis

This thesis revolves around four main contributions.

• We identified the need for accurate modeling of the causes of displacement expressed in the literature and proposed an optimal control formulation to respond to this need.

• We proposed a mathematical analysis of the optimal control problem, including existence of solutions and first-order optimality conditions.

• We implemented a modular adjoint method1 to solve registration problems, and validated our implementation using the Sparse Data Challenge dataset.

• We proposed an approach to evaluate a single net force applied by an instrument using our optimal control formulation.

The two last items are the object of the following conference article: Another article revolving around the mathematical analysis will soon be submitted to a peer-reviewed journal.

Chapter 2

State of the art and problem modelling 2.1 Problem overview

A few facts about the liver

Let us begin with a few figures about the liver itself. Located just below the diaphragm, the liver is the heaviest internal organ in the human body. It plays a role in a large number of vital functions, including detoxificating blood, secreting hormones and regulating glycogen storage. Several ligaments hold the liver in place, including the falciform ligament and the left and right triangular ligaments (see Figure 2.1). As a consequence, the liver is not rigidly connected to a bone or a rigid support.

The liver weighs in average 1.5 kg [START_REF] Gray | Anatomy of the human body[END_REF], and its dimensions are about 28 cm × 16 cm × 8 cm [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF]. The matter that composes the liver, called the parenchyma, enters the category of soft tissues. It is nearly incompressible and its Young modulus amounts to approximately 10 kPa [START_REF] Nikolaev | Estimation of boundary conditions for patientspecific liver simulation during augmented surgery[END_REF]. Another consequence of liver tissues containing water is a density of ρ = 1, 070 kg/m 3 [START_REF] Niehues | Liver volume measurement: reason of the difference between in vivo CT-volumetry and intraoperative ex vivo determination and how to cope it[END_REF]. The liver stiffness is increased by the presence of large blood vessels across the parenchyma.

The liver registration problem

A soft organ, the liver is constantly subject to large deformations. During a surgical operation, sources of such deformations include changes in the patient's position, breathing or organ manipulation by the surgeon. During minimally invasive surgery, the surgeon manipulates the liver using tools inserted through small incisions in the patient's abdomen and gets visual feedback on a screen using a laparoscopic camera. In this context, the constant motion and the limited field of view make it difficult for the surgeon to keep track of the position of the liver internal structures. Augmented reality systems have been designed to help the surgeon navigate through the surgical theater.

The objective of augmented reality systems is to track the deformations of the patient's liver during the surgical operation, and provide the medical staff with an updated view of the organ, along with its internal structures such as tumors or blood vessels. In this project, we follow the pipeline proposed by [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF]. In their pipeline, a pre-operative liver geometrical model and intra-operative images from the laparoscopic camera are combined through an elastic registration procedure to produce the updated liver tri-dimensional image.

On the one hand, pre-operative data consist of a liver model extracted from tomographic images. Such images accurately depict the internal structures of the liver. A tetrahedral mesh representing the liver is segmented from images and is embedded with an elastic model. On the other hand, intra-operative data are extracted from laparoscopic images provided in real time by the camera. A point cloud representing the part of the liver surface that is in the camera field of view is created. Figure 2.2 shows a pre-operative liver model, including tumors (green) and blood vessels (blue and red), along with a point cloud superimposed onto a laparoscopic intra-operative image. The augmented reality procedure follows three steps:

1. Intra-operative images are processed to extract a three-dimensional point cloud representing the current liver surface.

2. An initial rigid registration of the pre-operative model onto the point cloud is performed.

3. An elastic registration of the pre-operative biomechanical model is performed to In this manuscript, we are interested in the third part, which is the elastic registration procedure. It consists in computing a displacement field defined on the liver mesh, so that the deformed liver surface matches the intra-operative point cloud. Our purpose is to reconstruct the displacement field we need to apply to the organ to obtain the current configuration (or deformed configuration) from a reference configuration (or initial configuration). Those two configurations are known through pre-operative images and intra-operative measures, respectively.

The elastic registration procedure raises several challenges. First, the augmented reality pipeline is supposed to run in a few tenth of milliseconds in order to provide real-time updates. Though we are not expecting to meet the real-time requirement in this work, we should keep in mind that involved procedures should be as fast as possible. In addition, an intra-operative point cloud is not sufficient to identify a unique displacement field. Therefore, the registration problem has a lot of solutions, and most of them are far from being physically relevant. As a consequence, some other information, such as physical hypotheses, should be added before we start thinking about registration accuracy. The chosen model and methods should result from a trade-off between the physical acceptability of the solution and the efficiency required by real-time execution.

Note that, in this study, we do not consider or compare to registration approaches based on machine learning [START_REF] Brunet | Corotated meshless implicit dynamics for deformable bodies[END_REF][START_REF] Mendizabal | Developments and Novel Approaches in Biomechanics and Metamaterials[END_REF]. Though those modern methods are very efficient, we focus on getting some insight about the problem physics.

Elastic modelling of the liver

A key ingredient to include in a physics-based registration recipe is a good mechanical model. This model describes how the organ is deformed when forces are applied onto it. By providing an accurate description of the liver and its surroundings, one expects their model to behave like the real system in a direct simulation.

In this section, we describe a few models used in the literature to describe the liver. Beforehand, a short reminder about static elasticity problems seems in order. We first introduce the formulation of the elasticity system as a minimization problem, and then we describe a few common elastic models.

Deformation energy and forces

In this section, we give some context about static elasticity problems. We are only interested in static elasticity, which means that inertial effects are not taken into account. The phenomena we want to capture by performing liver registration are the deformations generated by external forces that change along time. Those external forces, whether they are created by the physician's gesture, contact with other organs of breathing movements, evolve very slowly compared to inertial forces. On the other hand, by considering a Young modulus E = 10 kPa, a density ρ = 1070 kg/m 3 and a Poisson ratio ν = 0.49, we obtain that the velocity of mechanical waves in the liver reads

c = E 2ρ(1 + ν) = 1.8 m/s.
As the average liver length is approximately 28 cm, the typical duration of dynamical effects is 0.16 s. Therefore, inertial effects, such as vibrations or mechanical waves propagating across the liver parenchyma, are too fast to be captured by the augmented reality system. An elastic material is a material that returns to its initial shape when left at rest. It is defined by an elastic energy that only depends on the current local deformation of the material (Ciarlet, 1988, Chapter 4). In particular, the elastic energy does not depend on the material velocity, the amplitude of its displacement or any past configuration. We consider a domain Ω ⊂ R 3 filled with an elastic material. The global deformation energy is obtained by accumulating the local energy associated to local deformation over the whole material. Namely, when a displacement field u ∈ H 1 (Ω, R 3 ) is applied to the material, the point x ∈ Ω is transported to x + u(x) and the global elastic energy reads

W (u) = ˆΩ w(∇u) dx,
where w ∈ C 2 (R 3×3 , R) represents the local deformation energy. Here w is a function of ∇u, which means that adding a translation displacement to u does not change the deformation energy. Actually, the deformation energy is insensitive to any rigid displacement, including rotations. We discuss this aspect below when describing usual elastic models. The deformation energy is minimal when u = 0, i.e. the material is in its rest configuration. The energy cost to bring the material from its rest shape to a given displacement u is W (u) -W (0).

External actions that apply onto the solid are represented by forces. Forces are a convenient way to provide a summary of various types of interaction by describing the first-order variation of the total energy around a given displacement u. Examples of forces include gravity forces, that derive from a gravitational potential, or contact forces, that enforce a non-interpenetration condition between two materials. In the context of static elasticity, a force is represented by a linear potential energy. For instance, a surface force distribution g ∈ L 2 (∂Ω, R 3 ) and a volume force distribution f ∈ L 2 (Ω, R 3 ) are represented by the energy terms

-⟨g, u⟩ = - ˆ∂Ω g • u ds and -⟨f, u⟩ = -ˆΩ f • u dx, respectively.
Stable equilibrium configurations correspond to displacement fields where the system energy is at a minimum. In other words, the displacement generated by f and g is a solution to the optimization problem

min u∈U W (u) -ˆΩ f • u dx - ˆ∂Ω g • u ds. (2.1)
In (2.1), U is the set of admissible displacements, or feasible set. Typical feasible sets involve a Dirichlet condition on a subset ∂Ω D ⊂ ∂Ω of the boundary to ensure the existence of a solution. From now on, we assume that a Dirichlet condition is applied, and we use the notation

U = H 1 D (Ω) = u ∈ H 1 (Ω, R 3 ) | u = 0 on ∂Ω D .
The strong formulation of the elasticity system can be derived (at least formally), from the first-order optimality conditions of (2.1). Note that, for v ∈ H 1 D (Ω) a small displacement,

⟨dW (u), v⟩ = ˆΩ w ′ (∇u) : ∇v dx = ˆ∂Ω (w ′ (∇u)n) • v ds -ˆΩ div(w ′ (∇u)) • v dx,
where the last expression is a consequence of Green's formula. Thus, an equilibrium displacement u satisfies the variational equality

∀v ∈ H 1 D (Ω) ˆΩ -div(w ′ (∇u)) -f • v dx + ˆ∂Ω (w ′ (∇u)n -g) • v ds = 0,
which yields the strong formulation

     -div(w ′ (∇u)) = f in Ω w ′ (∇u)n = g on ∂Ω N u = 0 on ∂Ω D , (2.2)
where ∂Ω N = ∂Ω \ ∂Ω D . In the context of elasticity, w ′ (∇u) is called the first Piola-Kirchhoff stress tensor.

Common elastic models

A wide variety of elastic models exists in the literature, with very different properties (see Wex et al., 2015, for examples in biomechanics). Though the elastic problems we consider can be formulated using (2.1) or (2.2), the numerical method used to compute a displacement depends on the model complexity. In this section, we introduce three of the most common elastic models. In those models, the material resistance to deformation and volume changes are characterized by the Lamé parameters, which are deduced from the Young modulus E and the Poisson ratio ν by

µ = E 2(1 + ν) and λ = Eν (1 + ν)(1 -2ν)
.

(2.3)

Saint Venant-Kirchhoff model

The simplest nonlinear elastic model is the Saint Venant-Kirchhoff model (St Venant, 1844;[START_REF] Kirchhoff | Über die Gleichungen des Gleichgewichtes eines elastischen Körpers bei nicht unendlich kleinen Verschiebungen seiner Theile[END_REF] which applies a quadratic penalty on the nonrigid part of the displacement.

Let us introduce a few notations. When a displacement u is applied to the elastic solid, a particle initially located at x ∈ Ω is transported to x + u(x). The gradient of Id +u, denoted by F = (I + ∇u), describes the local transformation around x, including nonrigid deformation and rotation. The tensor that is penalized in the Saint Venant-Kirchhoff energy is the Green-Saint Venant strain tensor

e(u) = F T F -I 2 = 1 2 ∇u + ∇u T + ∇u T ∇u .
Note that the rotation component of F is completely evacuated when F T F is considered, and thus only nonrigid deformations such as shear and volume changes are seen by the Green-Saint Venant tensor.

The local deformation energy is characterized by Hooke's tensor A, defined by Aε = 2µε + λ tr(ε)I, with µ and λ from (2.3). Its expression reads

w(∇u) = 1 2 Ae(u) : e(u) = µ|e(u)| 2 + λ 2 (tr e(u)) 2 .
Though it is a nonlinear model, the Saint Venant-Kirchhoff model involves a simple energy function and is not relevant for large deformations. In particular, it allows the deformation gradient F to be an indirect transformation, which translates as inverted triangles in a finite element simulation. On the other hand, solving a static problem requires to use a Newton method, which comes at a high computational cost. For these reasons, the Saint Venant-Kirchhoff model is not relevant for applications where fast computations are required or large deformations are involved.

Linear elasticity

The linear elastic model is obtained by performing a quadratic approximation of the Saint Venant-Kirchhoff energy around u = 0. Namely, the Green-Saint Venant tensor e(u) is replaced by the linearized strain tensor ε(u) = 1 2 ∇u + ∇u T and the linear elastic energy reads

w(∇u) = 1 2 Aε(u) : ε(u) = µ|ε(u)| 2 + λ 2 div(u) 2 .
When a Dirichlet condition is involved, the existence and uniqueness of the solution for smooth domains is a consequence of the Lax-Milgram theorem and Korn's inequality (see a detailed explanation in Allaire, 2007, Section 2.3). The solution of (2.1) solves the linear partial differential equation

     -div(Aε(u)) = f in Ω 0 Aε(u) • n = g on ∂Ω N u = 0 on ∂Ω D .
From a numerical point of view, the linear elastic model allows very fast computations, since a displacement is computed by solving a positive definite sparse linear system. Methods to solve a positive definite linear system include the Conjugate Gradient method [START_REF] Hestenes | Methods of Conjugate Gradients for Solving Linear Systems[END_REF] or direct methods involving the Cholesky factorization (see Nocedal and Wright, 2006, Appendix A). Direct methods are especially relevant when the elasticity system is solved several times, as the Cholesky factorization can be computed once and stored as long as needed.

Due to its definition, the linear elastic model only provides a first-order approximation of the displacement with respect to the force distribution. As a consequence it is only relevant for small displacements (which is a very subjective notion), typically when the material is very stiff compared to the forces applied onto it.

Neo-Hookean model

The Neo-Hookean model [START_REF] Rivlin | Large elastic deformations of isotropic materials IV. further developments of the general theory[END_REF] is more complex than the Saint Venant-Kirchhoff model, but it is far more relevant to describe soft materials undergoing large deformations. The local deformation energy penalizes the Green-Saint Venant deformation tensor e(u) and the Jacobian determinant J = det(F ). Though several formulations exist [START_REF] Smith | Stable Neo-Hookean Flesh Simulation[END_REF], we stick to a version previously used in the Mimesis team and software [START_REF] Brunet | Exploring new numerical methods for the simulation of soft tissue deformations in surgery assistance[END_REF],

w(∇u) = µ tr(e(u)) -µ ln(J) + λ 2 (ln J) 2 .
A key feature of Neo-Hookean models is the presence of terms involving ln(J) in the expression. These terms are defined only when the local deformation F is direct (i.e. det(F ) > 0) and grow toward infinity as F becomes more and more singular. As a consequence, resulting displacement fields only involve direct local deformations, and deformed meshes do not exhibit inverted triangles. In Figure 2.3, we apply a uniform volume force distribution f (x) ∝ (-4, 1) onto a two-dimensional rectangular beam, and we compute the resulting deformation using the three models described in this section. The color map represents the local Jacobian determinant of the transformation. Namely, inverted triangles are represented in red while non-inverted triangles are represented in blue. First, the linear deformation, which is just an amplification of an infinitesimal deformation, features a significant dilatation, which is common with linear deformations. In addition, rotations are penalized by the linear model, as the Green-Saint Venant tensor e(u) is replaced with its linearized version. As a consequence, the resulting displacement does not involve a sufficient rotation. Rotation is better represented in the Saint Venant-Kirchhoff displacement, especially at the end of the beam. However, the displacement also features inverted triangles close to the clamped boundary. Finally, the Neo-Hookean displacement looks more realistic as it features a large rotation, and contains no inverted triangles. Even though some overlapping appears between triangles, every triangle is transformed in a direct way. Thus, the Neo-Hookean model is relevant for large deformations without self-contact.

Solving numerically a neo-Hookean problem using a Newton method might be challenging, as the problem residual is not defined for non-feasible points (J < 0). As a consequence, the solver should only generate iterates that lie in the feasible region.

Remark 2.1. Let us add a word about the existence of solutions for nonlinear models. The neo-Hookean model, along with some more generic models, falls into the category of polyconvex stored energy functions [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. For such functions, the existence of a minimizer is guaranteed. In the case of Saint Venant-Kirchhoff materials, [START_REF] Raoult | Non-polyconvexity of the stored energy function of a Saint Venant-Kirchhoff material[END_REF] proves that the deformation energy function is not polyconvex, and proving the existence of minimizers is more difficult (see Ciarlet, 1988, Section 7.7 Remark (5) and Exercise 7.16). Though, it is possible to prove that the strong formulation (2.2) has solutions, using the implicit functions theorem (Ciarlet, 1988, Chapter 6).

Mechanical models in liver simulation

We now review a few models that were proposed in the computer-assisted intervention literature to model the liver and its surroundings. Such models should not only involve a constitutive law to describe the liver parenchyma, but also account for the interactions between the liver and its surroundings.

We first focus on constitutive laws used to model the liver parenchyma. Liver tissues, as well as most biological tissues, enter the category of soft materials and their simulation requires a hyperelastic model. In their extensive review [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF] indicate that, in general, the liver parenchyma exhibits an isotropic behavior, and thus can be represented using an isotropic model. Used models in the literature include the classic Saint Venant-Kirchhoff [START_REF] Delingette | Soft Tissue Modeling for Surgery Simulation[END_REF] and Neo-Hookean [START_REF] Miller | Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF] models, but also the more generic Ogden model [START_REF] Nikolaev | Estimation of boundary conditions for patientspecific liver simulation during augmented surgery[END_REF], defined by

w(∇u) = N p=1 µ p α p λ α p 1 + λ α p 2 + λ α p 3 -3 ,
where λ 1 , λ 2 , λ 3 are the principal stretches and N, α p , µ p are material parameters. However, hyperelastic models often require heavy calculations, too heavy for real-time registration. A popular compromise between model accuracy and computational efficiency, the linear co-rotational model [START_REF] Nesme | Efficient, Physically Plausible Finite Elements[END_REF] offers the possibility to simulate nonlinear deformations while solving linear systems. Adopters of the co-rotational model for liver registration include [START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF], [START_REF] Plantefève | Patient-Specific Biomechanical Modeling for Guidance During Minimally-Invasive Hepatic Surgery[END_REF] and [START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF], among others. The liver parenchyma is traversed by large blood vessels that are stiffer than liver tissues, resulting in an additional resistance to deformations. [START_REF] Peterlík | Modeling and Real-Time Simulation of a Vascularized Liver Tissue[END_REF] propose a liver biomechanical model including the mechanical properties of blood vessels. To meet the real-time requirement, the authors use the co-rotational model for the parenchyma, coupled with a beam model for blood vessels.

Another critical part of the mechanical model is the interaction between the liver and its surroundings. Though it is feasible to perform ex-vivo measures on the parenchyma to estimate mechanical parameters, estimating the parameters of ligaments and tissues that hold the liver is more difficult. As a consequence, boundary conditions applied to the liver mesh are often simple. To produce their synthetic validation data, [START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF] apply homogeneous Dirichlet conditions on the parts of the boundary attached to ligaments and at the entry of the portal and the hepatic veins. [START_REF] Ozkan | Compliance boundary conditions for patient-specific deformation simulation using the finite element method[END_REF] consider that the liver is embedded in a larger elastic continuum, and set the boundary conditions as an approximation of the behavior of surrounding tissues. [START_REF] Peterlík | Image-Driven Stochastic Identification of Boundary Conditions for Predictive Simulation[END_REF] represent ligaments surrounding the liver as linear springs holding the organ. They use a data assimilation procedure involving a Kalman filter to progressively estimate the associated stiffnesses. Their model was recently improved by [START_REF] Nikolaev | Estimation of boundary conditions for patientspecific liver simulation during augmented surgery[END_REF], who use nonlinear springs and a better initialization procedure.

A review of some registration methods

We now turn to elastic registration methods. Based on a mechanical model of the liver, these methods apply forces onto the liver virtual model to generate a deformation. The choice of a registration method is as critical as the choice of a mechanical model, as the forces applied by the method determine whether the resulting deformation is physically relevant or not.

Before we focus on registration methods for augmented surgery, we give a bit of context about image registration methods, as those methods have been around for several decades and have probably inspired the augmented surgery community.

Image registration methods based on a mechanical model

Image registration is a vast and well established domain of image processing. In the medical domain, image registration is used to account for anatomy variability between patients [START_REF] Christensen | Volumetric transformation of brain anatomy[END_REF], to combine images from different modalities [START_REF] D'agostino | A viscous fluid model for multimodal non-rigid image registration using mutual information[END_REF], or to match data with an atlas [START_REF] Wang | Physical model-based non-rigid registration incorporating statistical shape information[END_REF]. In this section, we describe some registration methods that involve a mechanical model, namely a solid or a fluid model. The review by [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] adopts a wider scope and includes methods based on interpolation and geometry.

Let us consider an image registration problem involving two liver images. Each image is represented by a density function (1 for bone, 0 for void) with compact support defined on R d , where d = 2 or 3 is the image dimension. The template image, represented by the distribution t, is for instance an atlas where each part of the liver is identified. The source image, represented by the distribution s, may be a patient-specific image. During the registration process, a displacement field u is computed, so that the point x in the template image is associated to the point x + u(x) in the source image. In other words, the template image is transported by a displacement field u to look like the source image. Usually, the displacement field u is the solution to a minimization problem, where the cost function F measures the discrepancy between t and s • (Id +u). A common choice for F (see for instance [START_REF] Rabbitt | Mapping of hyperelastic deformable templates using the finite element method[END_REF][START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF] is the least-squares error

F (u) = 1 2 ˆRd t(x) -s(x + u(x)) 2 dx.
While the choice of F depends on the available data and image modalities, the chosen physical model is related to the expected qualities of u, especially its smoothness and the possibility to obtain large deformations.

Models based on elasticity

In an elastic registration process, we consider that the domain is filled with an elastic material, so that applying a displacement field has a cost in terms of deformation energy.

The functional F defined in (2.3.1) plays the role of an external potential energy, which applies a loading onto the material to generate a deformation. The resulting displacement u is solution to the problem min

u F (u) + W (u), (2.4)
where W is the elastic deformation energy, which plays the role of a regularization term. First mentions of elastic registration appear in Broit's doctoral thesis (1981), followed by a series of articles involving his advisor where elastic registration is applied to match brain images with an atlas [START_REF] Bajcsy | Multiresolution elastic matching[END_REF][START_REF] Gee | Elastic matching: Continuum mechanical and probabilistic analysis[END_REF]. While these early works involve a linear elastic model, hyperelastic models are later adopted [START_REF] Rabbitt | Mapping of hyperelastic deformable templates using the finite element method[END_REF] to account for a wider range of deformations. Subsequent developments around hyperelastic models include struggling against inconsistent deformations, computing the reverse transformation along with the direct one, or writing more efficient numerical schemes.

Models based on viscosity

As they penalize the global deformation applied to the template image, elastic models may be too restrictive when large deformations are required. [START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF] prefer to use a model based on fluid mechanics. In a viscous model, the deformation energy is applied to the velocity field instead of the displacement field. As a consequence, the regularity of the final displacement field is preserved, but nothing prevents large deformations to occur.

The resulting registration is the consequence of an time-dependent displacement field u : [0, t f ] × Ω → R d . A particle initially at x is transported at x + u(t, x) at time t. The evolution of u is driven in a Eulerian fashion by a velocity field v : [0,

t f ] × Ω → R d such that ∀(t, x) ∈ [0, t f ] × Ω ∂ t u(t, x) = v(t, x + u(x, t)).
To enforce the smoothness of the final displacement u(t f , •), v is the solution to the elasticity problem min

v dF (u) • v + W (u).
Several modern shape-matching methods rely on velocity-based models. In the approach by de [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF], the authors use the shape optimization framework to register an unstructured mesh onto a fixed shape. Following the extension-regularization framework, they define the velocity field at each iteration as the solution to an elasticity problem. In the large deformation diffeomorphic metric mapping (LDDMM) framework [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF][START_REF] Beg | Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms[END_REF], still used nowadays, the velocity field is solution to a variational problem involving the whole time interval.

Organ registration in augmented surgery

Augmented surgery represents a particular application of registration, as pre-and intraoperative data both represent the same organ. In this context, the mechanical model has a real meaning, as the computed deformation should represent the true organ deformation. In contrast, mechanical models only play a role of regularizer in the image registration methods mentioned above.

Forces resulting from a fictitious agent

A large number of organ registration procedures are inspired from the famous Iterative Closest Point algorithm by [START_REF] Besl | Method for registration of 3-D shapes[END_REF]. Created for generic rigid registration, the Iterative Closest Point algorithm relies on orthogonal projections between the intra-operative point cloud Γ = {y 1 , y 2 , • • • } and the current liver mesh. At each iteration, the orthogonal projections a 1 , a 2 , • • • of points in Γ onto the liver mesh boundary are computed. Then, a rigid displacement u r , solution to the least-square problem

min u r 1 2 j ∥a j + u r (a j ) -y j ∥ 2 ,
is determined. This rigid displacement is applied to the liver mesh, projections onto the new mesh are computed and so on. The method stops when the least-square error does not decrease sufficiently between iterations. [START_REF] Clements | Robust surface registration using salient anatomical features for image-guided liver surgery: Algorithm and validation[END_REF] propose an improved version of the Iterative Closest Point algorithm, where landmark correspondence between the pre-and intra-operative data is taken into account. Adaptations of the Iterative Closest Point algorithm to elastic registration often consist in introducting fictitious forces between the observed point cloud and the liver model. For [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF], these forces are generated by linear springs (see Figure 2.4). Using the same notation as above, the static elasticity problem solved at each iteration reads

min u 1 2 j k j ∥a j + u(a j ) -y j ∥ 2 + W (u),
where k j is a stiffness parameter, which increases progressively along iterations to enforce the matching. Note how close to (2.4) this formulation is. [START_REF] Plantefève | Patient-Specific Biomechanical Modeling for Guidance During Minimally-Invasive Hepatic Surgery[END_REF] replace the linear springs with a nonlinear attractive force. They model the registration process using dynamic elasticity to ensure the existence of a solution without Dirichlet boundary conditions. In a similar fashion, [START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF] use an electrostatic attractive force between positively charged observed points and a negatively charged liver boundary. Artificial forces sometimes appear in a more implicit way, taking the form of an optimization constraint in the mechanical model. To perform the registration, [START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF] solve the optimization problem min u W (u) subject to constraint C(u) = 0, (2.5)

Figure 2.4: Observed data points (green) and their orthogonal projection onto the liver mesh [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF].

where the constraint C(u) = 0 enforces correspondence between the current displacement field u and observed data. In an earlier work [START_REF] Peterlík | Model-Based Identification of Anatomical Boundary Conditions in Living Tissues[END_REF], the authors use the resulting Lagrange multipliers to estimate the boundary condition necessary to perform a registration between two configurations of the same organ. Approaches using similar sliding constraints are used in the context of pelvic system simulation [START_REF] Courtecuisse | Three-dimensional physics-based registration of pelvic system using 2D dynamic magnetic resonance imaging slices[END_REF] and brain-shift compensation [START_REF] Morin | Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation[END_REF].

Forces resulting from an inverse problem

Whether they derive from an imaginary potential energy or play the role of constraintenforcing Lagrange multipliers, such fictitious forces do not correspond to a real physical phenomenon. They are created by data points, which do not really exist. As a consequence, in the methods above, forces that drive the registration are necessarily applied where data points are, without consideration for the true causes of displacement. If the registration problem involved landmarks located inside the organ, then registration forces would not even apply on the liver boundary. This often results in a poorly physically consistent displacement field, regardless of the direct model accuracy. [START_REF] Rucker | A Mechanics-Based Nonrigid Registration Method for Liver Surgery Using Sparse Intraoperative Data[END_REF] address this issue in the context of open liver surgery. They control the imposed displacement on the posterior face of the liver (in contact with the support), while the anterior surface is observed by the camera. [START_REF] Heiselman | Characterization and correction of intraoperative soft tissue deformation in image-guided laparoscopic liver surgery[END_REF] extend this work to laparoscopic surgery, by applying an imposed displacement on zones where ligaments hold the liver. Here, the authors use linear elasticity. In order to make computations faster, they create a basis of displacement fields associated to a range of non-homogeneous Dirichlet conditions. The authors use a similar displacement basis in a more recent work [START_REF] Heiselman | Intraoperative Correction of Liver Deformation Using Sparse Surface and Vascular Features via Linearized Iterative Boundary Reconstruction[END_REF] to register the liver model onto slices obtained from ultrasound intra-operative measures. Özgür et al. (2018) propose a method for initial pose reconstruction when depth information is not available in laparoscopic images. They compute a rest shape by compensating the effect of gravity in pre-operative images, and they simulate pneumoperitoneum pressure and gravity load in intra-operative conditions to determine the final liver shape. [START_REF] Xie | Constrained finite element method for runtime modeling of soft tissue deformation[END_REF] use surface-matching as a constraint in an optimization problem to update a force estimation between two successive frames. From a current force distribution g n , they estimate a new force distribution g n+1 by solving the problem min

g 1 2 ∥g -g n ∥ 2 subject to constraint C(u g ) = 0,
where C is defined as in (2.5) and u g is the displacement field generated by g. Note that, here, the constraint is not imposed in the direct problem, but in the inverse problem, so that no artificial protagonist is introduced into the direct problem.

Our contribution: an optimal control formulation

The approach studied in this manuscript falls in the category of registration methods based on an inverse problem. Methods described in the paragraph above lead to very accurate registration (see for instance the Vanderbilt result in Figure 6.8). However, these methods are often very integrated, specific to a given problem, and formulated in terms of matrix-vector operations. In addition, certain features, such as the linear displacement basis, are tailored for linear elasticity.

In this manuscript, we take a step backward and formulate the registration problem in the generic optimal control framework. Using an optimal control problem results in a very flexible formulation that can be easily modified to include new physical information or intra-operative data. The optimal control framework offers a wide range of generic tools to study and solve the problem.

Our contribution includes a mathematical study of the continuous registration problem and an implementation of an adjoint method to solve the problem numerically, even with a nonlinear elastic model. We also mention a possible application of our registration procedure to estimate a local force applied onto the liver.

In the remaining of this section, we introduce the notation we use all along the manuscript, and we introduce the optimal control formulation.

Mechanical model and data representation

The organ and its deformed configurations

We first define the notation for the liver constitutive law. Figure 2.5 illustrates the problem geometry. In its reference configuration, the liver is represented by an open domain Ω 0 ⊂ R d (with d = 3) filled with an elastic material. Due to its interactions with the environment, this continuum is subject to surface loadings, such as pressure or contact with surrounding tissues. As we neglect gravity, we do not consider body forces, though they may be easily accounted for as a parameter of the problem. The organ boundary ∂Ω 0 splits into two parts, ∂Ω D and ∂Ω N . On ∂Ω D , a homogeneous Dirichet boundary condition applies, whereas ∂Ω N is subject to a Neumann condition involving a surface force distribution g ∈ L 2 (∂Ω N , R d ). In our numerical experiments, Dirichlet boundary conditions are set at the hepatic vein entry and in the region where the liver is in contact with the falciform ligament. The space of displacement fields is defined by

Ω u = (Id +u)(Ω 0 ) Ω 0 Γ y p ∂Ω u (y) d ( y , ∂ Ω u )
H 1 D (Ω 0 ) = v ∈ H 1 (Ω 0 ) v = 0 on ∂Ω D ,
and we denote by u g ∈ H 1 D (Ω 0 , R d ) the displacement field generated by a given force distribution g ∈ L 2 (∂Ω N ). In the case of linear elasticity, u g solves the partial differential equation

     div(Aε(u)) = 0 in Ω 0 u = 0 on ∂Ω D Aε(u) • n = g on ∂Ω N .
(2.6)

Unless specified otherwise, we use a linear elastic model to describe the liver parenchyma. In addition, we do not model blood vessels in this study.

We denote by Ω u = {x + u(x), x ∈ Ω 0 } the volume occupied by the deformed organ. To be consistent with the numerical implementation, we do not expect (Id +u) to be a bijection. In the numerical framework, the liver is represented by a mesh and its boundary is represented by a set of triangles, but we do not take measures to detect or prevent self-intersection. When a boundary triangle is inside the deformed mesh, it is still considered to be a boundary element. For this reason, in the whole manuscript we always work with the surface (Id +u)(∂Ω 0 ), which is the boundary ∂Ω 0 transported by a displacement u, even when it does not match with the boundary of Ω u . However, the range of deformations expected in the application case usually does not include self-contact. In that case, both definitions coincide.

Data representation and objective function

While the reference configuration Ω 0 is known through pre-operative images, the current position shall be estimated with the help of data provided during the intervention by a measuring device. In the continuous problem, we represent observed data as a (d -1)dimensional manifold Γ (a surface, actually). This data represents the current location of a part of ∂Ω 0 .

Due to the liver geometry, it is very likely that the surgeon has a vague idea of which part of the liver surface appears on the camera, and therefore which part of the pre-operative liver boundary should match Γ. The part of the liver boundary that is candidate to match Γ is denoted by S 0 ⊂ ∂Ω 0 . From now on, we assume that S 0 is compact. When S 0 is transported by a displacement u, we use the notation

S u = (Id +u)(S 0 ) ⊂ (Id +u)(∂Ω 0 ).
In this context, the estimated displacement u should satisfy Γ ⊂ S u .

(2.7)

To enforce this condition, we compute u as the minimizer of a functional J : C(Ω 0 ) → R, which vanishes if and only if (2.7) holds. Moreover, it should not promote one part of S 0 in particular to match Γ. The design of our objective function is inspired by the work by de Buhan et al. ( 2016). In their procedure, a meshed set Ω is transported along a velocity field to match a target shape Ω T , known only through its signed distance function d Ω T , defined by

∀y ∈ R d , d Ω T (y) =      d(y, ∂Ω T ) if x / ∈ Ω T 0 if x ∈ ∂Ω T -d(y, ∂Ω T ) if x ∈ Ω T .
The proposed objective function reads

J(Ω) = ˆΩ d Ω T (x) dx.
(2.8)

It is minimized by using a gradient descent. The authors also considered the functional

P (Ω) = ˆ∂Ω d 2 Ω T (x) ds
to measure the discrepancy between the boundaries ∂Ω and ∂Ω T [START_REF] Nardoni | Mesh deformation strategies in shape optimization[END_REF]. However, this functional led to convergence difficulties, and for this reason they preferred to use (2.8). Unlike de Buhan et al., we are not matching a full shape, but only an open surface, and we cannot use a signed distance function. Moreover we cannot penalize the distance to Γ on the whole boundary ∂Ω 0 , as only an unknown subpart of it should match Γ.

Finally, we use a least-squares functional similar to that in previous registration methods. Inspired from the Iterative Closest Point algorithm, the functional J reads

J(u) = 1 2 ˆΓ d 2 (y, S u ) dy.
(2.9)

Note that, as d 2 (y, S u ) = min x∈S u ∥x -y∥ 2 , J(u) vanishes if and only if Γ ⊂ S u up to a zero Lebesgue measure set. The expression (2.9) does not seem to promote a particular part of S 0 as a better candidate to match Γ, at least among feasible solutions u such that J(u) = 0. Though, evaluating the distance function involves orthogonal projections, which promotes the part of S u that is closest to Γ.

In an improved version of the functional, pairings between points from Γ and points from S u may be performed using something else than orthogonal projections. Other surface-matching techniques are studied in the registration community, including functional maps [START_REF] Ovsjanikov | Functional Maps: A Flexible Representation of Maps between Shapes[END_REF] or varifolds [START_REF] Antonsanti | Partial Matching in the Space of Varifolds[END_REF]. For a review of surface-matching methods, see [START_REF] Sahillioğlu | Recent advances in shape correspondence[END_REF].

Remark 2.2. A possible minor improvement in the expression for J consists in giving a weight to points of Γ, depending on their measure uncertainty. Taking in account different noise variances for the points in Γ could lead to an objective function of the form

J σ (Ω) = 1 2 ˆΓ d 2 (y, S u ) σ 2 (y) dy, with 1/σ 2 ∈ L ∞ (Γ).

An optimal control formulation

While the role of J is to associate points from Γ and S u to compute a discrepancy, the role of the optimal control formulation is to make sure that the chosen minimizer of J satisfies the physical model. A simple version of the optimal control problem reads

inf g∈G J(u g ), (2.10)
where G = L 2 (∂Ω N ). The formulation (2.10) contains the minimal information necessary to compute a displacement field, and additional information should be added to reduce the set of solutions. For instance, new data, such as anatomical landmarks may be taken into account in the objective function, while information about the nature of forces that create the displacement should be included in the definition of G. In the following paragraph, we consider two options to promote regular displacement fields despite the presence of noise in observed data. At the end of the section, we mention another formulation with a constraint on the displacement field.

Penalty on the control norm

If we solve (2.10), the algorithm finds a control such that almost every point of Γ is included in S u g . In the medical application, Γ is extracted from laparoscopic images using depth estimation. As a consequence, the intra-operative surface observation is polluted by noise. An exact matching is not what we look for, as it would result in a very irregular displacement field. A better result would be to make the liver boundary pass through the point cloud without matching every single point. To enforce control smoothness, we consider two options, namely a regularization term and a constraint.

In both approaches, the goal is to prevent the control g from taking the large values it needs to reach every point in Γ.

In the regularized approach, a penalty term is included in the objective function, to prevent the control from oscillating too much. The regularized problem reads

inf g∈G J(u g ) + R(g), (2.11)
where R(g) is a penalization term, for instance

R(g) = α 2 ˆ∂Ω N ∥g∥ 2 ds or R(g) = α 2 ˆ∂Ω N ∥∇ ∂Ω N g∥ 2 ds,
and G is L 2 (∂Ω N ) or H 1 (∂Ω N ) depending on the chosen penalty. Though solving the regularized problem should result in a smoother control, it has a poor physical meaning, as the penalty is non-local. The algorithm could output a control with oscillations around Γ to save some cost on J, while creating an artificially smooth control on other parts of the boundary to save some cost on the penalty. Moreover, it is difficult to set a coefficient α based on physical considerations.

In the constrained approach, the feasible set is defined by

G M = g ∈ L ∞ (∂Ω N ) | ∥g∥ L ∞ (∂Ω N ) ⩽ M .
(2.12)

Admissible force distributions should not take pointwise intensities greater than a upper bound M > 0. Though managing a pointwise constraint is more difficult than evaluating a penalty term, it seems easier to find a parameter M with a physical meaning. Namely, M may be chosen so that pressures greater than M cannot be expected in the human body.

For future reference, we write here the problem that is considered in Chapter 4. The problem, which involves both a L 2 penalty term an a L ∞ constraint, reads

inf g∈G M J(u g ) + α 2 ∥g∥ 2 L 2 (∂Ω N ) .
(2.13)

Movement restriction on a subdomain

Information added into the optimal control formulation may also result from statistical observations. Assume that, in general, it is observed that the displacement field in a subdomain of the liver does not exceed a certain magnitude. Namely, a subdomain ω 0 ⊂ Ω 0 is expected to remain close to its initial configuration during the surgical intervention. This information may be included in the problem by adding a constraint on the displacement norm on ω 0 . The resulting problem reads

inf g∈G M J(u g ) subject to constraint 1 2 ˆω0 ∥u g ∥ 2 dx ⩽ U, (2.14)
where U > 0 is a parameter. Here, U can be estimated by measuring the average displacement amplitude in ω 0 during past experiments. Note that optimality conditions for problem (2.14) are derived in Chapter 4, but this approach is not implemented numerically, as we obtained satisfying results with the other approaches. However, it remains an interesting alternative to introduce physical information into the optimization problem.

Chapter 3

A few tools around shapes and mixed boundary conditions

One of our major contributions in this manuscript is to tackle the liver registration problem using a mathematical formalism. In this chapter, we introduce a few mathematical tools that should be of use to analyse the optimal control problem. First, as the initial momentum of the project revolved around shape optimization, we give some definitions around shapes, shape functionals and shape optimization. The mathematical definition of a shape is not really different from what anyone would call a shape. A shape functional is a functional whose parameter is a shape, and shape optimization consists in minimizing a shape functional. However, shape optimization requires some dedicated tools, as inner products and other features from Hilbert spaces are not available in a set of shapes.

Then, to obtain insight about the discrepancy functional J and approaches involving closest points, we state some results around the signed distance function and orthogonal projections. Regularity properties of J are very similar to those of the signed distance function to a shape. They both depend on the uniqueness of an orthogonal projection of a point onto a shape. In addition, those two functions share the properties of applications defined using the min operator.

Finally, we introduce the notion of Gröger-regular set, which is useful to study partial differential equations with mixed boundary conditions. The Gröger-regular property requires sufficient regularity for the domain boundary, but also for the boundary between regions with Dirichlet boundary conditions and regions with Neumann boundary conditions. Using this framework, we state a few regularity results found in the literature, for the elasticity system, but also for elliptic systems. The results mentioned in this section are used to study the existence of a solution to our optimal control problem.

Shapes and shape functionals

Surface registration is a very geometric operation, and in this work we need to manipulate functionals that depend on a deformed shape. In addition, the regularity of solutions to

x 0 B ψ(x 0 ) ψ(B) ψ ∈ C k (B, R d ) Figure 3.1: Local C k regularity for a shape. The C k -diffeomorphism ψ transforms the boundary part ∂Ω ∩ B into a part of the plane of equation x d = 0.
partial differential equations often depends on the domain regularity. For this, reason, we provide in this section some facts about shapes.

In this section, Ω is a domain, which means a bounded open subset of R d , and ∂Ω is its boundary. We might also use the term shape to emphasize the status of Ω as the variable of a shape functional, but it actually means the same as a domain. We first give a few definitions about shape regularity. Then, we introduce functionals that depend on a shape variable, and a tool to study their variations, Hadamard's boundary variation framework.

Regularity of a domain

In this section, we begin with some tools to describe a domain and obtain local information about its boundary. Then we briefly evoke the transposition of these notions in a numerical context.

Domains of class C k

The local regularity of a shape Ω is assessed by considering the part of its boundary ∂Ω that lies in the vicinity of any point x ∈ ∂Ω and answering the following question: how regular is a diffeomorphism which transforms this part of the boundary into something flat, namely a part of an affine subspace of dimension d -1? When, for instance, a diffeomorphism of class C k can be found for each point x of ∂Ω, then Ω is said to have a C k boundary (see Figure 3.1). A more formal definition (Gilbarg and Trudinger, 1977, Section 6.2) is given below.

Definition 3.1 (Domain of class C k ). Let 0 ⩽ k ⩽ +∞. A bounded domain Ω ⊂ R d and its boundary are of class C k when, at each point x ∈ ∂Ω, there is a ball B = B(x 0 ) and a C k -diffeomorphism ψ from B to D ⊂ R d such that ψ(B ∩ Ω) ⊂ R d + and ψ(B ∩ ∂Ω) ⊂ ∂R d + ,
where ,ℓ are defined the same way.

R d + = x ∈ R d | x d > 0 . For 0 < ℓ ⩽ 1, domains of class C k
The regularity of a domain Ω is useful to ensure regularity of solutions to partial differential equations defined on Ω, but also to study the variations of a shape functional when Ω is deformed. In the second situation, tools such as boundary normal or curvature sometimes need to be available to compute shape derivatives. These tools are defined right now.

We first define the boundary normal, which is characterized using the diffeomorphism ψ and its jacobian matrix [START_REF] Delfour | Advances in Design and Control. Metrics, analysis, differential calculus, and optimization[END_REF], Chapter 2, Remark 3.2).

n(x 0 ) = m(x 0 ) ∥m(x 0 )∥ where m(x 0 ) = -∇ψ(x 0 ) T e d ,
and ∇ψ is the Jacobian matrix of ψ.

Remark 3.1. To be convinced by the expression of m(x 0 ), note that the tangent space to ∂Ω at x 0 is defined by T x 0 = span(∇ϕ(0)e 1 , . . . , ∇ϕ(0)e d- 1), where ϕ = ψ -1 . Therefore we require m(x 0 ) to satisfy

∀1 ⩽ i ⩽ d -1 m(x 0 ) • ∇ϕ(0)e i = 0 m(x 0 ) • ∇ϕ(0)e d = -1.
In other words,

∇ϕ(0) T m(x 0 ) = -e d and m(x 0 ) = -∇ϕ(0) -T e d = -∇ψ(x 0 ) T e d .
To define curvature, it is easier to use the representation of C k domains as local epigraphs, introduced below (Delfour and Zolésio, 2011, Chapter 2, Definition 5.2 and Theorem 5.2).

Definition 3.3 (Local epigraph). A domain Ω

⊂ R d with nonempty boundary is said to be locally a C k -epigraph, k ⩾ 0, when for each x 0 ∈ ∂Ω, there is a ball B centered at 0, a hyperplane S ⊂ R d with normal a and a C k function φ : S ∩ B → R such that ∀v ∈ B x 0 + v ∈ Ω ⇔ v • a > φ(v S ) and x 0 + v ∈ ∂Ω ⇔ v • a = φ(v S ), where v S is the component of v along S. Proposition 3.1. For k ⩾ 1, a domain Ω has a C k boundary if and only if it is locally a C k -epigraph.
For k ⩾ 1 and 0 < ℓ ⩽ 1, domains of class C k,ℓ are characterized by their local epigraph representation. Note that the equivalence is not true for less regular domains, though. We can now state a definition for principal curvatures and mean curvature for C 2 domains. Definition 3.4 (Curvature). Let Ω a domain with C 2 boundary and x 0 ∈ ∂Ω. Without loss of generality, we assume that the local epigraph representation of Ω around x 0 involves the tangent space to ∂Ω at x 0 , i.e. S = T x 0 and the inner normal vector a = -n(x 0 ). Using the notation of Definition 3.3, the application φ ∈ C 2 (B ∩ T x 0 , R) satisfies φ(0) = 0 and ∇φ(0) = 0.

• The principal curvatures κ 1 , . . . , κ d-1 of ∂Ω at x 0 are the eigenvalues of the Hessian matrix ∇ 2 φ(0).

• The mean curvature of ∂Ω at x 0 is defined by

H(x 0 ) = d-1 i=1 κ i = ∆φ(0).
Remark 3.2. The curvatures are associated to the variations of the outer normal vector n(x) as x moves along the boundary ∂Ω. Namely, n and H satisfy

H = div ∂Ω (n), (3.1)
where div ∂Ω is the tangent divergence operator [START_REF] Henrot | Variation et optimisation de formes[END_REF], Definition 5.4.7).

Remark 3.3. The definition of boundary curvature is still valid for domains of class C 1,1 . In such domains, the normal vector n is a Lipschitz continuous function defined on ∂Ω. Due to Rademacher's theorem, the normal vector function has derivatives at almost every point in ∂Ω, and curvature can be defined at these points using (3.1).

Domains with Lipschitz boundary

Lipschitz regularity for sets is a weaker condition than C k regularity for k ⩾ 1. It is based on local epigraphs involving Lipschitz continuous functions instead of C k functions. We state a definition of Lipschitz sets in the case of domains with compact boundary (for instance bounded domains). When ∂Ω is not compact, several variants of the Lipschitz regularity exist, but they become all equivalent when ∂Ω is compact [START_REF] Delfour | Advances in Design and Control. Metrics, analysis, differential calculus, and optimization[END_REF], Chapter 2, Definition 5.1 and Theorem 5.1), which will always be the case in this work.

Definition 3.5 (Lipschitz domain). A domain Ω with compact boundary is said to have a Lipschitz boundary when there exists a constant c > 0 such that Ω is locally the epigraph of a Lipschitz continuous function with Lipschitz constant c.

In other words, a Lipschitz domain is locally a C 0,1 -epigraph involving functions φ that share a uniform Lipschitz constant. This condition is stronger than being a domain of class C 0,1 : A Lipschitz domain is of class C 0,1 but some domains of class C 0,1 are not Lipschitz. Grisvard (2011, Chapter 1, Section 1.2) gives more details about the differences between those definitions.

Lipschitz domains are omnipresent when it comes to solving real-world problems. For instance, mechanical or thermal simulations involving industrial pieces with sharp edges or corners require to solve partial differential equations on Lipschitz domains. In particular, the normal n is defined almost everywhere on a Lipschitz boundary, as a consequence of Rademacher's theorem, and Green's formula is valid on Lipschitz domains (Grisvard, 2011, Chapter 1, Theorem 1.5.3.1), which is a good start for solving partial differential equations.

Shape functional and shape derivative

In this section, we introduce the notion of shape functional. A shape functional is a real-valued function which depends on a variable shape. Though continuity for such functionals can be defined using a topology on shapes sets, the absence of Hilbert space structure on shape sets makes it more difficult to define differentiability with respect to a shape. We first provide a definition for shape functionals, and then we introduce Hadamard's framework for shape differentiation.

In this document, a shape functional is defined as follows (Delfour and Zolésio, 2011, Chapter 7, Definition 2.1). Definition 3.6 (Shape functional). Let D a nonempty subset of R d and P(D) = {Ω | Ω ⊂ D}. A shape functional is an application J : A → R, where A ⊂ P(D), such that, for Ω ∈ A, and T a homeomorphism of D such that T (Ω) = Ω, the equality J(Ω) = J(T (Ω)) holds.

Common choices for D include D = R d or D a compact subset of R d . The above definition emphasizes the independence of J with respect to continuous transformation that preserve Ω. In the case of open Lipschitz sets, it is common to apply a small deformation field onto the shape boundary to change the value of J.

Boundary variation framework

To study the variations of a shape functional J, it is convenient to use Hadamard's boundary variation method (see for instance [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF][START_REF] Henrot | Variation et optimisation de formes[END_REF], and references therein). To represent small variations of a bounded Lipschitz subset Ω of R d , we use the notation Ω V to denote the set Ω transported by a displacement field

V ∈ W 1,∞ (R d , R d ), namely Ω V = (Id +V )(Ω) = {x + V (x), x ∈ Ω} . (3.2)
The definition for shape differentiability we give is stated in Allaire (2007, Chapter 6, Definition 6.15).

Definition 3.7 (Shape differentiability).

Let Ω a bounded subset of R d with Lipschitz boundary.

The shape functional J is said to be shape differentiable in the direction

V ∈ W 1,∞ (R d , R d ) whenever the application t → J(Ω tV ) is differentiable at t = 0. If it exists, the directional shape derivative of J is defined by dJ(Ω)(V ) = lim t↘0 J(Ω tV ) -J(Ω) t . (3.3)
In addition, J is said to be Fréchet differentiable at Ω whenever the application

V → J (Ω V ) is Fréchet differentiable at 0 on W 1,∞ (R d , R d ).
This definition of shape differentiability, which involves Lipschitz continuous displacement fields, is appropriate for numerical computations, but several versions exist in the literature. Delfour and Zolésio (2011, Chapter 8, Definition 3.3) use a more general definition. They only consider C ∞ displacement fields with compact support and define the shape gradient as a distribution on R d . Sokolowski and Zolésio (1992, Chapter 2, Definition 2.19) state the definition with a continuous displacement field. Though, all these definitions consist in applying a displacement field to Ω.

We now state standard examples of shape differentiability results for shape functionals defined as integrals on a variable domain [START_REF] Henrot | Variation et optimisation de formes[END_REF], Theorems 5.2.2 and 5.4.17). In the second proposition, we use the notation

C 1,∞ = W 1,∞ ∩ C 1 . Proposition 3.2. Let Ω ⊂ R d a shape with Lipschitz boundary and f ∈ W 1,1 (R d , R d ).
The functional defined by

J : W 1,∞ (R d , R d ) → R, V → ´ΩV f dx is differentiable at 0 and its differential reads ∀V ∈ W 1,∞ (R d , R d ) ⟨dJ(0), V ⟩ = ˆ∂Ω (V • n)f ds. (3.4) Proposition 3.3. Let Ω ⊂ R d a shape with C 2 boundary and g ∈ W 2,1 (R d , R d ). The functional defined by J : C 1,∞ (R d , R d ) → R, V → ´∂Ω V g ds is differentiable at 0 and its differential reads ∀V ∈ C 1,∞ (R d , R d ) ⟨dJ(0), V ⟩ = ˆ∂Ω 0 (V • n) ∂g ∂n + Hg ds, (3.5)
where H is the mean curvature of ∂Ω.

A wide range of similar results exist in the literature. Standard shape functionals include functions based on the solution of a partial differential equation defined on Ω, or eigenvalues of differential operators such as the Laplacian. The interested reader is refered to Henrot and Pierre (2005, Chapter 5) and Delfour and Zolésio (2011, Chapter 9) for more information.

Structure theorem

Note that, in Proposition 3.2 and Proposition 3.3, the differential dJ(0) has its support included in ∂Ω. Moreover the shape derivatives depend only on the normal component of V on ∂Ω. Indeed, the variation of J only depends on the movements of the shape boundary ∂Ω. In addition, if Ω has a C 1 boundary and V only has a tangential component on ∂Ω, then Ω is not changed by V at first order. This property of shape derivatives is stated in the structure theorem below (Henrot and Pierre, 2005, Theorem 5.9.2).

Theorem 3.1 (Structure). Let k ⩽ 1, Ω ⊂ R d an open set with C k+1 boundary and J a shape functional. Assume that the application J : C k,∞ (R d , R d ) → R, V → J(Ω V ) is Fréchet differentiable at 0. Then there exists a linear form ℓ, continuous on C k (∂Ω), such that ∀V ∈ C k,∞ (R d , R d ) d J(0), V = ℓ(V | ∂Ω • n),
where n is the normal to ∂Ω.

Orthogonal projection and signed distance function

Given a domain Ω with Lipschitz boundary and a point y ∈ R d , one can wonder whether y is inside or outside Ω, or how far from Ω y is. The signed distance function and the orthogonal projection are two tools to answer such questions. For a detailed and exhaustive presentation, the reader may refer to Delfour and Zolésio (2011, Chapter 7). Unless specified otherwise, ∥ • ∥ and • denote the Euclidean norm and scalar product on R d , respectively.

As we are working with the signed distance to ∂Ω, we need a tool to differentiate functions based on a minimum. Such a theorem is provided by Danskin (1967, Chapter 3, Theorems I and II) in the case of functions defined on spaces of finite dimension.

Lemma 3.1 (Danskin's theorem). Let f : E × K → R a continuous function, with E a vector space of finite dimension and K a compact subset of a Euclidean space. Assume that (y, z) → ∂f ∂y (y, z) exists and is continuous, and define φ :

E → R, y → min z∈K f (y, z).
If y is a point of E and w ̸ = 0 is a direction, then φ has a directional derivative in the direction w at the point y, which reads

lim t↘0 φ(y + tw) -φ(y) t = min z∈K(y) ∂f ∂y (y, z), w ,
where

K(y) = {z ∈ k | f (y, z) = φ(y)} is the set of points in K which yield the minimum of f (y, •).
In addition, assume that K(y) is a singleton. Then φ is differentiable in the sense of Fréchet at y. Last but not least, if K(y) is a singleton for y in a neighborhood, then φ is C 1 on this neighborhood.

Generalities

Here we give a few definitions and properties about the signed distance function and the orthogonal projection.

∀y ∈ R d , d Ω (y) =      d(y, ∂Ω) if x / ∈ Ω 0 if x ∈ ∂Ω -d(y, ∂Ω) if x ∈ Ω,
where d(y, ∂Ω) = min x∈∂Ω ∥y -x∥ is the Euclidean distance between y and ∂Ω. Definition 3.9 (Orthogonal projection). Let Ω be an open domain. The set of orthogonal projections of y onto ∂Ω is defined by

Π ∂Ω (y) = {x ∈ ∂Ω | ∥x -y∥ = d(y, ∂Ω)} .
When Π ∂Ω (y) is a singleton, we denote by p ∂Ω (y) its element. The set of points which have several orthogonal projections onto ∂Ω is called the skeleton of ∂Ω and is denoted by Σ.

We state the following properties (Delfour and Zolésio, 2011, Chapter 7, Theorems 2.1 and 3.1) to better understand the relationship between the orthogonal projection and the signed distance function. To prove the result with the signed distance function, we consider two cases, depending on whether x and y are on the same side of ∂Ω. In the first case x, y ∈ Ω or x, y / ∈ Ω, and

|d Ω (x) -d Ω (y)| = |d(x, ∂Ω) -d(y, ∂Ω)| ⩽ ∥x -y∥.
In the other case x ∈ Ω and y /

∈ Ω, then [x, y] ∩ ∂Ω is nonempty. Considering z ∈ [x, y] ∩ ∂Ω, we obtain ∥y -x∥ = ∥y -z∥ + ∥z -x∥ ⩾ d(y, ∂Ω) + d(x, ∂Ω) = |d Ω (y) -d Ω (x)|
and d Ω is Lipschitz continuous on R d with Lipschitz constant at most 1. The differentiability is a consequence of Rademacher's theorem (Evans and Gariepy, 2015, Theorem 3.2).

2. Note that d 2 Ω (y) = min z∈∂Ω ∥z -y∥ 2 , with (z, y) → ∥z -y∥ 2 a C 1 function, and ∂Ω a compact set. We obtain (3.6) by applying Lemma 3.1.

3. Assume that y has a unique projection point p ∂Ω (y). Then Fréchet differentiability is also a consequence of Lemma 3.1, and we obtain ∇d 2 Ω (y) = 2(y -p ∂Ω (y)). Now, choose a point y ∈ R d which has several projection points, and let z 1 ̸ = z 2 ∈ Π ∂Ω (y). Note that the elements of Π ∂Ω (y) are equidistant from y. Equation (3.6) and the Cauchy-Schwarz inequality yield

dd 2 Ω (y)(z 1 -y) = min z∈Π ∂Ω (y) 2(y -z) • (z 1 -y) = -2∥y -z 1 ∥ 2 .
We compute the derivative in the opposite direction. As z 2 ̸ = z 1 and ∥z 1 -y∥ = ∥z 2 -y∥, the Cauchy-Schwarz inequality yields (y -z 2 ) • (y -z 1 ) < ∥z 1 -y∥ 2 , and the derivative reads

dd 2 Ω (y)(y -z 1 ) = min z∈Π ∂Ω (y) 2(y -z) • (y -z 1 ) ⩽ 2(y -z 2 ) • (y -z 1 ) < 2∥y -z 1 ∥ 2 .
We conclude that dd 2 Ω (y) is not a linear application and d 2 Ω is therefore not Fréchet differentiable at y. 4. This result is a consequence of points 1 and 3.

If ∂Ω is C 1 , each point y ∈ R d \ Σ can be associated to a unique pair of parameters consisting of its orthogonal projection onto ∂Ω and its signed distance to Ω. In this case, the normal vector is defined at every point of ∂Ω, and thus y is given by the formula

y = p ∂Ω (y) + d Ω (y) n (p ∂Ω (y)) ,
where n(x) is the outward normal vector at x ∈ ∂Ω. If ∂Ω is only Lipschitz, some pairs of the form (x, ℓ) ∈ ∂Ω × R might be associated with several points, especially if ∂Ω has an angular point at x. On Figure 3.2, y 2 and y ′ 2 are associated to the same signed distance value and projection.

The points of Σ should be handled with care in this parametrization, as they can be associated with several sets of parameters. However, when computing derivatives, it is possible to choose one orthogonal projection depending on the direction with respect to which the derivative is computed.

Ω Σ y 1 x 1 x ′ 1 y 2 y ′ 2 x 2 y 3 x 3 Figure 3
.2: Set (in gray) with its skeleton (dashed line) and the signed distance level curves (gray lines). At the bottom, y 1 ∈ Σ, has several projections. On the right, note that p ∂Ω (y 2 ) = p ∂Ω (y ′ 2 ) and d Ω (y 2 ) = d Ω (y ′ 2 ). The gradient of d Ω is defined at y 2 and y ′ 2 but not at x 2 . On the left, ∇ 2 d Ω is not defined at y 3 , as the curvature of ∂Ω is not defined at x 3 .

To handle shape functionals based on the signed distance function, we need to know how the orthogonal projection and the signed distance evolve when y moves or when Ω is subject to a deformation. The following propositions gather some information found in [START_REF] Henrot | Variation et optimisation de formes[END_REF], [START_REF] Delfour | Advances in Design and Control. Metrics, analysis, differential calculus, and optimization[END_REF] and [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF].

Pointwise derivatives

The regularity of d Ω away from Σ depends on the regularity of ∂Ω, and its derivatives are defined using the normal vector to ∂Ω and its curvature. The following proposition [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF] The signed distance is twice differentiable at y and

∇ 2 d Ω (y) = d-1 k=1 κ k 1 + κ k d Ω (y) e k ⊗ e k . (3.7)
Some examples of these conditions are illustrated on Figure 3.2. The signed distance function is not differentiable at y 1 , as y 1 ∈ Σ. Moreover, the signed distance Hessian at y 3 is not defined, because the curvature of ∂Ω is not defined at x 3 . Remark 3.4. If d Ω is differentiable at y but the normal is not defined at x, it is still possible to write ∇d Ω (y) = y-x d Ω (y) (see y 2 and y ′ 2 in Figure 3.2).

The following proposition (Delfour and Zolésio, 2011, Ch. 7, Th. 4.3) provides more specific information about the regularity of d Ω in a vicinity of ∂Ω. Proposition 3.6. Let Ω be a bounded domain.

1. If Ω is of class C 1,1 , then d Ω has C 1,1 regularity in a neighborhood W of ∂Ω. Moreover, ∀y ∈ W ∇d Ω (y) = n(p ∂Ω (y)) 2. If Ω is of class C k , with k ⩾ 2, then d Ω is C 2 in a neighborhood W of ∂Ω and d 2 Ω is C 1,1 in W .
An expression for the Hessian of d Ω is given by (3.7).

Working with a C 1,1 set, i.e. a set with bounded curvature, provides additional comfort as Σ is kept away from ∂Ω. For instance the set of equation y > x 4/3 has C 1 boundary but not bounded curvature (See Figure 3.3). Thus, d Ω is not differentiable on a neighborhood of ∂Ω.

We finally state a proposition (Dapogny, 2013, Lemma 4.2) about the pointwise derivative of the orthogonal projection operator at points far from the skeleton of Ω.

Proposition 3.7. Assume Ω has C 2 boundary, y ∈ R d \ Σ and x = p ∂Ω (y). Denote by (e 1 , • • • , e d-1
) the directions of principal curvatures at x and κ 1 , • • • , κ d-1 the associated principal curvatures. The orthogonal projection is Fréchet differentiable at y and its Jacobian matrix reads

∇p ∂Ω (y) = d-1 k=1 1 - κ k d Ω (y) 1 + κ k d Ω (y) e k ⊗ e k . (3.8)
Proof. According to Proposition 3.5, the orthogonal projection of y is unique and given by the formula p ∂Ω (y) = y -d Ω (y)∇d Ω (y).

As the hypotheses ensure the existence of the signed distance Hessian at y, the expression above is differentiable, with

∇p ∂Ω (y) = I -∇d Ω (y) ⊗ ∇d Ω (y) -d Ω (y) ∇ 2 d Ω (y).
To conclude, we use (3.7) and we notice that

I -∇d Ω (y) ⊗ ∇d Ω (y) = I -n(x) ⊗ n(x) = d-1 k=1 e k ⊗ e k .

Shape derivatives

In this section, we discuss shape differentiability for the signed distance function to a given set. To be more specific, we only consider the squared signed distance d 2 Ω , which is actually the squared distance from ∂Ω. We will see that there are some similarities between expressions for derivatives with respect to the shape and derivatives with respect to the point. Indeed, sometimes, deciding whether it is the point or the set that moves is only a question of reference frame.

In the following propositions (Dapogny, 2013, Section 4.2), we state Lipschitz continuity and directional differentiability with respect to the shape for d 2 Ω .

Proposition 3.8. Let y ∈ R d , Ω a bounded Lipschitz domain, and (Allaire, 2007, Chapter 6, Lemma 6.13). Assume d(y, ∂Ω w 2 ) ⩾ d(y, ∂Ω w 1 ) and choose

B W 1,∞ (0, 1) = w ∈ W 1,∞ (Ω) ∥w∥ W 1,∞ < 1 , the application B W 1,∞ (0, 1) → R + , w → 1 2 d 2 (y, ∂Ω w ) is Lipschitz continuous. Proof. Let w 1 , w 2 ∈ B W 1,∞ (0, 1). Because ∥w 1 ∥ W 1,∞ , ∥w 2 ∥ W 1,∞ < 1, (Id +w 1 ) and (Id +w 2 ) are W 1,∞ diffeomorphisms
x 1 ∈ Π ∂Ω w 1 (y). Then for x 2 ∈ ∂Ω w 2 d(y, ∂Ω w 2 ) -d(y, ∂Ω w 1 ) ⩽ ∥y -x 2 ∥ -∥y -x 1 ∥ ⩽ ∥x 2 -x 1 ∥.
The previous inequality holds for any x 2 ∈ ∂Ω w 2 . In particular, we may consider x = (Id +w 1 ) -1 (x 1 ) and set x 2 = x + w 2 (x). With these notations,

∥x 2 -x 1 ∥ = ∥x + w 2 (x) -x -w 1 (x)∥ = ∥w 2 (x) -w 1 (x)∥ ⩽ ∥w 2 -w 1 ∥ W 1,∞ .
Note also that w → d(y, ∂Ω w ) is bounded on B w 1,∞ (0, 1), with

d(y, ∂Ω w ) = min x∈∂Ω ∥x + w(x) -y∥ ⩽ min x∈∂Ω ∥x -y∥ + ∥w(x)∥ ⩽ d(y, ∂Ω) + 1.
We conclude the proof by noting that

1 2 d 2 (y, ∂Ω w 2 ) -1 2 d 2 (y, ∂Ω w 1 ) = 1 2 d(y, ∂Ω w 2 ) + d(y, ∂Ω w 1 ) d(y, ∂Ω w 2 ) -d(y, ∂Ω w 1 ) ⩽ d(y, ∂Ω) + 1 ∥w 2 -w 1 ∥ W 1,∞ . Proposition 3.9. Let y ∈ R d , and Ω a bounded Lipschitz domain. The application w → 1 2 d 2 (y, Ω w ) defined on W 1,∞ (Ω) has directional derivatives at 0 with lim t↘0 1 2 d 2 (y, Ω tw ) -1 2 d 2 (y, Ω) t = min z∈Π ∂Ω (y) w(z) • (z -y).
(3.9)

Proof. We use the notations

g y (t) = 1 2 d 2 (y, ∂Ω tw ) = min x+tw(x)∈∂Ω tw 1 2 ∥x + tw(x) -y∥ 2 ,
and, for t small enough for tw to be a diffeomorphism,

K(t) = x ∈ ∂Ω | 1 2 d 2 (y, ∂Ω tw ) = 1 2 ∥x + tw(x) -y∥ 2 = (I + tw) -1 (Π ∂Ω tw (y)). (3.10)
For t close to 0, K(t) is a nonempty subset of ∂Ω. Also, the application (t, x) → 1 2 ∥x + tw(x) -y∥ 2 has a continuous partial derivative with respect to t. The result is a consequence of Lemma 3.1.

Continuity and mixed boundary conditions

As mentioned earlier, in this document we are working with functionals based on geometric operations such as the orthogonal projection, and we consider shapes which are deformed using a displacement field. For this reason, it is important to give a meaning to the value of the displacement field at a single point. In other words, we need to manipulate displacement fields that are continuous and obtain convergences for the uniform convergence norm.

In this section, we enumerate a few continuity results for elliptic equations and for elasticity equations with mixed boundary conditions.

Some notation

From now on, we consider that Ω is an open domain with Lipschitz boundary ∂Ω. We denote by ∂Ω N the open part of the boundary where a Neumann condition applies, and ∂Ω D = ∂Ω \ ∂Ω N . The space C(Ω) contains functions that are continuous on Ω, endowed with the uniform convergence norm (the L ∞ norm). For κ > 0, we also define the space C 0,κ (Ω) of κ-Hölder continuous functions, endowed with the norm

∥u∥ C 0,κ (Ω) = ∥u∥ L ∞ (Ω) + sup x̸ =y |u(x) -u(y)| |x -y| κ ,
and for 1 ⩽ p ⩽ ∞, W 1,p (Ω) is the Sobolev space of L p functions whose first-order derivatives are also in L p (Ω), where things are measured with the norm

∥u∥ W 1,p (Ω) = ˆΩ |u| p + |∇u| p dx 1/p .
Those two spaces are interesting from the point of view of continuous solutions. To be more specific, for p > d, W 1,p (Ω) embeds continuously into C 0,κ (Ω) where κ = 1-d/p [START_REF] Brezis | Analyse fonctionnelle. Collection Mathématiques Appliquées pour la Maîtrise[END_REF], Theorem IX.12), and for κ > 0, the injection C 0,κ (Ω) → C(Ω) is compact (Brezis, 1983, Theorem IV.24), which means that a bounded sequence of C 0,κ (Ω) or W 1,p (Ω) (with p > d) converges uniformly, up to a subsequence. For this reason, we would like to see our solutions in such a space.

However, two main obstacles appear when it comes to obtaining sufficient regularity for solutions of elasticity systems with mixed boundary. First, the linear elastic system is not an elliptic equation, as its energy functional does not penalize rigid-body displacements. In addition, the presence of mixed boundary conditions degrades the regularity of solutions, even in C 1 domains. [START_REF] Gröger | A W 1,p -estimate for solutions to mixed boundary value problems for second order elliptic differential equations[END_REF] addressed the second difficulty by proposing a class of domains adapted to the study of partial differential equations with mixed boundary conditions. We now give a word about his framework.

Gröger regular sets for mixed-boundary conditions

To work with mixed boundary conditions, it is interesting to consider the subset G = Ω ∪ ∂Ω N of R d . The subset G contains all necessary information about the regularity of Ω and the distribution of Dirichlet and Neumann conditions on its boundary. Like in the definition of C k domains, the definition of Gröger regular sets involves diffeomorphisms from a part of the surface to the unit ball of R d . Before we give a formal definition, we define the following subsets of R d : Figure 3.4 shows the sets B + , B + ∪ D and B + ∪ D 0 . These three cases correspond to the case x ∈ ∂Ω \ ∂Ω N , x ∈ ∂Ω N and x ∈ ∂(∂Ω N ), respectively. A simpler way to express this definition is that a bounded subset G is Gröger-regular when ∂Ω N has a W 1,∞ boundary in ∂Ω [START_REF] Gröger | A W 1,p -estimate for solutions to mixed boundary value problems for second order elliptic differential equations[END_REF]. A very similar definition for Gröger-regular sets has been given independently by Droniou (2000, Condition (2.11)), where the author describes ∂Ω N and ∂Ω D as "well separated" when the regularity assumption is respected.

B = {x ∈ R d | ∥x∥ < 1} B + = {x ∈ B | x d > 0} D = {x ∈ B | x d = 0} D 0 = {x ∈ D | x 1 < 0}
Haller-Dintelmann et al. (2009, Section 5) give a few more details about Grögerregular sets in dimension 2 and 3. In particular, they note that in dimension 2, a domain G is Gröger-regular when ∂Ω is W 1,∞ , ∂Ω N is a finite union of open arcs in ∂Ω and ∂Ω D is a finite union of closed arcs, none of which is reduced to a single point. They also state the following characterization in dimension 3. Proposition 3.10. Assume Ω is a bounded domain of class C 0,1 and ∂Ω N is open in ∂Ω. Then G = Ω ∪ ∂Ω N is Gröger-regular if and only if the two conditions below hold:

1. ∂Ω D is the closure of its interior.

For

x ∈ ∂(∂Ω N ), there is a neighborhood U of x and a Lipschitz diffeomorphism Ψ such that Ψ(U ∩ ∂(∂Ω N )) = (-1, 1).
This characterization, which is not true in dimension larger than 3, implies that W 1,∞ polyhedrons of R 3 where ∂(∂Ω N ) is a finite union of line segments are Gröger regular sets.

Linear elliptic equations

Though the linear elasticity system does not fall into this category, we first present some regularity results in the case of linear elliptic equations, for which a wider literature exists. Here we only state results for scalar equations, keeping in mind that uncoupled systems of dimension d reduce to d scalar equations.

A convenient space to study problems with mixed boundary conditions is the Sobolev space W 1,p D (Ω), defined as the closure of the set

u ∈ C ∞ c (R d ) supp(u) ∩ ∂Ω D = ∅ (3.11)
for the usual norm in W 1,p (Ω). Functions contained in W 1,p D (Ω) have a zero trace on ∂Ω D but can take some nonzero values on ∂Ω N .

The topological dual of

W 1,p ′ D (Ω) is denoted by W -1,p D (Ω)
, where p ′ is the conjugate exponent to p, defined by 1/p + 1/p ′ = 1. This dual contains, among others, linear forms of type v → ´Ω f • v dx for a measurable function f , but also distributions supported by a submanifold of dimension d -1 included in G. In particular, it contains applications of the form v → ´∂Ω N g • v dx, which correspond to Neumann conditions. From now on, we consider that ∂Ω D has a nonzero relative Lebesgue measure in ∂Ω. We consider (scalar) linear partial differential equations of type 3.12) has a unique solution u ∈ H 1 D (Ω). In the context of Gröger-regular sets, [START_REF] Gröger | A W 1,p -estimate for solutions to mixed boundary value problems for second order elliptic differential equations[END_REF] studies (3.12) in terms of the operator A :

     -div(A∇u) + bu = f in Ω A∇u • n = g on ∂Ω N , u = 0 on ∂Ω D (3.12) where A ∈ L ∞ (Ω, M d (R)) satisfies the ellipticity condition ∀x ∈ Ω, ∀z ∈ R d m∥z∥ 2 ⩽ z • A(x)z ⩽ M ∥z∥ 2 (3.13) for some M > m > 0 and b ∈ L ∞ (Ω, R) is a nonnegative function. It is known, due to Lax-Milgram's theorem, that for f ∈ H -1 (Ω) and g ∈ H -1/2 (∂Ω N ), (
W 1,p D (Ω) → W -1,p D (Ω) defined by ⟨Au, v ⟩ = ˆΩ (A∇u • ∇v + buv) dx. (3.14)
His framework also includes the study of non-linear elliptic equations and the evolution of a regularity result when the domain is transformed using a Lipschitz diffeomorphism.

We state here a simplified corollary of his regularity result (Gröger, 1989, Theorem 1).

Proposition 3.11. Assume that Ω∪∂Ω N is Gröger-regular. If b > m almost everywhere in Ω, then there is a q > 2 such that for all p ∈ [2, q], the mapping A is a continuous isomorphism from W 1,p D (Ω) to W -1,p D (Ω).
In other words, if for such a p the application v → ´Ω f v dx + ´∂Ω N gv ds is in

W -1,p D (Ω), then the solution of (3.12) is in W 1,p D (Ω).
The author also provides for A -1 an estimate which depends on the ellipticity constant (1 -m/M ) 1/2 and the domain regularity (see also [START_REF] Gröger | Resolvent estimates in W -1,p for second order elliptic differential operators in case of mixed boundary conditions[END_REF]. Though this article introduces a framework for the study of mixed boundary conditions, the result in Proposition 3.11 is not fully satisfying. Indeed, in our application case, Ω represents an organ, and we cannot easily make sure that it is regular enough so that q > d.

Remark 3.5. Even in the case of Poisson's equation (A = I), the presence of mixed boundary conditions degrades the regularity of solutions. While solutions to the Dirichlet problem are in W 1,p for every p > 2, no regularity better than that of Proposition 3.11 can be expected.

A Hölder regularity result is proposed by Droniou (2000, Theorem 2.1). Though the author does not cite Gröger's work, the hypotheses he sets on the domain are similar to the definition of Gröger-regular sets. Once again, the author handles a general class of linear elliptic problems, involving for instance nonhomogeneous Dirichlet boundary conditions, a transport term and several configurations of boundary conditions. The regularity result adapted to our case is stated below. Proposition 3.12. Assume Ω ∪ ∂Ω N is Gröger-regular and choose p > d.

If f ∈ W 1,p ′ (Ω) ′ and g ∈ W 1-1/p ′ ,p ′ (∂Ω N ) ′
, then there exists κ ∈ (0, 1 -d/p] depending only on Ω, p, m, M such that the solution u of (3.12) is κ-Hölder continuous in Ω.

In addition, if there is a constant Λ > 0 such that

∥b∥ L dp/(d+p) (Ω) + ∥f ∥ W 1,p ′ (Ω) ′ + ∥g∥ W 1-1/p ′ ,p ′ (∂Ω N ) ′ ⩽ Λ, then there is a constant C = C(Ω, ∂Ω N , m, M, p, Λ) such that ∥u∥ C 0,κ (Ω) ⩽ C.
In Proposition 3.12, the Neumann data g is supposed to be a linear form which applies on W 1-1/p ′ ,p ′ (∂Ω N ). This space contains the traces of functions in W 1,p ′ D (Ω) on ∂Ω N . Therefore, f and g define elements of W -1,p D (Ω) and another way to formulate the result is that A -1 maps continuously W -1,p D (Ω) into C 0,κ (Ω). Years later, using Gröger's formalism, [START_REF] Haller-Dintelmann | Hölder continuity and optimal control for nonsmooth elliptic problems[END_REF] obtain a result similar to Proposition 3.12. In addition to the characterization of Gröger-regular sets mentioned above, they also give precisions about the interest of Hölder continuity results in the context of optimal control.

Elasticity equation

We now turn to the linear elasticity system. Though we already introduced the linear elasticity system in Section 2.2, we give some additional details here concerning the regularity of solutions. Unless stated otherwise, functions in this paragraph take their value in R d .

As opposed to elliptic equations, the elastic energy only penalizes the symmetric part of the displacement gradient (also known as the linearized strain tensor), denoted ε(u) = 1 2 (∇u T + ∇u). In this section A denotes the elasticity tensor, a linear operator which satisfies an ellipticity condition with respect to the strain tensor of the form

∀x ∈ Ω ∀ε ∈ M d (R) such that ε T = ε m∥ε∥ 2 ⩽ Aε : ε ⩽ M ∥ε∥ 2
for some M > m > 0. In the context of the standard elasticity system, the linear elasticity tensor is defined by Hooke's law

Aε(u) = 2µε(u) + λ tr(ε(u)).
Remember that the linear elasticity equation reads

     div(Aε(u)) = f in Ω 0 Aε(u) • n = g on ∂Ω N , u = 0 on ∂Ω D (3.15)
with g ∈ H -1/2 (∂Ω N ) a surface force field and f ∈ H -1 (Ω) stands for volume forces.

The weak formulation of (3.15) reads as follows:

find a displacement field u ∈ H 1 (Ω 0 ) such that ∀v ∈ H 1 D (Ω 0 ) ˆΩ0 Aε(u) : ε(v) dx = ˆΩ f • v dx + ˆ∂Ω N g • v ds.
A key result to establish the well-posedness of (3.15) is the equivalence between norms based on the strain tensor with the H 1 norm. We state it here as a lemma (Ciarlet, 1988, Chapter 6, Theorems 6.3-3 and 6.3-4).

Lemma 3.2. Denote the space

K D (Ω) = u ∈ L 2 (Ω) | ε(u) ∈ L 2 (Ω) and u = 0 on ∂Ω D .
The spaces K D (Ω) and H 1 D (Ω) are actually the same space. In addition, if ∂Ω D has positive measure in ∂Ω 0 , then the norm u → ∥ε(u)∥ L 2 (Ω) is equivalent to the H 1 norm.

Due to the presence of the symmetric gradient ε(u) instead of ∇u, the linear elastic system is not included in the framework of elliptic partial differential equations, and as a consequence, fewer works exist about it in the literature. We present a few regularity works where the authors adapt techniques used in the context of elliptic systems to the elasticity system.

An early work is proposed by [START_REF] Shi | Higher integrability of the gradient in linear elasticity[END_REF], where the authors adapt techniques used for elliptic equations to the linear elasticity system in a domain with C 1 boundary. In the case of data in W -1,p D (Ω), they obtain W 1,p -estimates for p in an interval (2, q) where q depends only on the domain and the elliptic constants of the system m and M .

An interesting result is proposed by [START_REF] Herzog | Integrability of displacement and stresses in linear and nonlinear elasticity with mixed boundary conditions[END_REF]. The authors adapt Gröger's framework to linear and non-linear equations that satisfy an ellipticity condition with respect to the strain tensor ε(u). We state their theorem in the case of linear elasticity. Proposition 3.13. Assume Ω ∪ ∂Ω N is Gröger-regular. There is a q > 2 such that for all p ∈ [2, q], the system (3.15) has a unique solution u ∈ W 1,p D (Ω) which depends Lipschitz continuously on the data, provided that f and g define elements of W -1,p D (Ω). Namely, there is a constant C such that

∥u∥ W 1,p (Ω) ⩽ C (∥f ∥ W -1,p (Ω) + ∥g∥ W -1,p (Ω) .
This result is the best we found in the literature. Though it ensures Hölder continuity of solutions in dimension 2, it is not sufficient in dimension 3.

We end this section by quoting a regularity result in dimension 3 from Ciarlet (1988, Chapter 6, Theorem 6.3-6) which requires more regular data. The theorem is stated for Dirichlet boundary conditions, but the author explains that the result still holds in the case of Neumann conditions, i.e. ∂Ω N = ∂Ω. Though, the result is not true in the case of mixed boundary conditions, as the solution might be irregular in a vicinity of the interface between ∂Ω N and ∂Ω D .

Solutions to a problem with Neumann boundary conditions are defined up to an infinitesimal rigid displacement. We denote by T (Ω) = {w | ε(w) = 0} the set of rigid translation in Ω. In particular, such displacements satisfy ´Ω0 Aε(w) : ε(v) dx = 0 for any test function v, which means that they do not cost any elastic energy. Solutions to the elastic problem are measured with the seminorm

∥u∥ W k,p (Ω)/T (Ω) = inf w∈T (Ω) ∥u + w∥ W k,p (Ω) .
(3.16)

As rigid translations are not penalized by the elastic deformation energy, they should not be encouraged by the external forces. In order for solutions to exist, the forces f and g should satisfy the following compatibility conditions

∀w ∈ T (Ω) ˆΩ f • w dx + ˆ∂Ω g • w ds = 0. (3.17)
Keeping this condition in mind, we state the result for Neumann boundary conditions.

Proposition 3.14. Assume ∂Ω ⊂ R 3 is of class C 2 with ∂Ω N = ∂Ω. Let p ∈ [6/5, ∞), f ∈ L p (Ω)
and g ∈ W 1-1/p,p (∂Ω) such that (3.17) holds, then the solutions of (3.15) are in W 2,p (Ω), and there is a constant C such that

∥u∥ W 2,p (Ω)/R(Ω) ⩽ C ∥f ∥ L p (Ω) + ∥g∥ W 1-1/p,p (Ω) .
Chapter 4

Existence of solutions and optimality conditions

Properties of the functional

To better understand the properties of the optimization problem, we first give a closer look at the functional (2.9). Remember that its expression reads

J(u) = ˆΓ j y (u) dy, where j y (u) = 1 2 d 2 (y, S u ). (4.1)
Here, j y denotes an elementary functional associated to a single point y ∈ Γ and S u denotes the part of the deformed boundary (Id +u)(∂Ω 0 ) that should match Γ. If Γ is a discrete set, for instance a point cloud, the integral in (4.1) is replaced by a sum. Remember also that, a priori, J(u) is defined for a continuous displacement field u ∈ C(Ω 0 ).

As J measures the discrepancy between the deformed surface S u and the observed data Γ, its properties depend on the geometric relationship between those two surfaces. In particular, the orthogonal projection plays a key role as a transformation between Γ and S u . On the one hand, the Gateaux differentiability of J depends on whether points from Γ have a unique projection point onto S u . On the other hand, the regularity of descent directions is degraded when several points from Γ project onto the same point in S u . In this section, we first discuss the functional differentiability, and then we mention the possibility to generate regular descent directions.

Directional differentiability and more

In this section, we study the differentiability of the functional J. In the context of optimization problems, computing derivatives is useful to find descent directions, i.e. directions that make the objective function decrease, and also to assess optimality. In our case, J is not always Gateaux differentiable, as for certain displacement fields u, J cannot be approximated around u using an affine function. However, we prove below that J has directional derivatives. In addition, the expression (4.1) involves a min operation (hidden in the distance term). We show below that the min operator also appears in the derivative expression, which suggests that it is easier to make J decrease than increase. In other words, non-differentiability results in an abundance rather than in a lack of descent directions.

We now derive an expression for the derivatives of J. In a first technical lemma, we focus on the elementary functional j y defined in (4.1). From there, we use the dominated convergence theorem to obtain a result for J.

As j y is expressed in terms of distance to S u , we compute derivatives using Danskin's theorem [START_REF] Danskin | The theory of max-min and its application to weapons allocation problems[END_REF], stated as Lemma 3.1 in this manuscript, to investigate directional differentiability. In his work, the author proposes conditions to obtain Fréchet differentiability and even C 1 regularity for the considered function, but those results are stated for finite-dimensional spaces. For this reason, we only use Lemma 3.1 for directional differentiability, by considering j y along a fixed direction. This result is very close to and inspired from Proposition 3.8 and Proposition 3.9, but here the derivative is expressed with respect to the displacement field instead of the shape.

Differentiability is related with the number of orthogonal projections of a given point y ∈ Γ onto the deformed surface S u . In this chapter, we use the notation

P y (u) = {x ∈ S 0 | ∥x + u(x) -y∥ = d(y, S u )} = (Id +u) -1 Π S u (y). (4.2)
Points in P y (u) belong to the original surface S 0 , but their images by Id +u are located at a distance d(y, S u ) from y. When P y (u) is a singleton, it means that y has a single projection point onto S u , and this projection point is the image of a single point from S 0 by Id +u.

Lemma 4.1. Assume Ω 0 ⊂ R d is a bounded set with Lipschitz boundary, and consider u a continuous displacement field and y a fixed point in Γ.

1. The application j y is locally Lipschitz continuous for the uniform convergence norm. The Lipschitz constant can be chosen independent from y ∈ Γ but still depends on Γ.

2. In addition, j y has directional derivatives at u. For v ∈ C(Ω 0 ), the directional derivative in the direction v reads

dj y (u)(v) = min x∈P y (u) v(x) • (x + u(x) -y).
where P y (u) is defined in (4.2).

Proof. Recalling that S u = (Id +u)(S 0 ) and that S 0 is compact, we write

j y (u) = min x∈S 0 1 2 ∥x + u(x) -y∥ 2 . (4.3) 1. Let u 1 , u 2 ∈ C(Ω 0 )
. Because S 0 is compact and u 1 is continuous, there exists x 1 ∈ S 0 such that j y (u 1 ) = 1 2 ∥x 1 + u 1 (x 1 ) -y∥ 2 . Now, using (4.3),

j y (u 2 ) -j y (u 1 ) ⩽ 1 2 ∥x 1 + u 2 (x 1 ) -y∥ 2 -1 2 ∥x 1 + u 1 (x 1 ) -y∥ 2 = 1 2 (2x 1 + u 2 (x 1 ) + u 1 (x 1 ) -2y) • (u 2 (x 1 ) -u 1 (x 1 )) ⩽ max x∈S 0 ∥x -y∥ + 1 2 ∥u 2 + u 1 ∥ L ∞ (Ω 0 ) ∥u 2 -u 1 ∥ L ∞ (Ω 0 ) .
As Γ is compact, the previous Lispchitz constant is controlled by

L = max y∈Γ max x∈S 0 ∥x -y∥ + 1 2 ∥u 2 + u 1 ∥ L ∞ (Ω 0 ) . 2. Let v ∈ C(Ω 0 ) a direction.
To compute the derivative of j y in the direction v, we apply Lemma 3.1 to the application

f : R + × S 0 → R (t, x) → 1 2 ∥x + u(x) + tv(x)
-y∥ 2 , and we use the notation

K(t) = {x ∈ S 0 | d(y, S u+tv ) = ∥x + u(x) + tv(x) -y∥} .
Note that for t > 0, K(t) is nonempty as S 0 is compact and u + tv is continuous. Moreover, the derivative ∂ t f exists and is continuous. As a consequence of Lemma 3.1 applied to the identity j y (u + tv) = min x∈S 0 f (t, x), j y is differentiable at u in the direction v, and its directional derivative reads

dj y (u)(v) = lim t↘0 j y (u + tv) -j y (u) t = min x∈K(0) ∂ t f (0, x),
hence the result with P y (u) = K(0).

Of course, if Γ contains a finite number of points, the derivative of J can be written using a simple sum. In the case where Γ is continuous and J is defined using an integral, we conclude using the dominated convergence theorem.

Proposition 4.1. Let Ω 0 ⊂ R d a Lipschitz bounded set, and u ∈ C(Ω 0 ). The functional (2.9) is directionally differentiable at u. Its derivative in the direction v ∈ C(Ω 0 ) is given by

dJ(u)(v) = ˆΓ min x∈P y (u) [v(x) • (x + u(x) -y)] dy, ( 4
j y (u + tv) -j y (u) t ⩽ L∥v∥ L ∞ (Ω 0 ) ⩽ Lr,
where L > 0 is independent from y ∈ Γ. Using the dominated convergence theorem, we pass to the limit in the expression

J(u + tv) -J(u) t = ˆΓ j y (u + tv) -j y (u) t dy --→ t↘0 ˆΓ dj y (u)(v) dy
to obtain the derivative. Now, if we assume that Γ has zero Lebesgue measure, formula (4.5) is easily derived from (4.4), and dJ(u) is obviously linear. Remembering that the application d(•, S u ) is continuous, we obtain

|⟨dJ(u), v⟩| ⩽ |Γ| max y∈Γ d(y, S u ) ∥v∥ L ∞ (Ω 0 ) ,
which shows that dJ is continuous.

As the expression of J is based on the distance with respect to S u , dJ is not always a linear form. Non-differentiability happens when P y (u) contains several points for a part of Γ with a positive relative Lebesgue measure. We see in the expression P y (u) = (Id +u) -1 [Π S u (y)] that, for a given y ∈ Γ, this non-uniqueness is either caused by the projection operator Π S u or by the deformation (Id +u). In the first situation, y belongs to the skeleton of S u and has several projections onto S u . In the second situation, y may have a unique projection p S u (y), but this projection is the image of several points of S 0 by the application (Id +u).

Figure 4.1 shows a configuration where Γ is included in the skeleton of S 0 (here u = 0). In this case, every point in Γ has two projections onto S 0 , while Id +u = Id is a bijection. If we consider the rigid translation v ≡ V with V ∈ R 2 , we note that dJ(0)(-v) = dJ(0)(v) < 0. Thus, dJ is not a linear application. In Figure 4.2, each point y ∈ Γ has a single projection onto S u , but this projection point is the image of both x 1 , x 2 ∈ S 0 by (Id +u). Once again, each leg of the deformed set can move closer to Γ with the same effect on J. However, if only one leg moves away, from Γ, the value of J does not change as the other leg is not moving, which is again a nonlinear behavior.

Let us look more closely at the example from Figure 4.1. If we only consider rigid deformations dJ(0)(v) and dJ(0)(-v) can be both negative but cannot be both positive, as y can move closer to x 1 or x 2 but cannot move away from x 1 and x 2 at the same time. In Figure 4.2, one leg moving closer to Γ is sufficient to make J decrease, while both legs need to move away from Γ to make J increase. Thus, it seems easier to make J decrease than increase. Due to the presence of a min operator in the derivative of J, we can generalize this observation in a tiny lemma. Γ

y x 1 x 2 V Figure 4
.1: Non-differentiable case with u = 0 and S 0 = ∂Ω 0 . The point y ∈ Γ is equidistant from x 1 and x 2 , so that P y (u) = {x 1 , x 2 }. The derivative of J has the same (negative) value in the directions

V 1 Ω 0 and -V 1 Ω 0 . Γ y a x 1 x 2 Figure 4
.2: Non-differentiable case with self-intersection and, again, S 0 = ∂Ω 0 . The dashed set is Ω 0 , while the gray set is Ω u . The displacement field u transports both x 1 and x 2 to a = p S u (y).

Lemma 4.2. Let u ∈ C(Ω 0 ) and v a direction. Then dJ(u)(-v) ⩽ -dJ(u)(v).

Proof. Let (x y ) y∈Γ a family such that x y ∈ P y (u) and

dJ(u)(v) = ˆΓ v(x y ) • (x y + u(x y ) -y) dy.
The directional derivative in the direction -v reads ˆΓ min

x∈P y (u) [-v(x) • (x + u(x) -y)] dy ⩽ ˆΓ -v(x y ) • (x y + u(x y ) -y) dy,
A simple consequence of Lemma 4.2 is that a non-differentiable point cannot be a local minimizer of J. Non-differentiable minimizers are a pet peeve of continuous optimization solvers, as it is difficult to measure optimality around such points. Fortunately, the following lemma evacuates this possibility.

Lemma 4.3. Let u ∈ C(Ω 0 ) a local minimizer of J. Then J is Gateaux differentiable at u with dJ(u) = 0.
Proof. As u is a local minimizer of J, every direction v satisfies dJ(u)(v) ⩾ 0. If v is a direction such that dJ(u)(v) > 0, then, due to Lemma 4.2, we obtain dJ(u)(-v) < 0, and -v is a descent direction. As it is in contradiction with u being a local minimizer, the only possibility is dJ(u)(v) = 0.

Last but not least, we try to obtain some insight about the first-order behavior of j y and J in the non-differentiable case. In Figure 4.1, y has two projections onto S u . If we consider all directions in the space of continuous displacement fields, we understand that, for a part of those directions, x 1 moves closer to y faster (or gets away from y slower) than x 2 , and the expression for dj y (0) involves x 1 . For the remaining bunch of directions, x 2 gets closer to y faster than x 1 and the expression for dj y (0) involves x 2 . We show below that those directions are organized in two closed cones, associated with x 1 and x 2 , respectively. On each cone, the function dj y (0) behaves like a linear application. In this very example, the cones are two half-spaces, separated by the hyperplane

H = v ∈ C(Ω 0 ) v(x 1 ) • (x 1 + u(x 1 ) -y) = v(x 2 ) • (x 2 + u(x 2 ) -y) .
The directions v ∈ H belong to both cones and leave y in the skeleton of S u+tv at first-order.

This behavior is similar for J. This time, each cone is associated with a family (x y ) y∈Γ , with x y ∈ P y (u), such that for each y, x y statisfies the min in (4.4). Such a cone may contain only one direction, or represent the whole space in the differentiable case. We state a more formal result in the following proposition.

Generate descent directions using extension-regularization

In this section, we only consider the Gateaux differentiable case. We discuss the possibility to compute descent directions to make the objective function decrease in the context of a displacement-based optimization procedure. In particular, we need descent directions to be continuous, so that the successive displacement fields in the optimization process remain in C(Ω 0 ). Steepest descent directions are usually computed by finding a representant of dJ(u) in a certain inner product. In infinite-dimensional spaces, such directions are not guaranteed to be as regular as the optimization variable u.

A classical technique to compute regular descent directions in shape optimization is the velocity extension-regularization method (see [START_REF] Doǧan | Discrete gradient flows for shape optimization and applications[END_REF]. When the objective function is Gateaux differentiable, this method consists in choosing the shape gradient in the H 1 inner product (or another elliptic inner product) as a descent direction. Due to elliptic regularity, resulting directions are more regular. We follow a similar approach and we describe several cases depending on the regularity of dJ(u).

We fix a displacement field u ∈ C(Ω 0 ). Following the extension-regularization framework, we look for descent directions

w ∈ C(Ω 0 ) ∩ H 1 D (Ω 0 ) such that ∀v ∈ C(Ω 0 ) ∩ H 1 D (Ω 0 ) ˆΩ0 A∇w : ∇v dx = -⟨dJ(u), v⟩, (4.6)
where A satisfies the ellipticity condition (3.13). In practice, the inner product associated to linear elasticity, defined by

⟨u, v⟩ el = ˆΩ0 Aε(u) : ε(v) dx,
is also used (see de [START_REF] De Buhan | An optimization method for elastic shape matching[END_REF]Dapogny et al., 2018, and associated codes).

Indeed, an elastic model is more likely to produce a deformation that does not require remeshing. Though, in this section, we stick to the elliptic case, in order to exploit theoretical results presented in Section 3.3.2.

As we are interested in continuous descent directions, we now check the hypotheses of Proposition 3.12. Hölder continuity of solutions is guaranteed when dJ(u) is an element of W -1,p D (Ω), for some p > d. In other words, dJ(u) should apply on functions in W 1,p ′ D (Ω), which is not always possible. Remember that p ′ is defined so that 1/p + 1/p ′ = 1, and thus satisfies p ′ < d/(d -1).

In order to better understand the differential dJ, we describe below a favorable case and a pathological case. The existence of different cases comes from the definition of dJ as an integral on Γ. Here, the key parameter is not the number of projections of a given point y ∈ Γ onto S u , but rather the number of points from Γ that project onto a given point a ∈ S u .

A favorable case

Though this case is based on a lot of hypotheses, it is quite common in practice. We assume that (Id +u) is a Lipschitz diffeomorphism from S 0 to S u . For instance, when Γ Π S 0 (Γ) Γ Π S 0 (Γ) Figure 4.3: Examples of favorable cases for u = 0. On the left, Γ is in bijection with Π S 0 (Γ). On the right, each connected part of Γ is in bijection with a subset of Π S 0 (Γ).

∥u∥ W 1,∞ < 1, (Id +u) is guaranteed to be a Lipschitz diffeomorphism [START_REF] Henrot | Variation et optimisation de formes[END_REF]. We also assume that the orthogonal projection operator p S u defines a Lipschitz diffeomorphism between Γ and S u . Finally, we define the Lipschitz diffeomorphism Ψ = (Id +u) -1 • p S u . In other words, there is a Lipschitz change of variable between Γ and Ψ(Γ) ⊂ S 0 . By doing the substitution in the expression of dJ(u), we obtain

⟨dJ(u), v⟩ = ˆΨ(Γ) v(x) • (x + u(x) -Ψ -1 (x))| det ∇Ψ -1 (x)| dx.
As Ψ is a Lipschitz diffeomorphism, det ∇Ψ -1 (x) is bounded, and we also note that

∥x + u(x) -Ψ -1 (x)∥ ⩽ max x∈S u ,y∈Γ ∥x -y∥. Therefore, if v ∈ W 1,1 D (Ω 0 ) and M = ∥ Id +u -Ψ -1 ∥ L ∞ ∥ det ∇Ψ -1 ∥ L ∞ , there is a con- stant C > 0 such that |⟨dJ(u), v⟩| ⩽ M ∥v∥ L 1 (∂Ω 0 ) ⩽ M C∥v∥ W 1,1 (Ω 0 ) .
Here we used the trace theorem (Evans, 2010, Section 5.5, Theorem 1) for the second inequality. Therefore, dJ(u) is represented by an element of W 1,1 D (Ω 0 ) ′ . As a consequence, it is possible to apply the regularity result and obtain Hölder continuity for descent directions.

Note that this calculation is still valid when Ψ is a W 1,p diffeomorphism with p > d. In this case, dJ(u) defines an element of L p (∂Ω 0 ). It is also possible to split Γ into a finite partition Γ 1 , Γ 2 , . . . , each element of which is in bijection with its image by Ψ, and to express dJ(u) as a sum of elements of L p (∂Ω 0 ). Figure 4.3 shows two examples of favorable cases with u = 0. In both cases, each connected part of Γ is diffeomorphic to its projection onto S u .

A pathological case

The regularity of dJ(u) is degraded when a subset of positive measure in Γ is associated by Ψ with a subset of dimension d -2 or less in ∂Ω 0 . Thus, the linear form dJ(u) can Γ Π S 0 (Γ) Γ Π S 0 (Γ) Figure 4.4: Examples of pathological cases for u = 0. In both cases, there is an accumulation point in Π S 0 , wich results in a loss of regularity for dJ(0). only apply to displacement fields for which the trace on a (d -2)-dimensional set has a meaning, for instance continuous displacement fields. An example is given in Figure 4.4. Note that this situation is likely to happen when the domain is a polyhedron, as corners tend to concentrate orthogonal projections.

An example of pathological case arises when Γ is a point cloud {y 1 , y 2 , . . . } and (Id +u) is a Lipschitz diffeomorphism. If we denote x j = Ψ(y j ), the expression of dJ(u) reads

⟨dJ(u), v⟩ = j v(x j ) • (x j + u(x j ) -y j ).
The linear form dJ applies to displacement fields v that can be evaluated at every point x j ∈ ∂Ω 0 . In other words, dJ only applies on continuous displacement fields. As displacement fields in W 1,p ′ D (Ω 0 ) are not guaranteed to be continuous, we conclude that dJ is not in W -1,p D (Ω 0 ). Therefore, the Hölder continuity result cannot be applied to a descent direction computed with the elliptic inner product.

In such a pathological case, descent directions that are generated using the extensionregularization framework are not in C(Ω 0 ), and thus cannot be used in the optimization process, at least from a continuous point of view. When solving the problem numerically in fixed dimension, regularity of descent directions is not a concern. However, if the mesh is refined, the discretized problem may approximate better the continuous problem and some irregularities may appear in the successive displacement fields generated by the optimization solver.

Existence of solutions

We now turn to the existence of solutions. Proving that solutions exist may be useful, even in a numerical context, to ensure that displacement fields returned by the registration algorithm are consistent with the continuous problem.

We use very standard tools to prove the existence of a solution, starting with a minimizing sequence in the space of controls, and then proving that the associated state sequence converges to a state that minimizes the functional J. In this section, the regularity of solutions to the elasticity system is critical to ensure the required convergences. Unfortunately, we did not find the relevant results in the literature, and we were not able to prove the existence of solutions in a general case.

We begin the section by considering a toy problem, obtained by replacing the elasticity system with an elliptic system, and we prove existence of solutions in a general case. Then, we turn to the problem involving the linear elasticity system, with a partial result. Finally, we consider a problem with stronger hypotheses, in order to obtain an existence result with the elasticity system. Remember that the feasible set G M is defined in (2.12) by

G M = g ∈ L ∞ (∂Ω N ) | ∥g∥ L ∞ (∂Ω N ) ⩽ M .
In this section, we also refer to Gröger-regular sets. This notion is defined and discussed in Section 3.3.

A full result with an elliptic equation

Before we consider the existence of solutions for the elastic model, we first consider a simpler toy problem. The main difficulty from this toy problem to the real problem (2.13) is to handle the regularity of solutions for the linear elasticity system. The toy problem is based on Poisson's equation and reads

min g∈G M J(u g ) s.t.      -∆u g + u g = 0 in Ω u g = 0 on ∂Ω D , ∂ n u g = g on ∂Ω N , (4.7)
The key result we use to prove the existence of solution is Proposition 3.12 which provides Hölder estimates for solutions of elliptic partial differential equations.

Theorem 4.1. Assume that Ω 0 is Gröger-regular. Then the toy problem (4.7) has at least one solution.

Proof. Consider a minimizing sequence (g j ) of elements of G M . As (g j ) is bounded in L ∞ (∂Ω N ), it converges ⋆-weakly in L ∞ (∂Ω N ), up to a subsequence. We denote its limit by g. Note that G M is closed for the ⋆-weak topology in L ∞ (∂Ω N ), and for this reason g ∈ G M . The proof falls into two parts. First, we prove that the sequence of associated states converges in some sense to the state associated to g . Then, we prove that the states actually converge uniformly to finish the proof.

Denote respectively by u j = u g j the state associated to g j and by u = u g the state associated to g. We first show that the sequence of states (u j ) converges toward u in some sense. The state equations for u and u j read

∀v ∈ H 1 D (Ω) ˆΩ0 ∇u j : ∇v + u j • v dx = ˆ∂Ω N g j • v ds ∀v ∈ H 1 D (Ω) ˆΩ0 (∇u : ∇v + u • v) dx = ˆ∂Ω N g • v ds.
Because (g j ) converges weakly in L 2 (∂Ω N ), we obtain

∀v ∈ H 1 D (Ω) ˆ∂Ω N g j • v ds ----→ j→+∞ ˆ∂Ω N g • v ds.
Now, using the state equations satisfied by u and u j results in

∀v ∈ H 1 D (Ω) ˆΩ0 ∇u j : ∇v + u j • v dx ----→ j→+∞ ˆΩ0 (∇u : ∇v + u • v) dx, i.e. u j ⇀ u in H 1 (Ω 0 ). In particular, u j ⇀ u in L 2 (Ω 0 )
. Now, we show that the convergence of (u j ) is actually uniform. Uniform convergence implies weak convergence in L 2 (Ω 0 ), which ensures that the uniform limit of (u j ) is u.

Droniou (2000, Remark 2.3) notes that L (d-1)p/d (∂Ω N ) ⊂ W 1-1/p ′ ,p ′ (∂Ω N ) ′ with continuous imbedding. As a consequence, the imbedding from

L ∞ (∂Ω N ) to W 1-1/p ′ ,p ′ (∂Ω N ) ′
is also continuous and (g j ) is bounded in W 1-1/p ′ ,p ′ (∂Ω N ) ′ . Due to Proposition 3.12, there exists κ > 0 such that u j is κ-Hölder continuous and there is a constant C > 0 such that

∀j ∥u j ∥ C 0,κ (Ω 0 ) < C.
In other words, the u j are bounded and equicontinuous. Due to Ascoli's theorem (Evans, 2010, Appendix C.7), up to a subsequence, (u j ) converges uniformly in C(Ω 0 ) toward u.

As J is continuous for the L ∞ (Ω) norm, J(u) is the limit of J(u j ) and u is a solution of (4.7).

Remark 4.2. When d = 2, it is also possible to prove the existence of solutions for the (unconstrained) penalized problem min g∈L 2 (∂Ω N ) J(u g ) + α 2 ∥g∥ 2 L 2 (∂Ω N ) , with α > 0. Indeed, taking p = 3, we obtain (d -1)p/d = 3/2 < 2, and thus

L 2 (∂Ω N ) ⊂ W 1-1/p ′ ,p ′ (∂Ω N ) ′ .
As the coerciveness of the penalization term makes sure that the minimizing sequence (g j ) is bounded in L 2 , Proposition 3.12 yields a uniform bound on the Hölder norm of the (u j ) and the existence of a solution follows.

Extension to the linear elastic system

Given the theorems introduced in Section 3.3.3, we cannot prove the existence of solutions for problem (2.13) with linear elasticity as state equation without any further assumption. In order to obtain Hölder estimates for solutions to the linear elastic system, we need to obtain q > d in Proposition 3.13. This is guaranteed in dimension 2, but not in dimension 3. In particular, q depends on Lamé coefficients and on the domain regularity. Since it seems restrictive to assume that the boundary of a liver is very regular, it is difficult to be sure that the domain is regular enough.

However, we can still have a look at the necessary modifications to adapt the existence proof to linear elasticity in dimension 2. We also highlight a specificity of dimension 2: the L ∞ constraint is not necessary in the proof provided that the penalty coefficient is positive (see Remark 4.2).

Theorem 4.2. Assume that d = 2 and Ω 0 is Gröger-regular. If α > 0 or M < ∞, then (2.13) has at least one solution.

Proof. As previously, consider a minimizing sequence (g j ) of points in L 2 (∂Ω N ) and denote by (u j ) the sequence of associated states. Either due to the L ∞ constraint (M < ∞) or to the penalty term coercivity (α > 0), (g j ) is bounded in L 2 (∂Ω N ) and converges weakly in L 2 toward g ∈ L 2 (∂Ω N ), up to a subsequence. As a consequence, we obtain

∀v ∈ H 1 D (Ω 0 ) ˆΩ0 Aε(u j ) : ε(v) dx → ˆΩ0 Aε(u) : ε(v) dx,
where u = u g is the state associated to g. As the linear elasticity bilinear form is a scalar product on H 1 D (Ω 0 ) (see Lemma 3.2), the previous equation means u j ⇀ u in H 1 (Ω 0 ). We now use a compact embedding to prove that (u j ) converges uniformly toward u. Using Proposition 3.13, there is a q > 2 such that g j ∈ W -1,p D (Ω 0 ) ⇒ u j ∈ W 1,p D (Ω 0 ) for any p ∈ [2, q]. We choose p = min(4, q), so that (d -1)p/d = p/2 ⩽ 2. As mentioned in Remark 4.2, the embedding L 2 (∂Ω N ) → W 1-1/p ′ ,p ′ (∂Ω N ) ′ is continuous in dimension 2. Therefore, (g j ) is uniformly bounded as a sequence of W -1,p D (Ω 0 ) and, using the estimate in Proposition 3.13, (u j ) is uniformly bounded as a sequence of W 1,p D (Ω 0 ). The Rellich-Kondrachov theorem (Adams and Fournier, 2003, Chapter 6, Theorem 6.3) states that the imbedding W 1,p (Ω 0 ) → C(Ω 0 ) is compact, and, up to a subsequence, (u j ) converges uniformly toward u.

As J is continuous for the uniform convergence norm, J(u) = lim J(u j ). In addition, (g j ) converges weakly in L 2 (∂Ω N ) and ∥g∥ L 2 (∂Ω N ) ⩽ lim inf ∥g j ∥ L 2 (∂Ω N ) . Finally, g is a solution of problem (2.13).

In dimension 3, the existence of solutions requires M < ∞ and is subject to the condition that the q in Proposition 3.13 satisfy q > 3. We did not find a sufficient condition in the literature to enforce this inequality.

A problem with stronger hypotheses

In this section, we transform problem (2.13) to apply Proposition 3.14. To this end, strong hypotheses are necessary. First, [START_REF] Ciarlet | Studies in Mathematics and its Applications. Three-dimensional elasticity[END_REF] explains that the result is not valid if the interface ∂Ω D ∩ ∂Ω N is not empty. For this reason, we consider the problem with Neumann conditions. For the sake of simplicity, we consider the case p = 2.

The surface loading field g must be the trace on ∂Ω 0 of a function in H 1 (Ω 0 ). It is measured using the norm

∥g∥ H 1/2 (∂Ω 0 ) = inf ∥w∥ H 1 (Ω 0 ) , w ∈ H 1 (Ω 0 ), w| ∂Ω 0 = g .
The problem we consider is

min g∈W M J(u g ) subject to constraint      div(Aε(u)) = 0 in Ω 0 Aε(u) • n = g on ∂Ω 0 , ´Ω0 u dx = 0 (4.8)
where Proof. We denote by (g j ) a minimizing sequence. As (g j ) is bounded in H 1/2 (∂Ω 0 ), it is bounded in L 2 (∂Ω 0 ) and converges weakly in L 2 (∂Ω 0 ), up to a subsequence, toward a limit g ∈ L 2 (∂Ω N ). Solutions to the elastic problem with Neumann boundary conditions are equal up to a rigid translation. The set of equivalence classes is denoted H 2 (Ω 0 )/T (Ω 0 ). By adding the condition ´Ω0 u dx = 0, we select the class representative with the least H 2 -norm.

W M = g ∈ H 1/2 (∂Ω 0 ) ∥g∥ H 1/2 (∂Ω 0 ) ⩽ M and ∀w ∈ R(Ω 0 ) ˆ∂Ω 0 g • w ds = 0 .
In other words, the state u j ∈ H 2 (Ω 0 ) associated to g j satisfies

∥u j ∥ H 2 (Ω 0 ) = min w∈R(Ω 0 ) ∥u j + w∥ H 2 (Ω 0 ) = ∥u j ∥ H 2 (Ω 0 )/R(Ω 0 ) , (4.9)
where notation is defined in (3.16). Due to Proposition 3.14, the sequence (u j ) is uniformly bounded in H 2 (Ω 0 ). Due to the continuous imbedding H 2 (Ω 0 ) → W 1,6 (Ω 0 ), (u j ) is also bounded in W 1,6 (Ω 0 ), and it converges uniformly, up to a subsequence, toward a limit u. Therefore J(u j ) ---→ j→∞ J(u).

We check that u is a solution associated to g. As the sequence (u j ) converges uniformly, it also converges weakly in L 2 (Ω 0 ). Now, as (u j ) is bounded in H 1 (Ω), it converges H 1 -weakly, up to a subsequence, toward u. This last remark, together with the elasticity equation and the weak convergence of (g j ) in L 2 (∂Ω 0 ), yields

∀v ∈ H 1 (Ω 0 ) ˆΩ0 Aε(v) : ε(u) dx = lim j→∞ ˆΩ0 Aε(v) : ε(u j ) dx = lim j→∞ ˆ∂Ω 0 v • g j ds = ˆ∂Ω 0 v • g ds.
Finally, u is a state associated to g and g is a solution of problem (4.8).

Optimality conditions

In this paragraph we derive optimality conditions for the optimal control problem (2.13).

When J is Gateaux differentiable, it is standard to use an adjoint state p to transform descent directions in the space of states into descent directions in the space of controls. The adjoint state associated with a state u where J is differentiable is solution to the adjoint problem

∀q ∈ H 1 D (Ω 0 ) ˆΩ0 Aε(p) : ε(q) dx = ⟨dJ(u), q⟩ . (4.10) If dJ(u) is in H -1 D (Ω 0 )
, then, due to the Lax-Milgram theorem, p is unique and belongs to H 1 D (Ω 0 ). In the following optimality conditions, we make the hypothesis dJ(u) ∈ H -1 D (Ω 0 ), keeping in mind that optimality conditions are useful for the numerical application, where everything is much simpler (and finite-dimensional).

In a similar way to Lemma 4.3, we consider the case of minimizers for the unconstrained problem (2.11), where the objective function is necessarily differentiable. The corresponding optimality conditions are stated in the following proposition.

Proposition 4.3. Let g ∈ L ∞ (∂Ω N ) a local (unconstrained) minimizer of the application F : g → J(u g ) + R(g). Then F is differentiable at g, and g satisfies the first-order optimality condition

∀h ∈ L ∞ (∂Ω N ) ˆ∂Ω N p • h ds + ⟨dR(g), h⟩ = 0,
where p is the adjoint state defined by (4.10).

Proof. For g ∈ L ∞ (∂Ω N ) a control, the directional derivative of F in a direction h ∈ L ∞ (∂Ω N ) reads dF (g)(h) = dJ(u g )(w h ) + ⟨dR(g), h⟩,
where

w h ∈ H 1 D (Ω 0 ) is the state associated to h, defined by ∀v ∈ H 1 D (Ω 0 ) ˆΩ0 Aε(w h ) : ε(v) dx = ˆ∂Ω N h • v ds. (4.11)
Now, if g is a minimizer of F , then for every direction h, g satisfies dF (g)(h) ⩾ 0. Assume that there is a direction h such that dF (g)(h) > 0. Then, using Lemma 4.2, we obtain dF (g)(-h) = dJ(u g )(-w h ) -⟨dR(g), h⟩ ⩽ -dJ(u g )(w h ) -⟨dR(g), h⟩ < 0.

In other words, -h is a descent direction, which contradicts g being a local minimizer. As a consequence, a local minimizer g satisfies dF (g)(h) = 0 for every direction h ∈ L ∞ (∂Ω N ). Therefore, F is differentiable at g with dF (g) ≡ 0.

In particular, J is differentiable at u g . We combine (4.11) and (4.10) to obtain an expression of dJ(u) as a function of h,

⟨dJ(u), w h ⟩ = ˆΩ0 Aε(p) : ε(w h ) dx = ˆ∂Ω N p • h ds.
Finally, the differential of F reads

⟨dF (g), h⟩ = ⟨dJ(u), w h ⟩ + ⟨dR(g), h⟩ = ˆ∂Ω N p • h ds + ⟨dR(g), h⟩ = 0,
hence the optimality condition.

The proof of Proposition 4.3 relies on the fact that every direction h ∈ L ∞ (∂Ω N ) is an admissible direction. In particular if h is an admissible direction, then -h is also an admissible direction. When the optimization problem involves a L ∞ constraint, the cone of admissible directions is not the whole space anymore, in particular when the constraint is active. As a consequence, we cannot expect J to be Gateaux differentiable at local minimizers for the constrained problem.

When J is not differentiable, J can be approximated at first order using several linear forms ℓ (defined below), associated to families (x y ) y∈Γ (see Section 4.1 for more details). Each linear form approximates J on a cone of directions. Those multiple linear forms are represented in the space of controls by as many adjoint states p ℓ , which satisfy

∀q ∈ H 1 D (Ω 0 ) ˆΩ0 Aε(p ℓ ) : ε(q) dx = ⟨ℓ, q⟩ = ˆΓ q(x y ) • (x y + u(x y ) -y) dy. (4.12)
For a given direction h in the space of controls, there is a linear form ℓ 0 associated to a family (x y,0 ), such that

dJ(u g )(w h ) = ⟨ℓ 0 , w h ⟩ = ˆ∂Ω N p ℓ 0 • h ds, (4.13)
where w h is the corresponding direction in the space of displacements.

In the following proposition, we consider the constrained problem (2.13), where R(g) = α 2 ∥g∥ 2 L Theorem 4.4. Let g ∈ G M a local minimizer of problem (2.13) and u g the associated state defined by (2.6). For (x y ) y∈Γ a family such that ∀y ∈ Γ, x y ∈ P y (u g ), we consider the linear form ℓ and the adjoint state p ℓ defined in (4.12). Then there exists a Lagrange multiplier

λ ℓ ∈ L 2 (∂Ω N , R) with for a.e. x ∈ ∂Ω N λ ℓ (x) = 0 if ∥g(x)∥ < M λ ℓ (x) ⩾ 0 if ∥g(x)∥ = M
such that g satisfies the first-order optimality condition for a.e. x ∈ ∂Ω N p ℓ (x) + (α + λ ℓ (x))g(x) = 0. (4.14)

Proof.

• We first give a justification for the definition of the adjoint state. We use the notation

F (g) = J(u g ) + α 2 ∥g∥ 2 L 2 (∂Ω N ) . For a direction h ∈ L ∞ (∂Ω N )
, we consider the associated displacement field w h , solution of (4.11). By combining (4.11) with (4.12), we obtain

ˆΓ w h (x y ) • (x y + u(x y ) -y) dy = ˆΩ0 Aε(p ℓ ) : ε(w h ) dx = ˆ∂Ω N p ℓ • h ds.
The derivative of F in the direction h satisfies dF (g)(h) = ˆΓ min

x∈P y (u) w h (x) • (x + u(x) -y) dy + α ˆ∂Ω N g • h ds ⩽ ˆΓ w h (x y ) • (x y + u(x y ) -y) dy + α ˆ∂Ω N g • h ds = ˆ∂Ω N (p ℓ + αg) • h ds.
• We now characterize the cone of admissible directions at g. Admissible directions h ∈ L ∞ (∂Ω N ) are only constrained at points where the L ∞ constraint is active. We split ∂Ω N into two disjoint subsets, the active set A and the inactive set I, so that A ∪ I = ∂Ω N . The two subsets are defined by

A = {x ∈ ∂Ω N | ∥g(x)∥ = M } and I = {x ∈ ∂Ω N | ∥g(x)∥ < M } . (4.15)
Admissible directions h satisfy the condition for a.e. x ∈ A h(x) • g(x) ⩽ 0.

To be more specific, the cone of admissible directions is the intersection L ∞ (∂Ω N )∩ K, where K is the closed cone

K = h ∈ L 2 (∂Ω N ) ∀λ ∈ L 2 (∂Ω N , R + ) ˆA λ(h • g) ds ⩽ 0 . Another expression for K ⊂ L 2 (∂Ω N ) reads K = h ∈ L 2 (∂Ω N ) ∀q ∈ Q ⟨h, q⟩ L 2 (∂Ω N ) ⩽ 0 , where Q = 1 A λg, λ ∈ L 2 (∂Ω N , R + ) .
Here, Q is the polar cone of K (Rockafellar, 1970, Section 14). As Q is convex, the orthogonal projection

s Q of -(p ℓ + αg) onto Q satisfies ∀q ∈ Q ⟨(p ℓ + αg) + s Q , q⟩ = -⟨-(p ℓ + αg) -s Q , q -s Q ⟩ ⩾ 0.
Therefore, -(p ℓ + αg) -s Q belongs to K, and we obtain the orthogonal decompo-

sition (p ℓ + αg) = -s Q -s K , with s Q ∈ Q and s K ∈ K.
• To derive the first-order optimality condition, we now consider the admissible direction

h = s K max(1, ∥s K ∥) ∈ L ∞ (∂Ω N ) ∩ K.
As g is a local minimizer, it results from Euler's inequation that

0 ⩽ dF (g)(h) ⩽ ˆ∂Ω N (p ℓ + αg) • h ds = - ˆ∂Ω N min ∥s K ∥, ∥s K ∥ 2 ds.
Finally, we obtain s K = 0, which yields the first-order optimality condition p ℓ + αg + s Q = 0.

We now turn to the problem with subdomain restriction. For the sake of simplicity, we assume that J is Gateaux differentiable at the considered local minimizer. However, we keep in mind that, in the case of a non-differentiable minimizer, the derived optimality condition applies to every adjoint state p ℓ . Problem (2.14) involves a differentiable constraint on the state. We manage this state constraint using another adjoint state z ∈ H 1 D (Ω 0 ), defined by

∀q ∈ H 1 D (Ω 0 ) ˆΩ0 Aε(z) : ε(q) dx = ˆω0 u • q dx. (4.16)
Though we derive optimality conditions for this problem, the state-constrained approach is not considered in the following chapters.

Theorem 4.5. Let g ∈ L ∞ (∂Ω N ) a local minimizer for problem (2.14) and let u g the associated state defined by (2.6). Assume also that J is differentiable at u g with dJ(u g ) ∈ H -1 D (Ω 0 ). Then there exists a Lagrange multiplier λ ∈ L 2 (∂Ω N , R) with for a.e.

x ∈ ∂Ω N λ(x) = 0 if ∥g(x)∥ < M λ(x) ⩾ 0 if ∥g(x)∥ = M
and a real multiplier η ∈ R + such that g satisfies the first-order optimality condition for a.e. x ∈ ∂Ω N p(x) + λ(x)g(x) + ηz(x) = 0. (4.17)

Moreover, one has either η = 0 or ´ω0 ∥u∥ 2 dx = U .

Proof. We use the notation F (g) = J(u g ), C(g) = ´ω0 ∥u g ∥ 2 dx, and for a direction h ∈ L ∞ (∂Ω N ), w h is the associated state defined by (4.11). Due to their definitions, the adjoint states p and z satisfy

∀h ∈ L ∞ (∂Ω N ) ⟨dF (g), h⟩ = ˆ∂Ω N p • h ds and ⟨dC(g), h⟩ = ˆ∂Ω N z • h ds
We also define A and I as previously.

If C(g) < U , i.e. the constraint on u g is inactive, the same procedure as above results in (4.14), which is (4.17) with η = 0. Now we consider the case

C(g) = U . Admissible directions h satisfy h • g ⩽ 0 in A and ´∂Ω N z • h ds ⩽ 0. It is easy to check that such directions h belong to L ∞ (∂Ω N ) ∩ K,
where the closed cone K ⊂ L 2 (∂Ω N ) is defined by

K = h ∈ L 2 (∂Ω N ) ∀λ ∈ L 2 (∂Ω N , R + ) ˆA λ(h • g) ds ⩽ 0 and ˆ∂Ω N z • h ds ⩽ 0 .
Let us check that the polar cone to K reads

Q = 1 A λg + ηz, λ ∈ L 2 (∂Ω N , R + ), η ∈ R + . Consider q = 1 A λg + ηz ∈ Q, with λ ∈ L 2 (∂Ω N , R + ) and η ∈ R + . Then for every h ∈ K, ⟨h, q⟩ L 2 (∂Ω N ) = ˆA λ(h • g) ds + η ˆ∂Ω N z • h ds ⩽ 0.
Now consider a direction h ∈ L 2 (∂Ω N ) such that ⟨q, h⟩ L 2 (∂Ω N ) ⩽ 0 for all q ∈ Q. By considering q 1 = 1 A λg ∈ Q and q 2 = z ∈ Q, we obtain h ∈ K. Therefore, for h ∈ L 2 (∂Ω N ), the equivalence

h ∈ K ⇔ ∀q ∈ Q ⟨h, q⟩ L 2 (∂Ω N ) ⩽ 0 holds and Q is the polar cone to K. As Q is convex, the orthogonal projection p Q of -p onto Q satisfies ∀q ∈ Q ⟨p + p Q , q⟩ = -⟨-p -p Q , q -p Q ⟩ ⩾ 0, which results in the unique orthogonal decomposition p = -p Q -p K with p Q ∈ Q and p K ∈ K. If we consider the admissible direction h = p K max(1,∥p K ∥) , Euler's inequation results in 0 ⩽ ˆ∂Ω N h • p ds = - ˆ∂Ω N min ∥p K ∥, ∥p K ∥ 2 ds.
As a consequence, p K = 0 and p + p Q = 0, which is the optimality condition we were looking for.

Chapter 5

An adjoint method to solve the registration problem

To solve the registration problem numerically, we discretize the optimal control formulation using the finite element method. The finite element method is relevant to handle domains with all kinds of shapes. For this reason it is widely used in the Mimesis team, especially through the home-brewed SOFA framework [START_REF] Allard | SOFA -an Open Source Framework for Medical Simulation[END_REF]. For all numerical investigations, we adopt a discretize-then-optimize approach, which means that the entire formulation is transformed into a finite-dimensional problem, and then solved using numerical tools for finite-dimensional optimization. The reasons for this choice are multiple. First, as we do not change or refine meshes during the online procedure, the problem dimension is fixed and going back to continuous considerations is not necessary. In addition, using expressions from the continuous framework may represent a lack of robustness, as the solver may fail because the discretized problem does not approximate the continuous problem well enough. We really need our numerical methods to be as simple as possible, and forgetting that the liver mesh actually represents a real liver seems to avoid unnecessary considerations, at the expense of continuous consistency.

In this chapter, we give some details about our implementation of an adjoint method. This choice results from our initial intention to use shape optimization. Adjoint methods are also relevant in optimal control when the direct problem is handled by a black-box software such as SOFA, and for this reason we kept using an adjoint method. Modularity is actually a key advantage of the adjoint method, as it can be broken down into several procedures, each of which can be handled by a separate protagonist. In our case, these procedures consist in evaluating the registration functional, managing the elastic model, and solving the optimization problem, respectively.

After defining notations to express the registration problem in terms of matrix and vector operations, we present an overview of the adjoint method and its structure. In the end of the chapter we consider separately each procedure involved in the adjoint pipeline. We first describe the chosen option to perform orthogonal projections onto the mesh boundary when evaluating the cost function. Then, we mention Newton methods that can be used to solve nonlinear elastic problems. Finally, we discuss the choice of an optimization solver and show some examples to illustrate their convergence properties.

Finite element discretization of the problem

In order to perform numerical computations, we derive a discretized version of the problem using the finite element method. The domain Ω 0 is represented by a tetrahedral mesh T 0 , and functions of interest are represented by P1 finite element functions defined on T 0 . When a deformation is applied to the mesh T 0 , the deformed mesh is known through the position of its vertices x 1 , • • • , x n ∈ R 3 . We denote by x = (x 1 , • • • , x n ) the vector containing the coordinates of all vertices of the deformed mesh, and x 0 the positions vector in the reference configuration. We also use u

= x -x 0 = (u 1 , • • • , u n )
to denote the displacement of vertices with respect to the initial configuration. With a slight abuse of notation, we use the same bold letter u to denote the P1 finite element function associated to the displacement field in T 0 . As a P1 finite element function is defined by its value at the mesh vertices, a finite element function and a vector represented by the same letter actually contain exactly the same information.

In the discrete framework, we write the linear elasticity equation (3.15) using the matrix formulation Au = Sg, (5.1)

where the force distribution g is a P1 finite element function defined only on the boundary elements of ∂Ω N . If we denote by I g the set of indices associated to vertices of ∂Ω N , one can write g = (g i ) i∈I g , resulting in a vector g smaller than u. Here, the stiffness matrix A stands for the linear elasticity inner product and satisfies

v T Au = ˆΩ0 Aε(u) : ε(v) dx,
while the matrix S represents the L 2 inner product on ∂Ω N , which satisfies the equality

u T Sg = ˆ∂Ω N u • g ds.
Here, A and S were obtained by expressing the involved bilinear forms in a basis of the finite element space. In the context of nonlinear elasticity, a similar development involving the same basis leads to a nonlinear system of the form

F(u) = Sg.
We give more details about this system at the end of the chapter. Though we produced an early implementation where we control the surface force distribution g, we proceed in a different way in our main implementation, for two reasons. First, choosing g as an optimization variable would require to perform several products between S and a vector. Namely, at least two products involving S or its adjoint per evaluation of the objective function would be necessary, whereas in a real-time context, it is preferable to make computations as light as possible. In addition, the SOFA framework does not provide an implementation of S. As a consequence, we control the nodal forces instead, which is consistent with the habits of the Mimesis team. The nodal forces are represented by the vector b

= Sg = (b 1 , • • • , b n ). An element b k ∈ R 3 can
be interpreted as the force (in Newtons) that applies on the k-th node, so that the work of the force distribution (also known as compliance), reads

ˆ∂Ω N u • g ds = u T b = n i=1 u i • b i .
In this context, the Neumann boundary ∂Ω N is defined in terms of nodes where we allow b k to be nonzero. In particular, b k is zero for every node that is not on the mesh boundary. The state equation in the linear case and in the nonlinear case reads

Au = b and F(u) = b,
respectively. From now on, in the manuscript, we use the b-formalism, even to describe optimization methods from the early implementation that are implemented using the g-formalism.

In the context of the adjoint method, the optimization procedure is performed in the space of nodal forces vectors b. The inner product used by generic optimization solvers to compute descent directions and measure optimality conditions is the standard inner product in R 3n defined by ⟨u, v⟩

= u T v = u i • v i .
As a consequence, in this chapter gradients and derivatives are always expressed using this inner product, with the notation

⟨df (u), v⟩ = ∇f (u) T v = n i=1 ∂ i f (u) • v i , (5.2)
where f is a scalar function.

A consequence of this choice is a less significant consistency with the continuous optimal control problem, as the space of nodal forces is not a finite element space, and the vector of nodal forces b has no counterpart in the continuous problem. In addition, as opposed to inner products in function spaces, the standard inner product in R 3n is completely blind to geometric information such as the size of triangles. In other words, we are not even trying to perform iterations that approximate iterations in the continuous case. However, every standard optimization solver uses the standard inner product in the discrete vector space. To be consistent with the continuous case, we would need to implement optimization solvers that work with the L 2 inner product, at the cost of several additional matrix-vector products and system inversions. Finally, working with nodal forces at the expense of consistency with the continuous framework is relevant in view of real-time requirements and available software.

Adjoint method

The adjoint method is a standard approach for solving optimal control problems, based on the optimality conditions derived in Section 4.3. In particular, an adjoint state p is used to compute the derivatives of the objective function with respect to the control. A general presentation of the adjoint method is given in Section 5.2.1.

In an adjoint method, the direct problem and the optimization problem are handled in separate procedures, which makes the resulting program very modular. In Section 5.2.2, we give a few details about our implementation of the adjoint method.

Overview of the method

The adjoint method is very specific to optimal control problems, where the entire system state is defined by the choice of a control. We consider the discrete version of the optimal control problem ( 2 (5.4)

Here, the elastic system is hidden in the objective function Φ. In particular, each evaluation of the objective value Φ(b) requires to solve the elastic system to compute u b . Now, to solve (5.4) using a first-order optimization method, computing the objective gradient ∇Φ(b) is also required. The derivation is similar to that in Section 4.3. If h is an admissible direction in the space of controls, the corresponding direction in the space of states w is solution of the tangent system F ′ (u b )w = h, where F ′ denotes the Jacobian of the elastic residual. With this notation, the first-order variation of Φ reads

∇Φ(b) T h = ∇J(u) T w + ∇R(b) T h.
The adjoint state p is used to transform the functional gradient ∇J(u b ) in the space of states into a gradient in the space of controls. It is defined as the solution to the adjoint system F ′ (u b ) T p = ∇J(u b ). Note that here the constraint Jacobian F ′ (u b ) is a symmetric matrix, as it represents the Hessian of the discretized elastic energy. We obtain

∇J(u b ) T w = p T F ′ (u b )w = p T h,
and, finally,

∇Φ(b) = p + ∇R(b).
This gradient can then be used by an optimization solver to minimize Φ. Algorithm 1 gives a summary of the adjoint procedure to compute ∇Φ(b). Problem (5.4) can be fed to generic first-order primal solvers, such as gradient descents or quasi-Newton solvers. By using solvers for constrained optimization problems, it is possible to handle constraints on b. If the optimization problem involves constraints on u, however, calculations similar to those above are necessary to transform those constraints into constraints on b. Remark 5.1. In the space of displacements, the first-order descent direction for J returned by the adjoint method reads

w = F ′ (u) -1 F ′ (u) -T ∇J(u),
Though it ensures the regularity and feasibility of displacement fields, this expression seems less natural than ∇J(u) as a descent direction for J, due to the change of metrics created by the elastic model. If the Jacobian F ′ (u) is poorly conditioned, using an adjoint method can degrade the numerical properties of the problem, making convergence more difficult.

As mentioned earlier, implementing an adjoint method is very convenient, as the direct solver and the optimization solver are separated from each other. In the next section, we briefly describe the software produced during the thesis.

A modular implementation of the adjoint method

We first produced an early implementation of the adjoint method using the FreeFem++ software [START_REF] Hecht | New development in FreeFem++[END_REF]. This implementation was used to test the adjoint method on twodimensional examples, with homemade procedures and optimization algorithms. Then we switched to Python, in order to benefit from the features of a modern language, such as object oriented programming. It was also the occasion to use state-of-the-art scientific libraries to solve larger-scale problems in a reasonable time. As mentioned above, this early implementation works with the surface force finite element function g as the control variable, while the Python implementation uses b as the control variable. For the sake of consistency, we always consider that b is the control in the description of the implemented procedures.

Our Python implementation involves three separate procedures that communicate together through simple interfaces, using the NumPy linear algebra library. Figure 5.1 shows an overview of the adjoint procedure to evaluate the shape functional and its gradient.

On the one hand, the shape functional procedure is responsible for evaluating the value and gradient of the functional J for a given displacement field u. It is expected to return an accurate gradient, as otherwise the optimization solver would fail to solve the problem with a tight tolerance. This part of the program is implemented using the Trimesh package and involves an exact projection onto the deformed triangular surface of the discretized organ. We give details about this procedure in Section 5.3.

On the other hand, the optimization solver is responsible for making the objective function decrease. In our implementation, only first-order methods are supported, as our shape functional procedure does not provide second-order information. We use the standard solvers available in the Scipy package. In particular, quasi-Newton solvers are the reference choice when only first-order information is available. The choice of an optimization solver is discussed in Section 5.5.

Remark 5.2. Note that using second order solvers such as the Newton method would require careful attention, as the objective Hessian ∇ 2 Φ is full due to the presence of the elastic problem. In addition, even if the Hessian is not assembled, performing Hessianvector products might require to solve several linear systems, resulting in a heavy computational cost.

Finally, the central element in the implementation is the finite element package, that manages the elastic problem. The finite element package is responsible for connecting the optimization solver, that works in the space of controls, with the shape functional procedure, that works in the space of states. This part of the program is critical for the success of the registration process. First, it should give accurate results for the same reasons as the shape functional procedure. But most of all, the direct solver should not fail, in particular when it involves an iterative algorithm. When a linear elastic model is used, the stiffness matrix is precomputed using Fenics or SOFA, then converted into a Scipy matrix and factorized before the registration starts. Thus, the direct and adjoint problem are solved in a fast and robust way. However, when the elastic model is nonlinear, the direct problem is solved iteratively using a Newton method. In this last case, enforcing convergence requires special attention, as discussed in Section 5.4.

Discretized shape functional

In this section, we describe a procedure to evaluate the shape functional J numerically. In the discrete framework, the deformed surface ∂Ω u is a triangular mesh, and defining S 0 consists in selecting a set of triangles that are used to evaluate J. The observation Γ is provided as the point cloud

Γ = {y 1 , • • • , y p } ⊂ R d .
Due to this lack of structure, the objective function cannot be defined as an integral on Γ. Thus, we use the expression

J(u) = 1 2p p i=1 d 2 (y i , S u ).
(5.5)

Optimization solver Finite element package Functional

New With this notation, 2J(u) is the quadratic mean of distances from a point in Γ to the deformed mesh surface.

Beside the fact that the integral is transformed into a sum, the discrete version of J is very close to the continuous one, as the points from Γ are projected exactly onto S u . The projection procedure exploits the triangular mesh structure of S u . It involves a nearest neighbor search that determines the closest triangle from y. As one nearest neighbor search is performed per point in Γ, its implementation should be efficient. We discuss such procedures at the end of this section.

First, we briefly explain how the derivative formula computed in Lemma 4.1 adapts to the finite element configuration.

Evaluation of the functional gradient

We now choose a point y ∈ Γ and we derive a practical formula to evaluate the gradient of the elementary application j : u → 1 2 d 2 (y, S u ).

(5.6)

The expression d(y, S u ) denotes the distance between y and its projection a onto S u . First, note that here the projection point a is supposed to be unique. Actually, even when y has several projection points onto S u , the nearest neighbor search algorithm is supposed to return only one of them. This detail may be a problem when it comes to obtaining reproducible results, as the projection point depends on the search procedure. However, cases with multiple projections are rare, as for almost every displacement field u, each point in Γ has a unique projection. In addition, a case with multiple projections means that several descent directions are available, and by choosing one of these projections, the algorithm selects one of those directions. Thus, the optimization procedure is not likely to fail because of that.

Since we are dealing with P1 finite element functions, differentiating j with respect to the displacement field u means differentiating j with respect to the displacements (u 1 , • • • , u n ) at the mesh vertices. Thus, the gradient ∇j(u) is defined by its vertex-wise components (∂ 1 j(u), • • • , ∂ n j(u)). We assume that the projection point a falls into the triangle x k 1 x k 2 x k 3 . Provided that a is the only projection point, the nodes x k 1 , x k 2 and x k 3 are the only ones whose displacement has an impact on the first-order variation of j. Indeed, if another mesh node moves, the triangle containing a is not affected and the value of j remains constant at first order. Therefore,

∀k / ∈ {k 1 , k 2 , k 3 } ∂ k j(u) = ∂j ∂u k (u) = 0.
We now apply Lemma 4.1 to determine the nonzero components of ∇j(u). We denote by θ 1 , θ 2 , θ 3 the barycentric coordinates of a in the triangle

x k 1 x k 2 x k 3 , so that a = θ 1 x k 1 + θ 2 x k 2 + θ 3 x k 3 , where θ 1 , θ 2 , θ 3 ⩾ 0 and θ 1 + θ 2 + θ 3 = 1. If v is a perturbation of the displacement field u, Lemma 4.1 yields ⟨dj(u), v⟩ = v(a) • (a -y) = θ 1 v k 1 + θ 2 v k 2 + θ 3 v k 3 • (a -y) .
This expression confirms that ∇j(u) has at most three nonzero components. Finally, the components of ∇j(u) read

∀i ∈ {1, 2, 3} ∂ k i j(u) = θ i (a -y) ∀k / ∈ {k 1 , k 2 , k 3 } ∂ k j(u) = 0.
(5.7) Formula (5.7) contains the cases where a falls onto an edge or a vertex. These cases appear when one or two of the barycentric coefficients are set to zero, respectively. Figure 5.2 illustrates the gradient formula in a two-dimensional case. On the left, y lies in the polar cone to S u at x k and its projection falls onto x k , resulting in a gradient fully supported by the node k. On the right, a belongs to a boundary edge and ∇j(u) has nonzero components at x k and x m , which are proportional to the corresponding barycentric coefficients. We just computed the gradient of the elementary function j i associated to a single point y i ∈ Γ. By summing all the elementary gradients, we obtain the global gradient of the functional

∇J(u) = 1 p p i=1 ∇j i (u).
(5.8)

Remark 5.3. Computing the first-order derivatives of j using automatic differentiation is very convenient, as the first-order derivatives of j are the same as the derivatives of the application

u → 1 2 ∥θ 1 x k 1 + θ 2 x k 2 + θ 3 x k 3 -y∥ 2 , -∂ k j(u) x k x ℓ x m y -∂ k j(u) -∂ m j(u) x k x ℓ x m y a Figure 5
.2: Illustration of ∇j(u) in two configurations. Points x ℓ , x k and x m are located on the mesh boundary, while the last vertex is inside the mesh. The red arrows represent the components of the descent direction -∇j(u). On the right, the equality y -a = -∂ k j(u) -∂ m j(u) holds.

with the convention x k = x 0 k + u k . In other words, the barycentric coefficients can be treated as constants in the gradient computation. For instance, assume that the nearest neighbor search procedure returns the indices k1, k2, k3 and the corresponding barycentric coordinates t1, t2, t3. Then, the elementary gradient ∇j(u) can be computed using, for instance the following lines in PyTorch: However, this remark is not true when it comes to evaluating the second-order derivatives of j. In this case, it is also necessary to save the operations performed on x to obtain θ 1 , θ 2 , θ 3 , as they depend on the displacement field.

Computation of orthogonal projections

We now discuss the nearest neighbor procedure that is used to find the nearest triangle from a point y. Computing J(u) requires to compute the orthogonal projection of each point of Γ onto S u , which means that the search procedure is run once per point in Γ. As a consequence, efficiency in this part of the program is critical for the global performance of the algorithm.

Though we did not investigate quantitatively the performance of orthogonal projections, we give a few details about the procedures involved in our code. A naive approach would consist in evaluating the distance between each point in Γ and each triangle in S u , which leads to a complexity of O(p × n tri ), where p = |Γ| and n tri is the number of triangles.

Except for the naive approach, we did not implement the procedures below ourselves. We used existing functions available in Python packages.

Signed distance field: a continuous approach We first considered an approach involving a signed distance field with respect to the deformed organ. Here, we only consider the case S 0 = ∂Ω 0 . Remember that, whenever the projection of a point y is unique, it can be expressed as a function of the signed distance field d Ω u by the expression

p ∂Ω u (y) = y -d Ω u (y)∇d Ω u (y).
(5.9)

We implemented this method using the MshDist software [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF], resulting in faster projections than the naive method, especially in three-dimensional cases. The procedure begins by evaluating the signed distance function d Ω u at the vertices of a background tetrahedral mesh, so that the signed distance and its gradient can be interpolated a each point y ∈ Γ. Then the coordinates of the orthogonal projection are computed using (5.9). Finally, some interpolation work is required to find the triangle containing the projection point. Figure 5.3 shows a signed distance field computed on a background mesh in a bi-dimensional case. An advantage of this approach is that the most expensive operation has to be executed once for all points in Γ. Though, the signed distance field has to be recomputed every time u changes. The signed distance approach comes with several other drawbacks. First, it is necessary to use a background mesh which should always contain the organ model while it is subject to deformations. If the deformed organ does not stay inside the background mesh, the procedure may fail. In addition, the signed distance field is approximated using a P1 finite element function. For the signed distance estimation to be accurate, the background mesh needs to be dense, which makes the computation costly. However, even with a dense background mesh, the signed distance gradient may be inaccurate, especially close to the skeleton of ∂Ω u . In the same vein, the signed distance value may be inaccurate close to ∂Ω u , making it impossible to solve the optimization problem with a tight tolerance. For those reasons, despite its apparent interest from a continuous point of view, an approach based on the signed distance is not adequate for our needs.

Spatial indexing structures

While the signed distance approach is based on the discretization of a continuous process, the solution we chose is based on the polyhedral structure of the discretized organ and explicitly consists in finding the closest triangle from y. The procedure is based on spatial indexing structures. Such structures are used for instance in computer graphics for collision detection [START_REF] James | BD-Tree: Output-Sensitive Collision Detection for Reduced Deformable Models[END_REF].

We use a procedure available in the Trimesh Python package1 , based on a R-tree [START_REF] Guttman | R-Trees: A Dynamic Index Structure for Spatial Searching[END_REF]. The functioning of such a structure is well beyond the scope of this Figure 5.3: Signed distance field with respect to a two-dimensional mesh. Projection directions are orthogonal to the level curves of the signed distance field.

manuscript, but here is the main idea. In an R-tree, triangles are gathered in small groups represented by their bounding boxes. Those small boxes are gathered in bigger boxes, which are gathered in even bigger boxes, and so on. This results in a hierarchy of boxes, represented by a tree structure. When searching for the closest triangle to a point y ∈ Γ, the algorithm does not need to explore smaller boxes if they are contained in a bigger box that is already too far from y. This results in an average query complexity of O(log n tri ). Figure 5.4 shows an example of such a box hierarchy on a surface mesh. The performance of an R-tree is determined by the insertion algorithm used to build the tree. More complex insertion algorithms are more costly but result in more efficient queries. A popular variant is the R*-tree [START_REF] Beckmann | The R*-Tree: An Efficient and Robust Access Method for Points and Rectangles[END_REF], where the area and overlapping of boxes is minimized. Later, packing algorithms [START_REF] Leutenegger | STR: a simple and efficient algorithm for R-tree packing[END_REF][START_REF] García | A Greedy Algorithm for Bulk Loading R-Trees[END_REF] have proven efficient to add several elements into the tree in one procedure.

The R-tree approach is particularly relevant in our application case. First, the procedure returns the triangle on which the projected point falls, along with barycentric coordinates. In addition, it is possible to select the triangles that are candidates to match the observed point cloud, when a signed distance field can only be computed with respect to the full shape.

Let us add a qualitative remark about the performance of the R-tree approach. In this manuscript, we use the algorithm from Trimesh as a black box but we are far from an efficient implementation. At each evaluation, a new R-tree is built and then one query per point in Γ is performed. As the mesh topology does not change along the deformation process, it is not necessary to build a new tree from zero before each evaluation of the Brute-force returns on the GPU Due to the real-time requirements, we are not expecting to deal with very dense meshes or point clouds. While spatial indexing structures are able to store millions of points, we only expect to compare at most a few thousands points with a few thousands triangles. In this case, evaluating the distances between every point and every triangle is still possible in a reasonable time, due to the highly parallel structure of graphic cards (or GPU for Graphics Processing Unit). In particular, this brute force approach is implemented in a function from the PyTorch3D library and might be considered for a GPU implementation of the registration procedure. Also, note that implementations of spatial indexing structures on the GPU exist in the literature [START_REF] Prasad | GPU-based Parallel R-tree Construction and Querying[END_REF].

Newton methods for static elasticity problem

We now turn to the procedure that manages the elastic model using a finite element discretization. Namely, this procedure is responsible for solving the direct problem to transform a force distribution b into a displacement field u, and to solve the adjoint problem to transform a gradient ∇J(u) in the space of displacements into an adjoint state p. When a linear model is used, the role of the finite element package boils down to solving linear systems involving the stiffness matrix. Even in the nonlinear case, the adjoint problem remains linear, and an adjoint states is computed by assembling and solving a linear system. For this reason, we focus here on the most challenging part, which is to solve the direct problem in the nonlinear case. This task is usually taken care of using a Newton method. Though, some precautions must be taken to ensure the convergence of the algorithm.

Invented in the 17th century, the Newton method (also known as the Newton-Raphson method) remains nowadays the go-to method to solve nonlinear variational problems such as static elasticity (Kelley, 1995, Chapter 5 and references therein). In a recent study, [START_REF] Morch | A review of higher order Newton type methods and the effect of numerical damping for the solution of an advanced coupled Lemaitre damage model[END_REF] compare the performance of Newton-type methods of order 2 and higher on a few selected problems, and evaluate the effect of damping on the method robustness. The most computationally efficient method in their benchmark turns out to be the standard Newton method, even though other variants exhibit a higher convergence rate in theory. Other approaches, such as quasi-Newton methods [START_REF] Gelin | Use of quasi-Newton methods for large strain elastic-plastic finite element computations[END_REF] show interesting performance in mechanical simulation [START_REF] Liu | Quasi-Newton Methods for Real-Time Simulation of Hyperelastic Materials[END_REF][START_REF] Yusa | Performance investigation of quasi-Newton-based parallel nonlinear FEM for large-deformation elastic-plastic analysis over 100 thousand degrees of freedom[END_REF]. However, we focus here on the Newton method and some of the remarks we formulate also apply to quasi-Newton methods.

In Section 2.2, we formulated the nonlinear elastic problem as the minimization problem (2.1). Here, we consider a discretized version of (2.1). After defining some notations, we describe the standard Newton algorithm and a variant adapted to optimization problems.

Remember that the discretized displacement field is represented by a vector u ∈ R 3n , where n is the number of vertices in the mesh, and forces are represented by b ∈ R 3n . The discrete optimization problem reads

min u∈U W (u) -b T u, (5.10)
where we keep using U to denote the feasible set in the discrete case, and we denote the first-order optimality conditions by F(u) -b = 0 where F = ∇W.

(5.11)

Newton method for variational problems

The most common approach to solving nonlinear elastic problems is to solve the optimality condition (5.11) as a variational problem. In this context, the Newton method takes the form of a fixed-point algorithm, where a linear approximation of (5.11) is solved at each iteration. The first-order development of the residual reads

F(u + w) = F(u) + K(u)w + o(∥u∥),
where K(u) = F ′ (u) denotes the residual Jacobian at u (also called the stiffness matrix). The Newton step w is the solution of the tangent system

K(u)w = b -F(u),
Algorithm 2: Standard Newton method for nonlinear variational systems. Data: Initial iterate u 0 , tolerance ϵ for k = 0, 1, 2, . . . do Evaluate the residual

r k = b -F(u k ) if ∥r k ∥ ⩽ ϵ then return u k end Compute a step w k by solving K(u k )w k = r k u k+1 ← u k + w k end
and the next iterate is defined as u + = u + w.

The main source of cost in the method is the inversion of a linear system involving K(u) at each iteration. This operation can either be taken care of by a direct solver, for instance involving a LU factorization (Nocedal and Wright, 2006, Appendix A), or by an iterative solver such as GMRES [START_REF] Saad | GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems[END_REF]. Though the principles of such solvers are out of the scope of this manuscript, we briefly describe the advantages of both families of solvers in Newton methods. Iterative solvers are more relevant for large-scale systems, as they do not require to assemble the system matrix. In addition, they can solve systems inexactly to gain time. For instance, the linear solver may be used with a large tolerance in the first Newton iterations and then the tolerance becomes tighter along iterations to enforce superlinear convergence of the Newton method (Nocedal and Wright, 2006, Chapter 11, Theorem 11.3). On the other hand, iterative solvers are sensitive to the system matrix conditioning and often need preconditioning to perform efficiently. Direct solvers, which rely on a factorization of the system matrix, are less sensitive to bad conditioning and can achieve higher accuracy in system solutions. Though matrix assembly and factorization represent a significant overhead, modern parallel implementations of direct solvers are very efficient. For these reasons, direct solvers are often preferred in our application domain, where linear systems keep reasonable sizes in order to keep up with the real-time requirement. The interested reader may refer to the book by [START_REF] Saad | Iterative methods for sparse linear systems[END_REF] and the review by [START_REF] Davis | A survey of direct methods for sparse linear systems[END_REF] for more information about iterative and direct solvers, respectively.

Algorithm 2 gives a summary of a standard Newton method. Newton methods are popular due to their fast convergence rate, provided that F is regular enough and K has good properties in a vicinity of the solution. The following theorem (Nocedal and Wright, 2006, Chapter 11, Theorem 11.2) states the theoretical convergence rate of the method.

Theorem 5.1. Suppose that F is continuously differentiable, let u ⋆ a solution of problem (5.11) such that K(u ⋆ ) is nonsingular, and let (u k ) the sequence of iterates generated by Algorithm 2. Then when u k is close enough to u ⋆ , the Newton method converges at a superlinear rate, i.e.

u k+1 -u ⋆ = o u k -u ⋆ . In addition, if F is Lipschitz continuously differentiable, then convergence is quadratic, i.e. u k+1 -u ⋆ = O u k -u ⋆ 2 .
Theorem 5.1 only states a local convergence result. In other words, the initial iterate should be close enough to the solution for the method to converge toward this solution. As a consequence, the Newton method suffers from a lack of robustness when iterations start far away from a solution.

In addition, monitoring convergence along iterations is sometimes difficult. The stopping criterion in Algorithm 2 is the residual norm ∥r k ∥. Though the residual norm decreases dramatically when getting close to the solution, it is not supposed to decrease in general and can even take very large values before convergence. As a consequence it is difficult to evaluate the progress of the method during execution.

In optimization algorithms based on the Newton method, some features were added to better enforce and monitor convergence, at the expense of efficiency.

Newton method for optimization problems

We now consider the elasticity problem under the form (5.10). Treating a nonlinear elastic problem as a minimization problem is a way to obtain useful insight, especially by looking at the properties of the objective function. In particular, the evolution of the objective value along iterations is convenient to monitor convergence in the first iterations. It is complementary to the residual norm, which is more meaningful in the last iterations. By requiring the objective function to decrease along iterations, optimization algorithms enforce global convergence toward a local minimizer of the problem. In this paragraph, we describe a trust-region Newton method that solves optimization problems. Trust-region methods have already been considered in the context of nonlinear mechanics [START_REF] Gross | On the Convergence of Recursive Trust-Region Methods for Multiscale Nonlinear Optimization and Applications to Nonlinear Mechanics[END_REF][START_REF] Youett | A Globally Convergent Filter-Trust-Region Method for Large Deformation Contact Problems[END_REF]. While we mainly refer to [START_REF] Nocedal | Numerical optimization[END_REF], the interested reader might find more specific information about trustregion methods in the book by [START_REF] Conn | Trust-region methods[END_REF].

Though it is based on the standard Newton method, the trust-region Newton method uses a slightly different formalism. During an iteration, the objective function u → W (u) -b T u is approximated around a given iterate u k by its second-order Taylor expansion. The local quadratic model reads

m k (w) = W (u k ) -b T u k objective value + (F(u k ) -b) objective gradient T w + 1 2 w T K(u k )
objective Hessian w.

(5.12) Now, the Jacobian K plays the role of the Hessian of W , and its properties depend on the convexity of W . In particular, it is positive definite whenever W is strictly convex.

Here we assume that K(u) is nonsingular everywhere. A standard Newton step, starting from u k , lands on the critical point of m k . However, this critical point is a minimizer of m k only when K(u k ) is positive definite. Otherwise, it may be a saddle point or a maximum of m k , while a minimum does not necessarily exist.

For this reason, the Newton step is reformulated into an optimization subproblem which aims to minimize the local model on a ball centered at u k . The ball, called the trustregion, defines a zone where the quadratic model m k provides a suitable approximation of the objective function, and enforces the existence of a solution to the subproblem. The trust-region subproblem reads min w m k (w) subject to ∥w∥ < ∆ k (5.13)

where ∆ k > 0 is the trust-region radius. Provided that m k is strictly convex and ∆ k is large enough, solving (5.13) boils down to performing a standard Newton step. This situation arises in particular when iterates approach a solution of the optimization problem, and thus the trust-region Newton method benefits the convergence rates stated in Theorem 5.1.

In practice, the trust-region problem is approximately solved using a truncated Conjugate Gradient method (Nocedal and Wright, 2006, Algorithm 7.2), which stops upon trust-region bound violation or generation of a direction of nonpositive curvature. In a recent paper, [START_REF] Dahito | The Conjugate Residual Method in Linesearch and Trust-Region Methods[END_REF] compare the performance of the Conjugate Gradient method and the Conjugate Residual method for trust-region problems.

Algorithm 3: Trust-region mechanism (Nocedal and Wright, 2006, Alg. 4.1).

Data:

Initial radius ∆ 0 , reject tolerance η ∈ [0, 1/4) for k = 0, 1, 2, . . . do Compute a step w by solving (5.13) Evaluate ρ k // Adjust trust-region radius if ρ k < 1/4 then Reduce trust-region: ∆ k+1 < ∆ k else if ρ k < 1/4 and ∥w∥ = ∆ k then Expand trust-region: ∆ k+1 > ∆ k else ∆ k+1 = ∆ k end // Accept or reject step if ρ k > η then Accept step: u k+1 = u k + w else Reject step: u k+1 = u k end end
The trust-region mechanism is summarized in Algorithm 3. To accept or reject a step w as well as to adjust the trust-region radius, the mechanism relies on the coefficient chosen optimization software should support sparse matrices and limited-memory quasi-Newton Hessian approximations.

ρ k = f (u k ) -f (u k + w) m k (0) -m k (w) ,
In this section, we use the notation g = ∇Φ(b) to denote the objective gradient. It has nothing to do with the surface force distribution.

Unconstrained problem

We first illustrate the adjoint method on a registration problem whithout L ∞ constraint. Dropping the L ∞ constraint makes sense in a finite dimensional context, as this constraint is not indispensable for solutions to exist. In addition, displacement regularity can be enforced by cheaper means such as penalty terms, which saves the cost of constraint management. As far as the optimization solver is concerned, the first-order optimality conditions for problem (5.4) in the unconstrained case read

∇Φ(b) = 0.

First-order methods

In the early FreeFem++ implementation, we implemented two simple first-order methods to validate the approach. We first used a standard gradient descent with backtracking linesearch [START_REF] Armijo | Minimization of functions having Lipschitz continuous first partial derivatives[END_REF]. Then we implemented an accelerated gradient descent [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], where some inertia is given to the sequence of iterates to avoid zigzagging phenomena (see also [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF].

In Figure 5.6, we show an example of registration involving a two-dimensional mesh. A point cloud is generated from a deformed configuration and used to perform the registration. As expected, the Nesterov accelerated gradient descent decreases the objective function faster than the standard gradient descent, especially in the first iterations (Figure 5.6c). However, note (Figure 5.6b) that the surface force distribution reconstructed by the Nesterov method is more irregular than the one reconstructed by the gradient descent. Indeed, accelerated and higher-order methods converge faster because it is easier for them to reach irregular surface forces distributions. On the other hand, the gradient descent method is very sensitive to the regularizing effet of the adjoint method, but it is also the reason why its convergence is slow.

An off-the-shelf quasi-Newton solver

In the Python implementation, we use a standard quasi-Newton solver available in the Scipy package, namely the L-BFGS-B algorithm [START_REF] Byrd | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF]. This limitedmemory BFGS method is normally tailored for bound-constrained problems, but remains relevant for our case. We create a toy problem by applying a surface loading onto a truncated sphere mesh of radius 1 (10,385 vertices) with a Dirichlet boundary condition applied on its base. Then, we sample the deformed surface to create a point cloud of 200 points. We perform the registration on a coarser mesh (907 vertices), resulting in a direct problem of size 2,721. More details about the toy problem are given in Section 6.1. For both data generation and reconstruction, a linear elastic model is used, with E = 1 and ν = 0.49. We use the adjoint method to estimate nodal forces on boundary nodes that are not subject to a Dirichlet boundary condition, resulting in an optimization problem with 418 × 3 = 1254 unknowns. We set the penalty term to zero, and we use a limited-memory BFGS method available in the Scipy package, with at most 200 stored updates. Figure 5.7a illustrates the data generation process and Figure 5.7b shows the reconstructed deformation along with the nodal forces on the coarser mesh. The reconstructed deformation is correctly fitting the pointcloud. Figure 5.7c shows the convergence statistics. After 184 iterations (4.5 s on our configuration), the objective gradient norm has decreased by a factor 10 5 . Meanwhile, the objective value, which is supposed to converge toward zero, has also decreased by 10 5 , which means that the average distance d(y, S u ) has decreased by a factor 10 2.5 ≈ 316.

Problem with pointwise constraint

In the discrete context, the L ∞ constraint takes the form of a pointwise (or should we say nodewise) constraint on the components b i of the control. Though the physical meaning is not exactly the same as the continuous formulation, the spirit is similar, namely avoiding too strong nodal forces. If we denote by I b the set of nodes where b i is nonzero, the constraint reads

∀i ∈ I b 1 2 ∥b i ∥ 2 ⩽ 1 2 M 2 , ( 5.14)
and the first-order optimality conditions read

∀i ∈ I b            ∂ i Φ(b) + µ i b i = 0 µ i ⩾ 0 1 2 ∥b i ∥ 2 -1 2 M 2 ⩽ 0 1 2 µ i ∥b i ∥ 2 -M 2 = 0, (5.15)
where µ i ∈ R + is the Lagrange multiplier associated with the i-th pointwise constraint.

Projected gradient and fixed-point iteration

Like previously, we preferred to use simple first-order algorithms to solve the pointwise constrained problem in the FreeFem++ implementation. As the pointwise constraint applies separately on each component b i , it is very easy to perform projections onto the feasible set

B = (b i ) i∈I b | ∀i ∈ I b ∥b i ∥ ⩽ M
by projecting each b i onto a ball of radius M . Thus, the projection operation reads

∀i ∈ I b Proj(b) i = b i if ∥b i ∥ ⩽ M M b i ∥b i ∥ if ∥b i ∥ > M.
We exploited those cheap projections in a projected gradient method [START_REF] Bertsekas | On the Goldstein-Levitin-Polyak gradient projection method[END_REF]). An overview of the projected gradient iteration is given in Algorithm 4. At each iteration, the step length is selected using a projected backtracking linesearch, and optimality is measured using the projected gradient norm ∥ Proj(bg) -b∥, where g = ∇Φ(b).

Algorithm 4: A projected gradient descent Data: Initial guess b 0 , tolerance ε for k = 0, 1, 2, . . . do Evaluate objective gradient:

g k = ∇Φ(b k ) if ∥ Proj(b k -g k ) -b k ∥ < ε then return b k end Select a step size β k by doing a linesearch Compute the next iterate b k+1 = Proj(b k -β k g k ) end
We also implemented a variant of the projected gradient descent, inspired from fixedpoint methods. When they converge, fixed-point methods based on optimality conditions may yield fast and robust convergence. In our variant, we replace the orthogonal projection by a so-called fixed-point step, derived from the optimality conditions. When R(b) = 0, the first line in (5.15) can be reformulated

∀i ∈ I b p i + µ i b i = 0.
As a consequence, at optimality, b i should be colinear to p i if µ i > 0. To enforce this condition, the fixed-point step reads

FPStep(b) i = b i if ∥b i ∥ ⩽ M -M p i ∥p i ∥ if ∥b i ∥ > M.
The FPStep operator is used in the last line of Algorithm 4 instead of the orthogonal projection. Note that the fixed-point step has a stronger effect than the projected gradient step. If we only consider components that are affected by the step, the fixed point step is equivalent to applying a projected gradient step with β k = +∞ to these components. This is the reason why the fixed-point method has good convergence properties in good conditions but might also show some instability. Figure 5.8 shows the convergence results on a simple two-dimensional problem. In the case of the projected gradient, the projected gradient norm decreases by a factor 10 4 in about 600 iterations, and a visually correct matching is achieved. The fixed-point method, however, seems to converge in a few iterations, but not toward a solution of the optimization problem. In particular, the resulting deformed mesh does not match the point cloud. Other tests with the fixed-point method exhibit behaviors where the mesh oscillates between two deformations, none of which is a solution of the optimization problem.

After several tentatives, we decided to drop the fixed-point approach due to its instability, and to keep using standard optimization methods. Figure 5.8: Results for the 2D problem after 1000 iterations of the projected gradient method and the fixed-point method. Constraint violation is not plotted, as projected method only generate feasible iterates.

A primal-dual algorithm with BFGS approximation

While using projected directions is very convenient in the context of gradient descents, it is way more difficult when using methods involving a change of metrics, such as Newton and quasi-Newton methods (Bertsekas, 1982). As a consequence, standard projected-directions algorithms are limited to problems with bound constraints [START_REF] Lin | Newton's method for large bound-constrained optimization problems[END_REF][START_REF] Byrd | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF], while our L ∞ constraint is nonlinear. For this reason, we turned to primal-dual methods in the Python implementation. Though the functioning of such methods is out of the scope of this manuscript, let us mention that they manage constraints using Lagrange multipliers (Nocedal and Wright, 2006, Chapters 17-19). Primal-dual methods include interior-points methods such as Ipopt [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], sequential quadratic programming algorithms like SLSQP [START_REF] Kraft | A software package for sequential quadratic programming[END_REF], or augmented Lagrangian methods [START_REF] Hestenes | Multiplier and gradient methods[END_REF]. We encountered difficulties to find an adequate optimization method, as most primaldual programs are Newton methods that require the objective Hessian, or only provide full BFGS approximations, which is not relevant for large problems. Though, we still present an example involving the same sphere problem as before, where we add constraint (5.14) with M = 10 -2 (The maximal nodal force in the solution to the unconstrained problem was ∥b∥ ∞ = 3 • 10 -2 ). We solve the optimization problem using the sequential quadratic programming algorithm SLSQP by [START_REF] Kraft | A software package for sequential quadratic programming[END_REF], available in the Scipy library.

Figure 5.9a shows the deformed sphere after the constrained reconstruction, with the same color scale as in Figure 5.7b for nodal forces. Note that nodal forces do not exceed the limit M in intensity. The convergence statistics are plotted in Figure 5.9b. As the internal variables of the optimization solver are not available from the Python interface, we used the projected gradient norm to measure optimality conditions. Note that, due to the constraint, the objective function does not decrease toward zero as previously, but converges toward a positive value. In addition, the forces intensities graph, along with the plotted reconstruction, suggest that every component of the nodal force distribution tends to reach the maximal allowed intensity M . Note that, as soon as all constraints are active (iteration 50), the projected gradient norm seems to decrease with more difficulty.

These examples shows the feasibility of a L ∞ constraint on the control. However, the L ∞ constraint is a nonlinear constraint that requires a primal-dual method to solve the optimization problem. Such methods are more relevant when second-order information about the objective function is available and are difficult to use with quasi-Newton approximations. On the other hand, assembling a Hessian matrix for the adjoint method is out of question due to the composite nature of the objective function Φ. For this reason, we do not investigate algorithms to deal with the constrained problem any further.

Remark 5.4. The best choice to manage the L ∞ constraint is probably the interiorpoints strategy, as it only generates that satisfy the constraint. Solving the elasticity problem with large forces using the Newton method requires many iterations, and in addition the optimization solver might converge toward a minimizer devoid of physical meaning. The interior-points strategy circumvents this possibility by preventing the magnitude of nodal forces to exceed M .

Numerical results

We now present a few applications of our approach in the context of augmented surgery. As a patient's life is at stake during a surgical operation, an augmented reality system for operation rooms is expected to be robust and accurate. In addition, such a method should be able to update the augmented view in real-time. With the following numerical tests, we try to evaluate the performances of our approach.

Preliminary investigations featuring the toy problems

Before we consider cases from the application domain, we present a couple of examples involving the toy problems used to illustrate the optimization methods in Section 5.5.

2D toy problem

In the case of the FreeFem++ implementation, the test case involves a two-dimensional truncated circular mesh of radius 1, where the distance from the disc center to the flat face is 1/2. The flat face of the mesh is subject to a homogeneous Dirichlet condition, while forces are applied on the circular boundary to generate a displacement. A point cloud with 22 points is created using a mesh with 690 nodes. An elastic deformation is applied to the mesh and a point cloud is extracted from the deformed boundary. For data generation as well as for reconstruction, a linear elastic model is used with E = 1 and ν = 0.49.

3D toy problem

The toy problem for the Python implementation is a three-dimensional version of the previous one. We create a mesh with the shape of a truncated sphere of radius 1, where the distance from the sphere center to the truncating plane is 1/2. The flat surface of the mesh is subject to a homogeneous Dirichlet condition, while forces are applied on the round surface to generate a displacement. After an elastic deformation is applied to the mesh, the round surface is sampled to create point clouds with different sizes. The registration procedure involves a similar mesh. The flat surface is still subject to a homogeneous Dirichlet condition, while the round surface is deformed to match the point cloud. In addition, nodal forces are allowed on every node which is not subject to the Dirichlet condition (i.e. every node from the round surface).

We generate data using a mesh with 10,385 nodes and 6,702 faces. Linear elasticity is used with E = 1 and ν = 0.49. Unless stated otherwise, the same model and parameters are used in the examples below. Figure 5.7a shows the initial mesh and the generated deformation. The mesh used for registration is generally coarser.

No accuracy guarantee, even in an ideal case

Let us put it straight from the start: our registration method does not come with any guarantee in terms of displacement accuracy. Compared to other methods, an increased displacement accuracy might result from a more rigorous physical modeling, but no theoretical result confirms that.

To illustrate this disclaimer, we evaluate the displacement accuracy of the procedure in an ideal case. First, the point cloud is very dense (10,000 points), contains no noise and therefore provides a good representation of the deformed surface. In addition, the same mesh is used for data creation and for reconstruction, which evacuates possible discrepancies between two meshes representing the same shape.

The evolution of the objective function is plotted in Figure 6.1, along with the error statistics in terms of displacement field and force distribution. The objective function keeps decreasing during the execution, and evaluates to 2 • 10 -7 , meaning that the quadratic mean of the distance between data points and the mesh approximately amounts to 6 • 10 -4 . Though, the displacement error seems to remain stable after 50 iterations. The displacement error settles between 0.006 and 0.3 for half of the points, and the maximal error is 0.6. These figures are to be compared to the sphere radius, which is 1. The lower part of Figure 6.1 shows the evolution of the magnitude and error distribution among nodal forces. These results give us a few insights about the registration procedure. First, the reconstructed nodal forces do not exceed the maximal intensity of nodal forces used to create data (see also Figure 6.2a). As the objective function only measures the discrepancy between the point cloud and the reconstructed surface, the algorithm tends to choose among the simplest deformations that make them match, regardless of the original deformation. Here, a simple deformation is one that requires low forces intensities. The (absolute) error on nodal forces seems to increase at the same rate as forces intensities along iterations, and both distributions remain in the same order of magnitude. As a consequence, we understand that, in general, it might be difficult to give a quantitative meaning to the computed distribution, especially when it contains a large number of degrees of freedom.

In Figure 6.2, the true and reconstructed deformations are compared. Green shapes represent the true deformation, while yellow shapes represent the reconstructed deformation. The true and reconstructed nodal forces distributions are plotted in Figure 6.2a. Despite the large error between individual nodal forces, the same accumulation of forces on top of the mesh appears in both distributions. In Figure 6 the individual position of each node was not retrieved by the procedure. In Figure 6.2c, we applied the same pattern to the true and reconstructed surfaces to highlight the difference of position between both surfaces. Though both surfaces are consistent with the point cloud, points from the initial sphere were not transported at the same place by the true and reconstructed displacement fields. This behavior is expected, as nothing prevents the reconstructed surface from sliding along the surface represented by the point cloud. In Figure 6.2d, we plotted the displacement error in the reconstructed mesh. The displacement error increases when we move away from the fixed boundary, and it reaches 0.6 at the top of the object (the sphere radius is 1).

To increase the accuracy of the method, more information should be included in the model, such as landmark correspondence, curvature information or distinguishing free and loaded surfaces.

Regularizing effect of the adjoint method

An advantage of the optimal control approach is that the mesh boundary nodes where forces apply can be chosen in advance, according to the physical model. In contrast, standard registration methods, by adding artificial forces depending on the data, impose that the forces that create the displacement be located close to data points. Using the 2D toy problem, we compare the registration result from the adjoint method with the result from a displacement-based method, which is similar to standard registration methods.

In the displacement-based method we implement, the optimization problem reads min u∈R 3n J(u).

Here, the optimization variable is u. To make the objective function decrease, u should move in the direction opposite to the gradient of J. Following the extension-regularization framework (discussed in Section 4.1.2), instead of moving in the standard gradient direction -∇J(u), we compute a gradient using the linear elastic inner product. This results in the scaled iteration

u + = u -βA -1 ∇J(u),
where A is the stiffness matrix. In other words, the opposite gradient direction -∇J(u) is transformed into a force that generates a direction descent in the space of displacements, which is exactly what is done in standard registration methods. As ∇J(u) is entirely supported by boundary nodes, the successive displacement fields are ensured to be elastic displacements created by a surface force distribution.

Remark 6.1. If ∇J(u) had some components inside the mesh, for instance if there were some landmarks inside the mesh, the resulting force distribution would also have components inside the mesh. This is how little control on the force distribution support we have when using a displacement-based method.

The displacement fields resulting from the adjoint and displacement-based methods are plotted in Figure 6.3a, along with the surface force distribution necessary to create the displacement. For the displacement-based method, these forces are located only on vertices close to data points, which supported the gradient ∇J(u) along iterations. The deformation computed using the adjoint method is more regular and the force distribution is spread all over the mesh boundary. Let us be more specific about the impact of the adjoint method. Both deformations are solutions to the optimal control problem, as the feasible set allows (but does not require) force distributions to be spread all over the mesh boundary. But the forces in the right figure actually spreading all over the mesh boundary is a consequence of the regularizing effect of the adjoint method. Another optimal control method would not necessarily return such a regular distribution without an additional regularization mechanism. Figure 6.3b shows the convergence and displacement error statistics. As the toy problem was generated using a regular force distribution, it could be expected that the volume registration error is lower in average for the adjoint method. The objective function graph illustrates the degradation of numerical properties mentioned in Remark 5.1. By adding a regularization layer when computing descent directions, the adjoint method makes also convergence more difficult. As the problem has a geometrical meaning, the space of displacements (in other words the geometrical space) is more relevant to solve it, and for this reason convergence is way faster for the displacement-based method.

Irregularity due to noise

In the context of augmented surgery, the point cloud is the result of an image processing pipeline, and as a consequence it comes with a certain amount of noise. In this case, the deformed mesh from the registration process should not meet every point in the cloud, as an irregular point cloud would yield an artificially irregular surface on the augmented reality view. Such numerical artefacts might be misleading for the medical staff. It is preferable that the mesh boundary pass through the point cloud without matching each point individually.

In this section we consider two simple approaches to regularize the computed displacement field in the mesh, in the context of the FreeFem++ implementation (i.e. where the control is g). Those two approaches are those suggested by the optimal control formulation. In the penalized approach, we set R(g) = α 2 ∥g∥ 2 L 2 , and the L ∞ constraint is dropped. In the constrained approach, the penalty function R is set to 0 while the L ∞ bound is set to a positive value M .

Comparison on a partial matching case

We experiment those two solutions on the two-dimensional toy problem. We add gaussian noise with standard deviation σ = 0.02 to the point cloud, and we solve the registration problem using Ipopt with a full BFGS Hessian approximation. The displacement error statistics for several parameters are displayed in Figure 6.4a. The median error seems constant across box plots, but the maximal displacement error depends on the chosen strategy. The penalty strategy did not manage to decrease the displacement error. Starting from the same distribution as the unregularized case for small coefficients, the displacement error does not reach a minimum before the penalty is too strong for the mesh to match the point cloud. On the other hand, a significant decrease is achieved with the constraint strategy, with a minimum for M = 0.4. Note that the final force distribution in the noise-free case has a maximal intensity of 0.5.

In Figure 6.4b, we plotted the resulting force distribution for the unregularized case, the constrained case that achieves the minimal displacement error (M = 0.4), and the first penalized case where the error distribution is impacted (α = 10 -1 ). The constraint and penalty strategies both result in more regular surface forces, but the constrained result still features an irregular surface. However, the case M = 0.4 did a satisfying job. The corresponding displacement errors are plotted in Figure 6.4c. It is noteworthy that the main source of error is not the noise in Γ, but rather the point cloud not covering the whole mesh boundary, resulting in large displacement error far from the points. This phenomenon appears in the unregularized and the penalty case. Finally, the lower error in the constrained case is just an uncontrolled consequence of the mesh being pushed inwards to compensate for the intensity bound. On the other hand, the penalty approach just resulted in a more regular surface without impacting the remaining of the mesh, which is no more than what we expect from a regularization strategy. Far from breaking the tie between the two strategies, this example shows that using a regularization has an impact on the region close to the point cloud as well as on regions far from the point cloud.

Comparison on a full matching case

For the sake of curiosity, we also compare the two regularization strategies on a full surface-matching case. Now, the point cloud represents the whole deformed surface, so that the error is less likely to come from the limited point cloud coverage. The error statistics are plotted in Figure 6.5, along with the force distribution and displacement error field for a few selected cases. It seems that constrained approaches managed to reduce the displacement error where it is the largest (close to certain points). This time, the penalty strategy was effective to decrease the displacement error. In particular, when α = 10 -3 , the maximal displacement error is three times smaller than in the unregularized case.

The conclusion of this section is certainly not that one of those two strategies is better than the other in terms of noise regularization. The choice of a strategy should take account of other arguments in favor of one or the other option. Those arguments include ease of implementation for the penalty and the possibility to reject excessive forces for the constraint. Remark 6.2. In the early implementation, we also experimented the H 1 penalty, defined by R(g) = α 2 ∥∇g∥ 2 L 2 . However, when the penalty coefficient α was sufficiently large to regularize the mesh boundary, it also resulted in a quasi-uniform force distribution. In addition, this penalty was difficult to use in the three-dimensional cases, where the optimization variables are the nodal forces b. For these reasons, we quickly dropped this option.

Sparse Data Challenge dataset

The accuracy of a registration method is critical for the augmented surgery application. A good accuracy allows the surgeon to target small structures inside the liver while avoiding blood vessels. In practice, clinical collaborators of the Mimesis team usually require an accuracy of 5 millimeters. Though we did not state any theoretical result about accuracy, here we evaluate this criterion experimentally.

We use the adjoint method on a liver registration problem involving the Sparse Data Challenge1 dataset. The dataset consists of one tetrahedral mesh representing a liver phantom in its initial configuration and 112 point clouds acquired from deformed configurations of the same phantom [START_REF] Collins | Improving Registration Robustness for Image-Guided Liver Surgery in a Novel Human-to-Phantom Data Framework[END_REF][START_REF] Brewer | The image-to-physical liver registration sparse data challenge[END_REF]. To generate data, the challenge organizers used a silicon liver phantom. They created deformations by laying the phantom on irregular supports on its posterior face, and acquired pictures of the anterior face with a camera to produce the point clouds (Rucker et al., 2014, Figure 2). The phantom also contains 159 targets whose position is measured in reference and deformed configurations to establish ground truth data. The position of targets remains unknown to us, and the registration error is computed by the challenge website after we upload the final mesh position associated with each point cloud. By keeping the ground truth hidden from participants, the challenge organizers ensure that there is no bias associated with knowledge of the ground truth.

Due to the acquisition process, the forces that create the deformation are only applied on the posterior surface of the liver, while the surface represented by the point cloud is the anterior surface. As a consequence, we label the posterior surface as the loaded surface, where the force distribution is allowed to take nonzero values, while the anterior surface is labeled as the matching surface, which is used to evaluate the objective function. Figure 6.6 illustrates the definition of the matching and loaded surfaces.

We first perform a rigid alignment between the matching surface and the point cloud using the standard Iterative Closest Point method [START_REF] Besl | Method for registration of 3-D shapes[END_REF]. The provided point clouds are split into four parts which correspond to parts of the liver that are easily identifiable, namely, the right and left ridges, the falciform region and the general surface. This information is used to initialize the rigid registration process. The point cloud and the mesh after rigid registration are shown in Figure 6.6.

We now turn to the elastic registration process. To ensure uniqueness of the solution to the direct elastic problem, we set a fixed boundary condition on a small zone of the posterior face. As the liver mesh represents a phantom, no information about blood vessels is available and for this reason we just choose six adjacent triangles close to the center of the posterior face to apply the fixed displacement constraint. We solve problem (5.3) (with R(b) = 0) using the adjoint method. As only small deformations are involved in this dataset, we use a linear elastic model for the liver phantom, with E = 1 and ν = 0.4. The procedure is stopped after 200 iterations (80s on our configuration).

Let us add a word about the choice of elasticity parameters, which are not specified by the challenge organizers. Actually, a small part of the ground truth is available to participants in order to calibrate some parameters. This sample consists of the position of 35 targets (out of 159) after deformation for 4 cases (out of 112). The partial ground truth was used to determine the Poisson ratio ν. Concerning the Young modulus E, it has no impact in this study, and for this reason we set it to 1. Indeed, in this example we retrieve a displacement, and we are not interested in the force intensities. Even for nonlinear models, E is a multiplicative scaling parameter in the elastic energy expression, which mean that the force distribution necessary to generate a given displacement field is proportional to E. As a consequence, changing the Young modulus just means changing the intensity of estimated forces, while the computed displacement remains the same. Figure 6.7 shows the target registration errors associated with the partial ground truth data provided by the challenge authors (position of 35 targets for 4 point clouds). The box plot shows that significant improvement is achieved by the elastic registration step with respect to the rigid registration result. In average, the elastic registration improves the target registration error by 46%. However, these error results are very partial, while the only error statistics for the full dataset are those displayed by the challenge website after submission of the results. In Table 6.1, we reported the target registration error statistics for all datasets. We compare these results with other submissions to the challenge in Figure 6.8. We obtained the second-best result displayed on the challenge website. In particular, our performance is very close to that of the Vanderbilt team (the challenge organizers), who also use approaches based on a mechanical model [START_REF] Heiselman | Intraoperative Correction of Liver Deformation Using Sparse Surface and Vascular Features via Linearized Iterative Boundary Reconstruction[END_REF]. This confirms that approaches that care about respecting the direct model are relevant in augmented surgery.

Local force estimation

Is it possible to give a meaning to the force distribution reconstructed by the optimal control method? If it applies onto a large part of the liver boundary, probably not. However, when the distribution is very local, estimating the resulting net force may be considered. In this section we consider an example associated to robotic surgery. In remote laparoscopic surgery, the surgeon manipulates the patient's liver by controlling a robotic arm. In this context, estimating the force applied by the robotic instrument is necessary, as an excessive force applied onto the organ may cause serious damage to living tissues and threaten the intervention outcome. On the other hand, it is difficult to embed force sensors on instruments that evolve in an environment like the human body, mostly due to geometry, biocompatibility and asepsis constraints [START_REF] Marbán | Haptic Feedback in Surgical Robotics: Still a Challenge[END_REF]. For this reason, certain standard surgical systems, such as the famous daVinci robot are not equipped with force sensors. Indirect methods based on image processing have been proposed to keep track of the applied force without a sensing system inside the patient's body [START_REF] Nazari | Image-Based Force Estimation in Medical Applications: A Review[END_REF].

Here we estimate a force in a context similar to that in [START_REF] Haouchine | Vision-Based Force Feedback Estimation for Robot-Assisted Surgery Using Instrument-Constrained Biomechanical Three-Dimensional Maps[END_REF], where the authors monitor the evolution of the force applied by an instrument using a sequence of point clouds extracted from laparoscopic frames. We generate synthetic test cases using a liver mesh of 3,046 vertices and a linear elastic model (E = 20, 000 Pa, ν = 0.45). Dirichlet conditions are applied in the region where the hepatic vein enters the liver and in the falciform region. To simulate an instrument manipulating the liver, we apply a very local force onto the mesh. We then sample the deformed mesh around the force application point to create a point cloud. Figure 6.9 shows the displacement generated by synthetic local forces, along with the reconstructed deformation using the sampled point cloud of 500 points.

Single estimation

We begin by estimating the single local force from the case illustrated in Figure 6.9. To do so, we allow the force distribution b to be nonzero only in a small zone (about 50 vertices) surrounding the triangles concerned by the traction force. We solve the optimization problem with a relative tolerance of 5 • 10 -4 on the objective gradient norm, and we compute the net force estimation as the sum f est = i b i ∈ R 3 of all the nodal forces of the reconstructed distribution b. The estimated force is compared to the ground truth f true by computing the relative error ∥f est -f true ∥/∥f true ∥. Note that working with nodal forces makes it easy to work with several meshes, as we do not have to take the area of surface triangles into account.

In figure Figure 6.9, the same mesh is used to generate the data and to perform the registration. As a consequence, performing an accurate registration is easy, as the Figure 6.9: Left: synthetic deformation generated by a local force. The black wireframe represents the reference configuration. Right: reconstructed deformation using the point cloud. Nodal forces are summed to produce a net force estimation.

shape of tetrahedrons appears in the point cloud. The reconstructed force distribution and the estimation error are represented in Figure 6.10, on the left. At the end of the optimization process, the relative error on the net force is below 1%. Actually, we noticed that, with a denser point cloud (5,000 points), the reconstructed force seems to converge toward the ground truth force. This bias appears when the same discretization is used for data generation and reconstruction. On the right part of Figure 6.10, we show the result of registration performed with a different mesh (3,829 vertices), obtained by remeshing the original mesh. The force distribution is noisier and nodal forces take larger values, as tetrahedrons must be stretched to match the point cloud. Though, the noise is somehow filtered by summing all nodal forces, leaving a net force estimation with an error below 10%.

In Figure 6.11, we present the same effect with a refined mesh, (19,234 nodes) obtained by splitting the original mesh. The L ∞ norm of nodal forces, plotted in the top right graph, increases along iterations. The bottom images show the force distribution at iterations 57 (when the estimation error is minimal) and 293 (last iteration). Though the force distribution at the last iteration is very noisy and takes large values, the estimation error remains stable around 8%.

Estimation of a variable force

To mimic the action of a robotic tool manipulating the liver in the field of view of a laparoscopic camera, we create sequences of point clouds by applying a time-dependent force to a small spot on the liver surface. For each frame in the sequence, the resulting displacement is computed from the current force and a part of the deformed boundary is sampled to create a point cloud of 500 points. Each sequence consists of 50 forces, with a displacement of about 1 mm in a random direction between two successive forces. We create 5 sequences, where the force is applied on different points.

The registration is performed using the remeshed mesh of 3,829 vertices, already used in Figure 6.10 on the right. For each sequence, the observed point clouds are successively fed to the procedure to update the force estimation by solving the optimization problem.

The optimization algorithm is initialized with the last reconstructed distribution, so that only a few iterations are required for the update. Table 6.2 shows the average error obtained for each sequence, together with the average execution time of the updates and the number of evaluations of the objective function. The noise in the reconstructed distribution results in an intensity smaller than the reference, which is the main cause of an overall error of 10 % in average. In all cases, the update time is below 1.5 s in average, which is an encouraging result to show the feasibility of our approach in a real-time context. Note that this test case is only a toy problem, which is far from representing clinical conditions. In our result, the main source of error is the use of two different meshes for data generation and registration. However, other sources of error arise in a clinical context, as the physical model only provides a rough approximation of the real system. For instance, the estimated net force is proportional to the Young modulus used in the liver model. Patient-specific parameters may be estimated using elastography, but nowadays this procedure is not part of a standard clinical routine. According to [START_REF] Oudry | Comparison of four different techniques to evaluate the elastic properties of phantom in elastography: is there a gold standard?[END_REF], an error of 20 % is typical for clinical elastographic measurements. The elastographic error, combined with other sources of error, might amount to an overall error that prohibits the estimation of the robot force. Though, our results are at least encouraging in terms of performance. In a clinical software, it is not expected to use meshes much finer than those used in this study, in particular because of memory requirements due to computing on graphic processors.

Toward force estimation with a hyperelastic model

All along the project, we encountered difficulties to use the adjoint method with a nonlinear hyperelastic model, because of the instability of standard Newton methods available in most finite element packages. This point is discussed in Section 5.4. However, we are now in a position to present an early convergence result involving a Neo-Hookean model. Elastic parameters are the same as previously.

Like previously, we created synthetic deformations by applying a local force onto the liver mesh. As hyperelastic models are less rigid than the linear elastic model, applying nodal forces onto one or two triangles only results in a very irregular displacement field.

To produce more realistic deformations, we spread the synthetic force on a slightly larger zone. The reconstructed force distribution is also spread on a zone larger than before.

In this example, the finite element model is managed using the Dolfin package [START_REF] Logg | DOLFIN: Automated Finite Element Computing[END_REF][START_REF] Logg | DOLFIN: a C++/Python finite element library[END_REF]. To solve the direct problem, we combine both Newton methods described in Section 5.4. For the first iterations, we use a trust-region Newton method available in the Scipy package. This Newton method is guaranteed to converge, but it is also very slow. As a conjugate gradient method is used to solve linear subproblems, the trust-region Newton method can only solve the nonlinear problem up to a residual norm of 10 -9 . To obtain a solution with a tighter tolerance, we perform a few iterations with the standard Newton method available in Dolfin to reach a residual norm of 10 -14 . Figure 6.12 shows convergence statistics and resulting deformations for two test cases. Here the same mesh is used for data generation and reconstruction. In both cases, we performed as many adjoint iterations as possible, letting the solver fail when the direct solver accuracy is not sufficient to compute descent directions. Note that this failure did not happen after the same number of iterations for both examples. For case 1, iterations stopped after 10 min 49 s, while for case 2 they stopped after 1 h 4 min 11 s (on a single thread). However, in both cases the adjoint method managed to decrease the gradient norm by approximately four orders of magnitude. This result shows the feasibility of the adjoint method with a nonlinear elastic model. From the point of view of the inverse problem, however, the very noisy force distributions illustrate the difficulty associated to an enlarged support for nodal forces.

Future developments and conclusion

This project results from discussions between Yannick Privat and Stéphane Cotin, about the possibility to use a mathematical point of view to tackle the liver registration problem. The early work was quite prospective, and the optimal control approach we came up with is the result of numerous discussions, negociations and also false leads. While some choices result from tests and comparison between several methods, other choices were made by default in order to move forward in the project. The first category includes the choice of tree-based procedures to perform orthogonal projections, while the second category includes the choice to use an adjoint method or to set Dirichlet boundary conditions. In this chapter, we discuss a few possible enhancements for our work. We did not implement those options by lack of time, but they could improve computational performances or remove difficulties from the problem.

Avoid unnecessary difficulties by dropping Dirichlet boundary conditions

All along this manuscript, we only considered static elasticity problems involving a Dirichlet boundary condition. The main role of the Dirichlet boundary condition is to ensure the existence of a solution displacement u for a given force distribution g. Its role is also, to a lesser extent, to take into account anatomic factors that limit the liver motion. In liver surgery, the Dirichlet condition is often set at the hepatic vein entry, or sometimes at places where ligaments hold the liver (e.g. the falciform region). Note that the Dirichlet boundary is necessary because a static elastic model is used. When a time-dependent dynamic model is used, existence of an equilibrium configuration is not required. However, using a Dirichlet condition is far from being the best option, and it was even a source of difficulty in certain parts of the project. From a theoretical point of view, the presence of mixed boundary conditions degrades the regularity of solutions to the elasticity system. In particular, it is common to see stress accumulation at the interface between Dirichlet and Neumann boundary conditions. Due to this lack of regularity, we failed to prove the existence of a solution to the optimization problem in the general case. From a numerical point of view, setting a Dirichlet boundary condition before elastic registration sets the initial rigid registration result in stone. If the rigid registration went wrong, this cannot be corrected during elastic registration. Finally, from the modeling point of view, a hard Dirichlet condition does not even model accurately liver/ligaments. Ligaments, as well as main blood vessels, limit the organ motion but do not prevent it completely. For all these reasons, getting rid of the Dirichlet boundary condition should remove a difficulty from the problem. Though, an alternative mechanism should be used to ensure the uniqueness of solutions to the direct problem.

Replacing Dirichlet condition with stiffnesses

A first possibility is to relax Dirichlet boundary conditions into Robin boundary conditions. Instead of being fixed, the corresponding surface ∂Ω D is subject to a force that attracts it toward a target position u D . This model seems particulary relevant when ligaments are modeled as linear or nonlinear springs. For instance, [START_REF] Nikolaev | Estimation of boundary conditions for patientspecific liver simulation during augmented surgery[END_REF] model ligaments that hold the liver as a neo-Hookean material, whose stiffness is estimated using a Kalman filter. To keep things simple, assume that interactions with ligaments are represented using a cubic stiffness. The elasticity problem then reads

min u W (u) + κ 3 ˆ∂Ω D ∥u -u D ∥ 3 ds - ˆ∂Ω N g • u ds, (7.1)
where κ > 0 is the spring stiffness. With the same notation as (2.2), the resulting strong formulation reads

     -div(w ′ (∇u)) = 0 in Ω 0 w ′ (∇u)n = g on ∂Ω N w ′ (∇u)n = κ∥u D -u∥(u D -u) on ∂Ω D . (7.2)
Using stiffnesses instead of Dirichlet boundary conditions may result in a more accurate modeling of the liver mechanics as well as in more flexibility to adjust rigid transformation in the shape matching problem. In addition, provided that forces due to ligaments are not concentrated at one point, the resulting displacement is solution to a Neumann problem with w ′ (∇u)n ∈ L ∞ (∂Ω 0 ), which might ensure better regularity than mixed boundary conditions.

Note that rigid displacements should be penalized enough by the stiffnesses to avoid unwanted large rotations of the liver during registration. In particular, the stiffness coefficients κ should be large enough, and the restrained surfaces ∂Ω D should be spread all around the surface so that the liver can resist against torques. Last but not least, this strategy requires to know the elastic parameters of ligaments holding the liver.

A pure traction elastic problem

Another option consists in using only Neumann boundary conditions and work in a subset of H 1 (Ω 0 , R 3 ) to ensure uniqueness of solutions to the direct problem. In this situation, all forces that apply on the liver boundary are represented by the unknown Neumann condition. Constraints and penalties on the force distribution g may be specialized, depending on types of forces that are expected on each part of the liver boundary. In particular, larger intensities may be allowed on regions that were previously subject to Dirichlet conditions, as thoses forces are responsible for balancing the forces that drive the registration.

Let us have a look at the linear elasticity case. The former Dirichlet region is now included in the Neumann region, and we denote by ∂Ω F the region where a free boundary condition must be applied. It is known that the null space of the bilinear form a(u, v) =

ˆΩ0

Aε(u) : ε(v) dx is the space of uniform translations T = {1 Ω 0 x, x ∈ R 3 }. Thus, the space H 1 (Ω 0 , R 3 ) can be decomposed into the direct sum T ⊕ T ⊥ (either for the L 2 or H 1 inner product), where

T ⊥ = v ∈ H 1 (Ω 0 , R 3 ) ˆΩ0 v dx = 0 .
Provided that the surface force distribution belongs to the zero-mean subspace

S ⊥ = {1 ∂Ω N x, x ∈ R 3 } ⊥ = g ∈ L 2 (∂Ω D )
ˆ∂Ω N g dx = 0 , then solutions in H 1 (Ω 0 , R 3 ) to the system

    
-div(Aε(u)) = 0 in Ω 0 Aε(u)n = g on ∂Ω N Aε(u)n = 0 on ∂Ω F .

(7.3) are equal up to a uniform translation. In particular, there is a unique solution u g that belongs to T ⊥ . Note that in this case, u is the least-norm representant of solutions to (7.3). Finally, to be consistent with the numerical development, we say a word about the case g / ∈ S ⊥ , where (7. Though the study of the continuous problem needs to be consolidated, the full-Neumann configuration can be translated in our numerical framework without increasing too much the size of the direct system. We now derive a linear elasticity operator adapted to this new situation. When it does not include Dirichlet conditions, the stiffness matrix A does not have full rank. With the same notation as before, its kernel and range read

T = ker A = {(x, x, • • • , x), x ∈ R 3 } and T ⊥ = im A = v ∈ R 3n n i=1
v n = 0 .

To construct the direct problem matrix, we also introduce the matrices S and S b defined by

∀x ∈ R 3 Sx = (x, • • • , x), [S b x] i = x if i ∈ I b 0 otherwise; ∀v ∈ R 3n S T v = n i=1 v i , S T b v = i∈I b v i .
The matrix S b has the same effect as S, but its action only affects nodes corresponding to ∂Ω N and has no effect on nodes corresponding to ∂Ω F or nodes inside the domain.

The direct problem can be formulated in a compact way using the system Coupled with a good shape-matching functional, those two solutions should be good candidates to replace the Dirichlet condition. Compared to the stiffness-restricted problem, the restriction-free problem seems more relevant when little information from the direct model is available, e.g. ligament stiffness is unknown. However, it can result in too many degrees of freedom, especially if the translational component of the displacement is also chosen as an optimization variable instead of being set to zero.

Jump to lightspeed with neural networks

After looking at update time results in Section 6.3 it is not clear that achieving real-time registration using the adjoint method is feasible, even with a better implementation of the adjoint procedure and a good trade-off between accuracy and computation cost. The adjoint method shares processor time with other parts of the pipeline, including point cloud generation from laparoscopic images and displaying the augmented view. This constatation holds when a linear model is used and the stiffness matrix is factorized once when the procedure begins. When a nonlinear model is used, however, solving the direct problem requires an iterative method, which means several linear systems to assemble and solve. In that case, at least, it seems clear that real-time computation is not feasible as is.

A neural network to predict elastic displacements

Like biomechanical simulation, deep learning is a central research interest in the Mimesis team, and for this reason we considered replacing the iterative nonlinear solver with a neural network. The network, developed in the team by [START_REF] Odot | DeepPhysics: A physics aware deep learning framework for real-time simulation[END_REF], is a fullyconnected multi-layer perceptron with four layers of size 3n. In such a structure, each neuron in a given layer is connected to every neuron from the previous layer and to every neuron from the next layer. Figure 7.1 shows the network structure. The result z i of the i-th layer is obtained from the previous layer using the formula

∀ 1 ⩽ i ⩽ 4 z i = f i (W i z i-1 + c i ), (7.6)
where parameters W i and c i result from the network training. The function f i applies pointwise an activation function to each component of its input. Here the chosen activation function is the Parametric Rectified Linear Unit, defined by

PRelu(x) = x if x ⩾ 0 ax if x ⩽ 0,
where a ⩾ 0 is a learned parameter specific to each neuron. The network input receives the nodal force distribution z 0 = b and the displacement field u = z 4 is read at the other end of the network. The network is trained by showing it a wide range of force distributions, along with the corresponding displacement fields. As neural networks mainly involve matrix-vector products and pointwise operations, implementations take advantage of the architecture of graphic processors to perform very fast predictions. A prediction operation, denoted u = N(b), is performed in approximately five milliseconds on a modern computer.

Optimization involving the network

The development of the adjoint method is carried out jointly with Alban Odot, while the network design and training is taken care of exclusively by Alban Odot. The structure

Input tensor

Output tensor

W 1 f 1 , c 1 W 2 f 2 , c 2 W 3 f 3 , c 3 W 4 f 4 , c 4
3n neurons of the adjoint method involving the neural network is shown in Figure 7.2. In our current implementation, the same optimization solver as before is used. However, network operations are run on the graphic card, using the PyTorch framework. The functional evaluation is performed on the GPU as well, by using the brute-force projection function from PyTorch3D mentioned in Section 5.3. As the optimization solver does not run on the GPU, this implementation requires a lot of exchange between the central processor and the GPU, which is a known source of overhead. In a future version, quasi-Newton operations should run on the GPU to minimize inter-device exchanges.

The whole backward chain is taken care of automatically by the automatic differentiation module in PyTorch. Yet, let us add a word about backward propagation in the neural network. This operation replaces the adjoint problem and involves the adjoint operator of the neural network. In Figure 7.1, information now circulates from the right to the left. The rightmost tensor receives the gradient s 4 = ∇J(u) in the space of displacements, while the adjoint state p = s 0 is read at the left extremity. The operation to progress a layer backward is the adjoint operation to (7.6), namely

∀ 1 ⩽ i ⩽ 4 s i-1 = W T i ∇f i (W i z i-1 + c i )s i ,
where ∇f i (W i z i-1 + c i ) is a diagonal matrix. In other words, the information circulates in the same wires and using the same weights as forward propagation. Each neuron, instead of applying the activation function, multiplies its input by the activation function derivative that has been saved during forward propagation. We obtained early results using a parallelepipedic beam mesh with a Dirichlet condition at one extremity. A synthetic deformation is created by applying a force distribution on the beam surface and its surface is sampled to create a point cloud. A reconstructed deformation is obtained in approximately 300 ms on our configuration. Figure 7.3 shows the synthetic and reconstructed beam deformation, along with convergence statistics. Beside the fact that the reconstructed beam is twisted the wrong way, the pointcloud matching seems visually acceptable. Despite the use of a L 2 penalty, high-frequency oscillations appear in the reconstructed deformation but keep a reasonable amplitude. Note that the Dirichlet condition is not enforced but learned by the network, and Dirichlet nodes may be subject to a small displacement in certain cases.

The beam result is encouraging, but our network-based optimization approach is far from being operational. In the next step of our investigations, we intend to perform nonlinear registration and net force estimation using a liver mesh.

Improve robustness with primal-dual PDE management

The whole numerical implementation presented in this manuscript revolves around an adjoint method, where the only variable seen by the optimization solver is the control b. To converge up to a sufficiently tight tolerance, the optimization solver needs the adjoint method to return an objective value and gradient with high precision. In other words, the whole adjoint pipeline, including the direct solver, the functional evalutation procedure and the adjoint solver (see Figure 5.1) should return a result with an error as close as possible to machine precision. For this reason, the adjoint method is relevant when the elasticity problem is solved using a direct procedure, such as a direct linear solver for the linear case or the neural network in the nonlinear case. However, when the elastic problem is solved using an iterative method, expecting a result with a closeto-machine precision at each iteration is very costly and unrealistic at the same time. In particular, if a variation of the nodal forces is too small, the solver returns the previous displacement field as a good enough solution and the whole process fails because the objective function does not decrease. In this section, we discuss an alternative to the adjoint method where elasticity constraint is managed using Lagrange multipliers. We first mention two primal-dual algorithms, the Augmented Lagrangian method and the Sequential Quadratic Programming method, that are good candidates for our application. As both methods are often used in the context of second-order methods, we also discuss the feasibility of evaluating a numerical Hessian for the shape functional J.

Primal-dual solvers

In a primal-dual algorithm, u and b are both optimization variables in (5.3), while the adjoint state p is used as a Lagrange multiplier to manage the constraint F(u) = b. The min-max problem associated to (5. During an iteration of the Augmented Lagrangian method [START_REF] Hestenes | Multiplier and gradient methods[END_REF], the primal and dual variables are updated alternatively. First, a minimization subproblem involving a variant of the Lagrangian is solved with a fixed Lagrange multiplier p. The minimization subproblem reads Compared to the standard Lagrangian, the augmented Lagrangian involves an additional penalty term to better enforce the constraint. In the second part of the iteration, the multiplier is updated using

p + = p -µ(F(u) -b),
so that, along iterations, the sequence of solutions to (7.7) converges toward the solution to (5.3). A simple summary of the Augmented Lagrangian iteration is given in Algorithm 5. The reader interested in more details about the tuning of parameters may consult Nocedal and Wright (2006, Chapter 17). The Augmented Lagrangian method is flexible, as the user can choose how to solve the minimization subproblem. In particular, constraints on u or b can either be handled by the inner subproblem solver, or managed by the augmented Lagrangian method by introducing additional Lagrange multipliers.

Algorithm 5: Augmented Lagrangian strategy. (Nocedal and Wright, 2006, Alg.17 Another method for equality-constrained problems is the Sequential quadratic programming method (Nocedal and Wright, 2006, Chapter 18). The main idea of this algorithm is to apply the Newton method to solve the first-order optimality conditions where ∇ 2 F(u) T p = ∂ ∂u (∇F(u) T p). Of course, sequential quadratic programming algorithms feature mechanisms to ensure convergence, such as a line search, trust-region or filter procedure. Though we do not give more details about sequential quadratic programming methods, let us add that sequential quadratic programming methods are usually expected to perform better on optimization problems with nonlinear constraints.

As they involve both u and b as optimization variables, primal-dual methods often lead to large problems and require to solve large linear systems. However, in our case, the objective function is splitted between a function of u and a function of b, leading to sparse systems. This splitting may even be exploited by using methods such as the Alternating direction method of multipliers (see Ouyang et al., 2015, and references therein). Finally, primal dual methods might be interesting if computations are sped up using reduced bases. In that case, an iterative procedure would still be necessary to perform the optimization, but the optimization problem would keep a reasonable size.

Computing a Hessian for J

In primal-dual methods, the displacement field u and the functional J are manipulated directly by the optimization solver instead of being hidden in the depths of the adjoint procedure. As a consequence, it makes sense to consider second-order information for J. Sequential Quadratic Programming methods require to solve a quadratic subproblem involving the Hessian matrix ∇ 2 J(u). Here, using a BFGS approximation is tricky, as it would deteriorate the sparsity structure of the system matrix. Concerning the Augmented Lagrangian method, solving the minimization subproblem with an increasing accuracy using a first-order or quasi-newton method would result in very long iterations. In both case, being able to evaluate a numerical Hessian for J seems profitable.

The strategy to compute numerically the Hessian of J is the same as the one to evaluate the gradient of J. For a given point y ∈ Γ, the projection point p ∂Ω u (y) lies either on a vertex, inside an edge or inside a triangle of the deformed mesh. In addition, in the majority of cases, when a small perturbation v is applied to u, the projection point p ∂Ω u+v (y) is expected to remain on the same element. As a consequence, the elementary function j y is locally equal to the squared distance between y and a point, a line or a plane. Let us take an example based on Figure 5.2. On the left, j y is locally equal to u → 1 2 ∥y -x k ∥ 2 , and the only nonzero term in the Hessian matrix ∇ 2 j y (u) is the diagonal block ∂ 2 kk j y (u) = I corresponding to x k . On the right, however, j y is locally equal to The expression for θ is obtained by minimizing ∥v + θu∥, while the expression for f (u, v) is obtained by noticing that v + θu and θu are orthogonal and applying the Pythagorean theorem (see notations in Figure 7.4). Here, θ depends on u and v and its derivatives read

∂θ ∂v (u, v) = - u ∥u∥ 2 ∂θ ∂u (u, v) = - v + 2θu ∥u∥ 2 .
It can be verified that the first-order derivatives of f , From this expression, it is difficult to obtain visual insight about j y , but it shows the feasibility of computing a Hessian matrix. In the case of a triangle element, computing second-order derivatives is much more involved, as our change of variable is not applicable anymore. However, a numerical implementation should take advantage of automatic differentiation to compute Hessian coefficients.

It might seem strange to compute the second-order derivatives of a function which is not Gateaux differentiable everywhere. However, as the optimization method converges, displacement steps get smaller and pairings between the mesh boundary elements and observed points are not expected to change too much. On the other hand, second-order information can be helpful in the first iterations of the optimization solver, to avoid taking too large steps.

Conclusion

In this manuscript, we presented a registration procedure based on an optimal control formulation of the registration problem. We not only proposed a structure for the implementation of a registration method, but also some mathematical insight that might be of use to understand the behavior of the discretized procedure. Of course, our approach is far from being mature, and many other choices could have been made, such as a more clever functional J, a more accurate physical model or more efficient numerical methods. However, new mathematically sound approaches in augmented surgery might be derived by adapting tools from the generic optimal control framework. Proposing open and easily adaptable procedures should help democratize computer-assisted surgery, letting more people around the world benefit from safer care and go back home from the hospital in good health.

Figure 2 :

 2 Figure 2 : Géométrie du problème. La distance entre un point y ∈ Γ et sa projection orthogonale sur ∂Ω u s'écrit d(y, ∂Ω u ).

  p.p. x ∈ ∂Ω N λ(x) = 0 si ∥g(x)∥ < M λ(x) ⩾ 0 si ∥g(x)∥ = M, tel que g vérifie la condition d'optimalité de premier ordre p.p. x ∈ ∂Ω N p(x) + (α + λ(x))g(x) = 0. Chapitre 5. Méthode d'adjoint pour résoudre le problème de recalage Nous discrétisons le problème par la méthode des éléments finis. Le problème de contrôle optimal discrétisé = J(u b ) + R(b), où b est le vecteur des forces nodales s'appliquant sur les sommets du maillage, et u b est le champ de déplacement solution du système élastique discrétisé F(u b ) = b. Dans la formulation discrète, la seule variable sur laquelle le solveur d'optimisation agit est le contrôle b. Le calcul des dérivées de J(u b ) par rapport à b nécessite l'utilisation d'un état adjoint. Plus précisément, le gradient de la fonction objectif s'écrit ∇Φ(b) = p + ∇R(b), où l'état adjoint p est la solution du problème adjoint F ′ (u b ) T p = ∇J(u b ).

  fonction objectif ∇Φ(b) = p + ∇R(b) Résolution du problème adjoint F ′ (u) T p = ∇J(u)

Figure 1

 1 Figure 1.2: Problem geometry. The distance between a point y ∈ Γ and its orthogonal projection onto Ω u is denoted by d(y, ∂Ω u ).

Figure 1 . 3 :

 13 Figure 1.3: Structure of the adjoint procedure. The forward chain (top) aims to evaluate the functional for a given control b, while the backward chain (bottom) transforms the gradient ∇J(u) into a gradient in the space of controls.

Figure 1 .

 1 Figure 1.4: Left: synthetic deformation generated by a local force. Right: reconstructed deformation using the point cloud.

Figure 2

 2 Figure 2.2: Left: liver geometry extracted from pre-operative tomographic images. Right: point cloud representing a part of the deformed liver surface.

Figure 2 . 3 :

 23 Figure 2.3: Resulting deformations for the uniform force distribution, using the linear model (bottom left), the Saint Venant-Kirchhoff model (top left) and the Neo-Hookean model (right).

  Figure 2.5: Problem geometry. The distance between a point y ∈ Γ and its orthogonal projection onto ∂Ω u is denoted by d(y, ∂Ω u ).

Definition 3. 2 (

 2 Outer normal). Let Ω a domain with C 1 boundary, x 0 ∈ ∂Ω and ψ the diffeomorphism from Definition 3.1. Without loss of generality we assume ψ(x 0 ) = 0 R d , and we denote by (e 1 , • • • , e d ) the canonical base of R d . The outer unit normal to ∂Ω at x 0 is defined by

  Definition 3.8 (Signed distance function). Let Ω ⊂ R d be an open set with Lipschitz boundary. The signed distance function with respect to Ω is defined by

  Figure 3.3: Set of equation y < x 4/3 . The signed distance function is not C 1 at (0, 0).

Figure 3

 3 Figure 3.4: The elementary sets B + , B + ∪ D and B + ∪ D 0 . A full line indicates that the boundary is included in the set. A set is Gröger-regular when its boundary locally looks like one of these sets, up to a Lipschitz diffeomorphism.

  .4) where P y (u) is defined in (4.2). In addition, denote by Γ ⊂ Γ the set where P y (u) is not a singleton. For y ∈ Γ \ Γ, denote by x y the single element of P y (u). If Γ has zero relative Lebesgue measure in Γ, dJ(u) is a linear continuous application on C(Ω 0 ) defined by⟨dJ(u), v⟩ = ˆΓ\ Γ v(x y ) • (x y + u(x y ) -y) dy. (4.5) Proof. Let B = B C(Ω 0 ) (0,r) a small ball of radius r > 0. Due to Lemma 4.1, for y ∈ Γ, j y is Lipschitz continuous on u + B. Therefore, for v ∈ B ∀t ∈ (0, 1)

Remark 4. 1 . 1

 11 Haller-Dintelmann et al. (2009, Remark 4.2) note that Proposition 3.11 can be stated with p in an interval of type [2 -δ, 2 + δ]. As a consequence, if Proposition 3.11 yields a δ > 0 such that 2 -δ < d/(d -1) (which is true in dimension 2), then we have dJ(u) ∈ C(Ω 0 ) ′ ⊂ W -1,2-δ D (Ω 0 ), and resulting descent directions are in W

Theorem 4. 3 .

 3 Assume that d = 3, Ω 0 has a C 2 boundary and p = 2. Then problem (4.8) has at least one solution.

  .11), featuring the nodal force vector b in the role of the control, and the displacement field u as the state. A possible formulation reads min u∈R 3n ,b∈B J(u) + R(b) subject to constraint F(u) = b, (5.3) where B is the set of admissible controls. Though the elastic system F(u) = b is denoted as an optimization constraint, it is a particular constraint, as a given control b yields a unique state u b so that the pair (u b , b) is feasible. In an adjoint method, everything is expressed as a function of b, making b the only variable controlled by the optimization solver. The resulting optimization problem reads min b∈B Φ(b) where Φ(b) = J(u b ) + R(b).

Algorithm 1 :

 1 Computation of the objective gradient using an adjoint method. Data: Current iterate b Compute the displacement u by solving F(u) = b Evaluate J(u) and ∇J(u) Compute the adjoint state p by solving F ′ (u) T p = ∇J(u) Result: ∇Φ(b) = p + ∇R(b)

Figure 5 . 1 :

 51 Figure 5.1: Structure of the objective evaluation procedure. The forward chain (top) aims to evaluate the functional for a given control b, while the backward chain (bottom) transforms the gradient ∇J(u) into a gradient in the space of controls.

  # x is the vector of node positions , y is the point from Gamma . k1 , k2 , k3 , t1 , t2 , t3 = n e a r e s t _ n e i g h b o r _ p r o c e d u r e (x , y ) x . requires_grad = True # Save operations performed on x from now on . j = 0.5 * norm ( t1 * x [ k1 ] + t2 * x [ k2 ] + t3 * x [ k3 ] -y )**2 # Functional value j . backward () # Backpropagation nabla_j = x . grad # Functional gradient

Figure 5

 5 Figure 5.4: Bounding box hierarchy around a triangular surface mesh. Small boxes are the bounding boxes of single triangles, while the biggest box is the bounding box of the whole surface. Red boxes represent an intermediate level of boxes.

Figure 5 . 5 :

 55 Figure5.5: Resulting deformation and convergence statistics of the standard Newton method and the trust-region method (including rejected steps for the trust-region method). Missing points in the Neo-Hookean case correspond to infeasible points where the objective function evaluates to +∞.

  (a) Initial (left) and deformed (right) configurations with generated point cloud. (b) Resulting deformations for the standard (left) and accelerated (right) gradient descents.

Figure 5 . 6 :

 56 Figure 5.6: Results for the 2D unconstrained problem after 500 iterations of the standard and accelerated gradient methods.

  (a) Initial sphere mesh and deformed configuration with generated point cloud. (b) Reconstructed deformation and nodal forces that generated the optimal displacement.

Figure 5 . 7 :

 57 Figure 5.7: Reconstruction results for the 3D unconstrained problem.

  (a) Resulting deformations for the fixed-point (left) and projected gradient (right) algorithms.

  (a) Reconstructed deformation and nodal forces. The color scale is the same as in Figure5.7b. Convergence statistics for the quasi-Newton method.

Figure 5 . 9 :

 59 Figure 5.9: Reconstruction results for the 3D problem with pointwise constraint.

Figure 6 . 1 :

 61 Figure 6.1: Top: objective value and displacement error at mesh nodes. Bottom: Intensity and error statistics for nodal forces.

  (a) Comparison between true (left) and reconstructed (right) forces. (b) Point cloud and superposition of the true and reconstructed surface.

  (c) Discrepancy between the true (left) and reconstructed (right) surface. (d) Displacement error.

Figure 6 . 2 :

 62 Figure 6.2: Registration results in the ideal case.

  (a) Displacement-based (left) and adjoint-based (middle) registration results. Evolution of the objective function and node displacement errors.

Figure 6 . 3 :

 63 Figure 6.3: Comparison between the displacement-based approach and the adjoint method.

  Nodal registration error for several penalty coefficients and constraint bounds.

  Resulting force distributions and displacement error for the case without regularization (top), the constrained case with M = 0.5 (middle) and the penalized case with α = 10 -3 (bottom).

Figure 6 . 5 :

 65 Figure 6.5: Comparison of several kinds of regularization for the full matching case.

Figure 6 . 6 :

 66 Figure 6.6: Left: Initial liver mesh with matching surface on top and loaded surface at the bottom. Right: Point cloud and liver mesh after rigid registration. Point colors correspond to the left ridge (blue), the right ridge (green), the falciform region (magenta), and remaining points (yellow), respectively.

Figure 6 .

 6 Figure 6.10: Reconstructed force distribution and evolution of the force estimation error for the original mesh (left) and the remeshed mesh (right). Nodal forces are summed to produce a net force estimation.

Figure 6 .

 6 Figure 6.11: Top: Evolution of the force estimation error and the maximal nodal force for the refined mesh. Bottom: reconstructed force distribution at iterations 57 (left) and 293 (right).

Figure 6 .

 6 Figure 6.12: Convergence statistics and resulting deformations for two test cases with the nonlinear model. The true deformation and force distribution are plotted on the left, while the reconstructed deformation and force distributions are plotted on the right.

  g∥ 2 dx s.c -div(Aε(u)) = 0 in Ω 0 Aε(u)n = 0 on ∂Ω F . (7.4) Of course, as S ⊥ is the orthogonal of S = {1 ∂Ω N x, x ∈ R 3 } in L 2 (∂Ω D ), solutions to (7.4) satisfy Aε(u)n = g -1 ∂Ω N |∂Ω N | ˆ∂Ω N g dx.

  only has nonzero coefficients at nodes in I b . Let us check that this system does what is expected from it. The first line reads Au = b -S b λ. A necessary and sufficient condition for solutions to this line to exist is to set λ as the mean of (b i ) i∈I b , so that b -S b λ lies in the range of A. Here, by using S b instead of S, we correct b with a uniform vector field defined only on the nodes of I b and avoid to artificially add forces to other nodes. It also ensures that, the adjoint state p resulting from the adjoint system solution satisfies i∈I b p i = 0, which is consistent with the requirement b ∈ T ⊥ . Now, with the first line only, u is defined up to an element of ker A. By ensuring u ∈ T ⊥ , the second line S T u = 0 enforces the uniqueness of u, which is, again, the least-norm solution of the first line. Thus, dropping the Dirichlet boundary only adds 3 lines to the direct problem.

Figure 7 . 1 :

 71 Figure 7.1: Structure of the multi-layer perceptron. The W i , c i and f i represent the parameters associated to each step of the network.

Figure 7 . 3 :

 73 Figure 7.3: Results of the network-based optimization procedure.

  , b, p) where L(u, b, p) = J(u) + R(b) -p T (F(u) -b).

  ) + R(b) -p T (F(u) -b) + µ 2 ∥F(u) -b∥ 2 .(7.7)

  .3) Data: Initial multiplier p 0 , initial penalty coefficient µ 0 , constraint violation tolerance ηfor k = 0, 1, 2, . . . do Find a pair (u k , b k ) that minimizes L(u, b, p 0 ) + µ k 2 ∥F(u) -b∥ 2 if ∥F(u) -b∥ ⩽ η then return (u k , b k ) end Adjust Lagrange multiplier: p k+1 = p k -µ k (F(u k ) -b k ) Choose a new penalty coefficient µ k+1 ⩾ µ k end

  ) -F ′ (u) T p = 0 ∇R(b) + p = 0 b -F(u) = 0.Namely, the quadratic subproblem that is assembled at each iteration to compute a step (w, (u) -F ′′ (u) T p ∇ 2 R(b) w h subject to constraint F(u) + F ′ (u)w = b + h,

  1 2 d 2 (y, (x k x m )), and the Hessian matrix should contain 4 nonzero blocks, at the lines and columns corresponding to x k and x m .To compute those Hessian coefficients, we use the change of variablev = x k + x m 2 -y and u = x k -x m 2 ,and we define the functionf (u, v) = 1 2 d 2 (y, (x k x m )) = 1 2 ∥v∥ 2 -θ 2 ∥u∥ 2 , where θ = θ(u, v) = -⟨u, v⟩ ∥u∥ 2 .

Figure 7 1

 71 Figure 7.4: Notations for computing the Hessian of j y .

  [START_REF] Mestdagh | An Optimal Control Problem for Elastic Registration and Force Estimation in Augmented Surgery[END_REF]. « An Optimal Control Problem for Elastic Registration and Force Estimation in Augmented Surgery ». In : Medical Image Computing and Computer Assisted Intervention -MICCAI 2022. Sous la dir. de Linwei Wang et al. Cham : Springer Nature Switzerland, p. 74-83. isbn : 978-3-031-16449-1. doi : 10.1007/978-3-031-16449-1_8

  ).

	Sequence N. eval. Update Time Relative error
	Case 1	9.2	1.42 s	8.9 %
	Case 2	5.6	0.85 s	16.2 %
	Case 3	5.5	0.84 s	5.7 %
	Case 4	5.4	0.82 s	4.4 %
	Case 5	6.2	0.96 s	15.0 %
	Table 6.2: Average number of objective evaluations, execution time per update and
	relative error for each sequence.

Code available at https://github.com/gmestdagh/adjoint-elastic-registration.

(∂Ω N ) , and we do not assume that F is differentiable. As a consequence, we prove an optimality condition that should hold for every adjoint state p ℓ at this point. When F turns out to be differentiable at g, there is only one adjoint state p that should satisfy the first-order optimality condition below.

https://trimsh.org/

Challenge website: sparsedatachallenge.org.

Remerciements

Proposition 4.2. Let u a continuous displacement field, and v ∈ C(Ω 0 ) a direction. We denote by (x y ) y∈Γ a family such that x y ∈ P y (u) and dJ(u)(v) = ˆΓ v(x y ) • (x y + u(x y ) -y) dy.

If w ∈ C(Ω 0 ) is another direction associated with the same family (x y ) y∈Γ , then for α and β two nonnegative coefficients, dJ(u)(αv + βw) = ˆΓ(αv + βw)(x y ) • (x y + u(x y ) -y) dy.

In other words, each direction in the cone generated by v and w is associated with the same family (x y ) y∈Γ .

Proof. First, note that for y ∈ Γ, x y ∈ P y (u), and as a consequence ˆΓ min

and, as α, β ⩾ 0, ˆΓ(αv

Finally, we obtain ˆΓ min

and αv + βw is also associated with the family (x y ).

Instead of one single linear form dJ(u) in the differentiable case, non-differentiable points yield as many linear forms as there are families of type (x y ) y∈Γ . Each one of these linear forms ℓ satisfies

As a consequence, a descent direction for one linear form is also a descent direction for J. In the numerical implementation, the projection algorithm only return one family (x y ), that is used to generate a descent direction.

where f (u) = W (u) -b T u denotes the objective value. The coefficient ρ k compares the model reduction and the actual reduction allowed by the step w. In particular, when a step w does not decrease the objective function as much as predicted by the model, it is rejected and the trust-region radius is decreased. Though rejecting steps represent an additional computational cost, it ensures the convergence of the trust-region Newton method, as stated in the following theorem (Nocedal and Wright, 2006, Chapter 4, Theorem 4.1).

Theorem 5.2. Assume η > 0 in Algorithm 3. Assume also that ∥K∥ is uniformly bounded, that the objective function W (u)-b T u is bounded below and that F is Lipschitz continuous. Let (u k ) be the iterates generated by the trust-region Newton method. Then

In Figure 5.5, we applied random nodal forces to the 2D beam from Figure 2.3, resulting in a very irregular displacement field. We plotted the convergence statistics of the standard and trust-region Newton method, in the case of a Saint Venant-Kirchhoff model and in the case of a Neo-Hookean model. Due to the good properties of the Saint Venant-Kirchhoff energy functional, the standard Newton method behaves well and converges much faster than the trust-region method. On the other hand, the trustregion method is the only one to converge in the more difficult case of the Neo-Hookean model. While the standard Newton methods fails as soon as an iterates lands onto a non-feasible point, the trust-region mechanism allows the trust-region method to take shorter steps when necessary.

The very aggressive standard Newton method and the more careful trust-region Newton method are only two examples among a full spectrum of variants with various levels of convergence enforcement. Strategies we did not mention include line search algorithms (Nocedal and Wright, 2006, Chapter 3), which are relevant for convex problems and quasi-Newton methods. The energy-minimization nature of static elasticity problems gives access to a large choice of strategies to enforce convergence. Compromises between all those strategies should be found to combine efficiency requirements with the robustness necessary in sensitive applications such as augmented surgery.

Optimization procedure

In this section, we discuss the choice of an appropriate optimization algorithm to solve the optimal control problem in the unconstrained case and in the pointwise-constrained case. Note that in this manuscript we do not solve the state-constrained problem (2.14) numerically. The chosen algorithms should satisfy two main requirements. First, secondorder information is not available for the objective function Φ, thus Newton methods are excluded from possible choices. For this reason, we rather turn toward first-order methods and quasi-Newton algorithms. Second, we are dealing with optimization problems that can be large, and thus storing full matrices is often impossible. For this reason, 44,57,67,84). As these datasets were used to calibrate the method, results are better for these sets than in average.

20-28%

28-36% 36-44% All Data Sets Surface Coverage