Stéphane Cotin
email: stephane.cotin@inria.fr

Guillaume Mestdagh
email: guillaume.mestdagh@inria.fr

Yannick Privat
email: yannick.privat@unistra.fr

Organ registration from partial surface data in augmented surgery from an optimal control perspective

Keywords: Shape registration, Augmented Surgery, Optimal Control

We address the problem of organ registration in augmented surgery, where the deformation of the patient's organ is reconstructed in real-time from a partial observation of its surface. Physicsbased registration methods rely on adding artificial forces to drive the registration, which may result in implausible displacement fields. In this paper, we look at this inverse problem through the lens of optimal control, in an attempt to reconstruct a physically consistent surface load. The resulting optimization problem features an elastic model, a least-squares data attachment term based on orthogonal projections, and an admissible set of surface loads defined prior to reconstruction in the mechanical model. After a discussion about the existence of solutions, we analyze the necessary optimality conditions and use them to derive a suitable optimization algorithm. We implement an adjoint method and we test our approach on multiple examples, including the so-called Sparse Data Challenge. We obtain very promising results, that illustrate the feasibility of our approach with linear and nonlinear models.

Introduction

In this paper, we consider an elastic solid, known in its rest configuration, but subject to an unknown deformation, that we need to estimate. The deformation is created by forces applied to the solid boundary. In addition to material elastic properties, the information available to reconstruct the displacement field inside the solid is an observation of a part of the deformed surface. In particular, correspondence between this observed surface and the initial surface is not known.

This inverse problem arises in computer-assisted surgery, where the elastic solid is, in fact, the patient's liver [START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF]. To help the medical staff navigate the operation theatre, an augmented reality system is used to visualize the motion of the liver internal structures (blood vessels, tumours, etc.), which requires to reconstruct the liver displacement. Tomographic images, taken before the surgical intervention, provide a detailed view of the liver in rest configuration, along with its internal structures. A three-dimensional pre-operative mesh is created from these images. On the other hand, during a minimally-invasive intervention, the liver is not directly visible, and all the available intra-operative information is a video feed from a camera inserted in the patient's abdomen. A point cloud, representing the visible part of the liver surface, is extracted from the current frame of the video and is used as an observation. Though this is not our goal in this paper, the liver displacement should be updated in real-time, as each frame from the camera brings a new observation.

Several methods for anatomical registration focus on the surface-matching problem [START_REF] Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF][START_REF] Charon | Metamorphoses of functional shapes in Sobolev spaces[END_REF][START_REF] Bauer | Chapter 16 -Metric registration of curves and surfaces using optimal control[END_REF][START_REF] Antonsanti | Partial Matching in the Space of Varifolds[END_REF][START_REF] Charon | Shape Spaces: From Geometry to Biological Plausibility[END_REF]. They aim to compute a diffeomorphism between the initial liver boundary and the observed surface, and their performance is measured in terms of surface similarity between the source and target surfaces. Pure surface-matching may be relevant for our problem when the observation represents the whole deformed boundary, or at least a large part [START_REF] Antonsanti | Partial Matching in the Space of Varifolds[END_REF], so that a volume displacement field can be reasonably computed from the reconstructed surface correspondence. In the augmented surgery setting, though, only one side of the liver appears in the camera field of view, and the observed point cloud typically covers less than 50% of the deformed boundary. Even with a perfect surface registration, this is too small to determine a liver volume displacement. To compensate for the lack of observed data, most intra-operative registration methods involve an elastic model of the liver. Liver mechanics has been much studied, with attention paid to the stress-strain law [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF], additional stiffness due to blood vessels [START_REF] Peterlík | Modeling and Real-Time Simulation of a Vascularized Liver Tissue[END_REF][START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF], and boundary conditions [START_REF] Ozkan | Compliance boundary conditions for patient-specific deformation simulation using the finite element method[END_REF][START_REF] Nikolaev | Estimation of boundary conditions for patient-specific liver simulation during augmented surgery[END_REF].

A whole category of existing methods for intra-operative registration inherits from so-called elastic image registration methods [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]Section II.A.1], where a cost function is minimized to enforce landmark or surface correspondence, while an elastic model is used to regularize the displacement field [START_REF] Burger | A hyperelastic regularization energy for image registration[END_REF][START_REF] De Buhan | An optimization method for elastic shape matching[END_REF][START_REF] Bauer | A numerical framework for elastic surface matching, comparison, and interpolation[END_REF][START_REF] Kurtek | Elastic Shape Analysis of Surfaces and Images[END_REF][START_REF] Hartman | Elastic Shape Analysis of Surfaces with Second-Order Sobolev Metrics: A Comprehensive Numerical Framework[END_REF]. These methods simulate fictitious forces that attract the liver toward the observed surface, in the same fashion as the Iterative Closest Point algorithm [START_REF] Besl | Method for registration of 3-D shapes[END_REF]. Artificial forces applied to the organ surface include electrostatic forces [START_REF] Suwelack | Physics-based shape matching for intraoperative image guidance[END_REF], linear and nonlinear springs [START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF][START_REF] Plantefève | Patient-Specific Biomechanical Modeling for Guidance During Minimally-Invasive Hepatic Surgery[END_REF], or Lagrange multipliers that enforce a so-called sliding constraint between the deformed surface and the observed point cloud [START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF][START_REF] Courtecuisse | Threedimensional physics-based registration of pelvic system using 2D dynamic magnetic resonance imaging slices[END_REF][START_REF] Morin | Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation[END_REF].

Unfortunately, such forces, that apply where the surface-matching occurs, do not reflect the real causes of displacement (namely, contact forces with other tissues). Whether they derive from an imaginary potential or act as constraint-enforcing Lagrange multipliers, they are created by the point cloud, which does not really exist. This results in unrealistic displacement fields, regardless of the elastic model elaborateness [START_REF] Rucker | A Mechanics-Based Nonrigid Registration Method for Liver Surgery Using Sparse Intraoperative Data[END_REF]. A more reasonable physical consideration is that the liver surface represented by the point cloud is visible, and therefore not subject to any mechanical contact force. In other words, using an elastic model does not guarantee that the resulting deformation is physically plausible.

Several methods prefer to solve an inverse problem, where a range of admissible deformations is defined based on physical hypotheses. Using linear elasticity, the authors in [START_REF] Heiselman | Characterization and correction of intraoperative soft tissue deformation in image-guided laparoscopic liver surgery[END_REF][START_REF] Heiselman | Intraoperative Correction of Liver Deformation Using Sparse Surface and Vascular Features via Linearized Iterative Boundary Reconstruction[END_REF] precompute a basis of displacements corresponding to forces applied in zones where ligaments hold the liver, while they impose a free boundary condition on the remaining surface. In [START_REF] Özgür | Preoperative liver registration for augmented monocular laparoscopy using backward-forward biomechanical simulation[END_REF], the initial liver pose is estimated by taking into account the effects of gravity and gas insufflation in the patient's abdomen. In [START_REF] Xie | Constrained finite element method for runtime modeling of soft tissue deformation[END_REF], forces that create the displacement are restricted to certain nodes of the liver mesh. Approaches based on an inverse problem often exhibit an improved registration accuracy and produce more physically relevant displacement fields, at the expense of execution time.

In this paper, we look at the registration inverse problem from the perspective of optimal control. While most registration methods are very specific to a given registration scenario, expressed in terms of finite-dimensional matrix-vector operations, and, often tailored for linear elasticity, we take a step backward and consider the continuous problem in a more general setting. Using the optimal control framework results in a very flexible formulation, yielding at the same time physically relevant displacement fields, and where it is easy to include new physical hypotheses or take additional intra-operative data into account. The wide range of tools provided by the optimal control framework may be of help to obtain precious mathematical insight about the registration problem. Last but not least, we expect our approach to lead to new numerical methods inspired from generic numerical optimization algorithms.

The optimal control problem studied in this problem was already introduced in our conference paper [START_REF] Mestdagh | An Optimal Control Problem for Elastic Registration and Force Estimation in Augmented Surgery[END_REF], which includes numerical examples associated to the augmented surgery domain. Outside the augmented surgery domain, other approaches based on the physical causes of displacement involve an optimal control problem. In [START_REF] Kaltenmark | Estimation of a growth development with partial diffeomorphic mappings[END_REF], the authors control the accretion process leading to the shape of a horn. In [START_REF] Hsieh | Mechanistic modeling of longitudinal shape changes: equations of motion and inverse problems[END_REF], the forces driving the motion of a viscoplastic material are estimated. The authors in [START_REF] Hogea | An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects[END_REF] estimate the parameters of a tumor growth model using a series of tomographic images of the same patient. The same parameter estimation problem is treated together with an atlas-to-patient registration problem in [START_REF] Gooya | GLISTR: Glioma Image Segmentation and Registration[END_REF][START_REF] Scheufele | Coupling brain-tumor biophysical models and diffeomorphic image registration[END_REF]. All these works revolve around a time-dependent direct model, while our direct model is an elliptic equation. In our case, the optimal control problem must be solved again and again each time the observed surface is updated. The authors in [START_REF] Sundar | Biomechanically-Constrained 4D Estimation of Myocardial Motion[END_REF] retrieve the muscle fiber tension in a biomechanical heart model to reconstruct the heart motion along a series of MRI images.

A static version of the same model is used in [START_REF] Ortigosa | Optimal Control of Soft Materials Using a Hausdorff Distance Functional[END_REF]. The authors control a soft robot to reach a target shape, that is fully known.

In the remaining of this paper, we begin by stating the physical model, along with the optimization problem we consider (Section 2). Then, we propose a mathematical analysis of the optimal control problem, including the existence of solutions and optimality conditions (Section 3). Finally, we show some numerical tests to evaluate the performances of our approach (Section 4). In particular, in the last test case, we extend our numerical approach to a nonlinear (Neo-Hookean) elastic model.

Towards an optimal control formulation

In our approach, we determine a surface force distribution on the organ boundary that generates a displacement compatible with observed data. The reconstructed force distribution is chosen among a set of admissible force distributions, which reflects a range of physical causes. Of course, we do not expect the registration procedure to accurately estimate the forces that created the observed deformation. The formulation we propose nevertheless allows us to get close to the physics of the problem in a satisfactory way, while leading to efficient algorithmic methods that are easy to implement. We illustrate our approach with physical examples.

Pre-operative biomechanical model

Figure 1 shows a sketch of the system we consider. In its reference configuration, the organ is represented by an open domain Ω 0 ⊂ R 3 with Lipschitz boundary, filled with an elastic material. Due to its interactions with its environment, the organ is subject to surface loadings, such as pressure or contact with surrounding tissues. Volume forces, such as gravity, are not considered in this study. The boundary ∂Ω 0 falls into two parts, ∂Ω D and ∂Ω N . On ∂Ω D , a homogeneous Dirichlet boundary condition applies, whereas the Neumann boundary ∂Ω N is subject to a surface force distribution g ∈ L 2 (∂Ω N , R 3 ). A homogeneous Dirichlet boundary condition is typically set in zones where main blood vessels enter the organ. The space of displacement fields associated with this partition reads

H 1 D (Ω 0 ) = H 1 D (Ω 0 , R 3 ) = u ∈ H 1 (Ω 0 , R 3 ) u = 0 on ∂Ω D ,
and we denote by u g ∈ H 1 D (Ω 0 ) the elastic displacement field generated by a given surface force distribution g on ∂Ω N .

Elastic stress-strain laws used in the literature to describe the behaviour of living tissues range from the linear elastic model [START_REF] Heiselman | Characterization and correction of intraoperative soft tissue deformation in image-guided laparoscopic liver surgery[END_REF] to hyperelastic models such as Neo-Hookean [START_REF] Miller | Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation[END_REF] or Ogden [START_REF] Nikolaev | Estimation of boundary conditions for patient-specific liver simulation during augmented surgery[END_REF] models (see also [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF] and references therein). When computation time is a constraint, the linear co-rotational model is sometimes preferred [START_REF] Nesme | Efficient, Physically Plausible Finite Elements[END_REF][START_REF] Plantefève | Patient-Specific Biomechanical Modeling for Guidance During Minimally-Invasive Hepatic Surgery[END_REF], as it generates nonlinear deformations at the cost of a linear system inversion. In this paper, we only consider the linear elastic system, except for a brief nonlinear example in the results section. In the linear case, the generated displacement u g solves the partial differential

equation (PDE)    -div(Aε(u)) = 0 in Ω 0 u = 0 on ∂Ω D Aε(u) • n = g on ∂Ω N , (1)
where the linearized strain tensor and the Hooke tensor are defined by We denote by Ω u = (Id +u)(Ω 0 ) the volume occupied by the deformed organ. Note that we do not assume (Id +u) to be a bijection in our theoretical study. This is consistent with our numerical framework, where, given the range of deformations expected in the surgery application, we do not need a mechanism to detect self-intersection of the deformed solid.

ε(u) =

Intra-operative data and surface-matching functional

While the reference configuration Ω 0 is known through pre-operative images, the current position shall be estimated from an intra-operative observation. In our continuous model, the observed data take the form of a compact two-dimensional surface Γ ⊂ R 3 . We denote by S 0 ⊂ ∂Ω 0 the matching surface, i.e. a part of the initial organ boundary that is supposed to partially match with Γ as the registration ends. We assume that S 0 is a compact subset of ∂Ω 0 known by the user, and we will use the notation

S u = (Id +u)(S 0 ) ⊂ (Id +u)(∂Ω 0 )
to denote the image of S 0 under a displacement u. The surface S 0 is chosen before the operation starts, as the surgeon has a vague idea of which side of the organ should be visible on the camera. Note that S 0 only acts as a cue for the algorithm, and only a part of it is supposed to match with Γ. If no information is available, it is still possible to set S 0 = ∂Ω 0 .

To measure the discrepancy between a given displacement and the observed data, we introduce the least-squares function J : C(Ω 0 ) → R, defined by

J(u) = 1 2 Γ d 2 (y, S u ) dy, (2)
where d(y, S u ) 2 = min x∈Su x -y 2 denote the Euclidean distance between y ∈ Γ and S u . Therefore, J(u) vanishes whenever Γ ⊂ S u , up to a zero Lebesgue measure set. There is a whole zoo of cost functions for partial shape-matching (see the reviews [START_REF] Van Kaick | A Survey on Shape Correspondence[END_REF][START_REF] Sahillioğlu | Recent advances in shape correspondence[END_REF] and reference therein), sometimes involving frameworks such as varifolds [START_REF] Kaltenmark | A General Framework for Curve and Surface Comparison and Registration With Oriented Varifolds[END_REF][START_REF] Antonsanti | Partial Matching in the Space of Varifolds[END_REF] or functional maps [START_REF] Rodolá | Partial Functional Correspondence[END_REF]. Most of them fit into our optimal control formulation, and some of them are even differentiable. However, they sometimes require to compute normal vectors or process the point cloud before the registration begins. On the other hand, in the context of real-time simulation, it is more appropriate for the cost functional to be as computationally simple as possible. Also, the elastic model already enforces displacement regularity, and an elaborate surface discrepancy functional is not really needed here. For these reasons, we fall back to [START_REF] Arguillère | Shape deformation analysis from the optimal control viewpoint[END_REF], which is a continuous version of functions previously used in augmented surgery [START_REF] Rucker | A Mechanics-Based Nonrigid Registration Method for Liver Surgery Using Sparse Intraoperative Data[END_REF][START_REF] Peterlík | Fast elastic registration of soft tissues under large deformations[END_REF]. Although it is not Gateau-differentiable and requires a special treatment when computing optimality conditions, this function can be quickly evaluated when Γ is a point cloud, and therefore is more relevant in the context of augmented surgery. We consider the constrained optimal control problem

min g∈G M Φ(g) where Φ(g) = J(u g ), (3)
where g and u g play the role of the control and the state, respectively. The feasible set reads

G M = g ∈ L ∞ (∂Ω N , R 3 ) g L ∞ (∂Ω N ) < M .
Here, the definition of G M is based on the physical hypothesis that the magnitude of surfaces forces that are expected in the human body does not exceed the positive constant M > 0. Note that, in our numerical examples, the definition of admissible forces also includes information concerning the support of g, i.e. the zones where surface forces apply on the organ boundary.

Setting an upper bound on the pointwise magnitude of g is also useful to handle noise in the observed surface Γ. In the medical application, Γ is a point cloud provided by an image processing pipeline. It is likely to include a certain level of noise. In this context, an exact matching between Γ and the organ boundary would result in a very irregular displacement field u g and oscillations in the surface force distribution. A better outcome would consist in the organ boundary passing through the point cloud without meeting every point individually. Using the constraint prevents the control g from taking the large values it needs to reach every point in Γ.

Analysis of the optimal control problem

To begin with, an analysis of the optimization problem should help us obtain some mathematical insight. After we discuss the existence of solutions to the continuous problem, we derive first-order optimality conditions. Optimality conditions will be useful when it comes to implementing a numerical method to solve the problem in silico. We take care of that part in the end of the section.

For q > 0, we introduce the notation W 1,q D (Ω 0 , R 3 ) to denote the Sobolev space

W 1,q D (Ω 0 , R 3 ) = u ∈ W 1,q (Ω 0 , R 3 ) u = 0 on ∂Ω D ,
and q > 0 is the conjugate exponent to q, given by 1/q + 1/q = 1. Then, the topological dual of W 1,q D (Ω 0 , R 3 ) with respect to the pivot space L 2 (Ω 0 , R 3 ) is denoted by W -1,q

D

(Ω 0 , R 3 ).

Well-posedness

Existence of solutions for Problem (3) is critical to guarantee the stability properties of the discrete problem with respect to the mesh size, or which regularity can be be expected for numerical solutions.

The following results are based on the work by Gröger [START_REF] Gröger | A W 1,p -estimate for solutions to mixed boundary value problems for second order elliptic differential equations[END_REF] around the regularity of solutions to partial differential equations with mixed boundary conditions. The definition of a Gröger-regular set is based on the following subsets of R 3 (with the notation x = (x 1 , x 2 , x 3 )):

B = {x ∈ R 3 | x < 1}, B + = {x ∈ B | x 3 > 0}, D = {x ∈ B | x 3 = 0}, D 0 = {x ∈ D | x 1 < 0}. Definition 3.1.
Let Ω 0 be a bounded open subset of R 3 and let ∂Ω N be a relatively open part of its boundary. The set G = Ω 0 ∪ ∂Ω N is said to be regular in the sense of Gröger if, for each x ∈ ∂Ω 0 , there is a neighbourhood U of x and a Lipschitz diffeomorphism

Ψ : U → B such that Ψ(U ∩ G) either coincides with B + , or B + ∪ D, or B + ∪ D 0 . Remark 3.1. A simpler definition, specific to the three-dimensional case, is provided in [44, section 5]: G is Gröger-regular when Ω 0 is a bounded domain of class W 1,∞ , ∂Ω D = ∂Ω 0 \ ∂Ω N is closed, and the boundary ∂(∂Ω N ), seen as a subset of ∂Ω 0 , is W 1,∞ .
Let us mention [START_REF] Gröger | A W 1,p -estimate for solutions to mixed boundary value problems for second order elliptic differential equations[END_REF][START_REF] Gröger | Resolvent estimates in W -1,p for second order elliptic differential operators in case of mixed boundary conditions[END_REF][START_REF] Droniou | Solving convection-diffusion equations with mixed, Neumann and Fourier boundary conditions and measures as data, by a duality method[END_REF][START_REF] Haller-Dintelmann | Hölder continuity and optimal control for nonsmooth elliptic problems[END_REF], investigating the W 1,q or C 0,α regularity of solutions to linear and nonlinear elliptic PDEs with mixed boundary conditions. Unfortunately, these results do not apply to the linear elasticity system, for which the literature is sparser. In [START_REF] Shi | Higher integrability of the gradient in linear elasticity[END_REF], W 1,q regularity is obtained for solutions to the linear elastic system on a domain of class C 1 . In [START_REF] Herzog | Integrability of displacement and stresses in linear and nonlinear elasticity with mixed boundary conditions[END_REF], the authors extend Gröger's framework to a class of linear and nonlinear elastic systems that satisfy an ellipticity condition with respect to the linear strain tensor ε(u). We state their result in the specific case of the linear elasticity system. Proposition 3.1. [48, Theorem 1.1] Assume that Ω 0 ∪ ∂Ω N is Gröger-regular, and that the Lamé coefficients µ, λ satisfy µ > 0 and 2µ + 3λ > 0. Then there exists a q 0 > 2 such that, for all q ∈ [2, q 0 ], the solution to

(1) is in W 1,q D (Ω 0 , R 3 ) provided that g defines an element of W -1,q D (Ω 0 , R 3
). In addition, there exists C > 0, depending only on Ω 0 and q, such that

u W 1,q (Ω 0 ) C g W -1,q (Ω 0 ) .
Our existence result relies on the compact embedding W 1,q (Ω 0 ) → C(Ω 0 ) for q > dim R 3 [49]. In this regard, Proposition 3.1 is not fully satisfying, as it only guarantees q 0 > 2. In the context of elliptic systems, the articles [START_REF] Krumbiegel | Second Order Sufficient Optimality Conditions for Parabolic Optimal Control Problems with Pointwise State Constraints[END_REF]Appendix] and [START_REF] Disser | Optimal Sobolev Regularity for Linear Second-Order Divergence Elliptic Operators Occurring in Real-World Problems[END_REF] contain some sufficient conditions for the condition q 0 > 3 to hold, but it is not clear that their conditions are relevant in the case of elasticity. For this reason, we keep the condition q 0 > 3 as an assumption in our existence result, stated below. Note that, if we consider the same problem in R 2 , such an assumption is not necessary as q 0 > dim R 2 is already guaranteed by Proposition 3.1. Also, if one considers a simpler toy problem where the elastic system is replaced with an elliptic system, the existence result remains valid in any dimension without any assumption, using a C 0,α result from [START_REF] Droniou | Solving convection-diffusion equations with mixed, Neumann and Fourier boundary conditions and measures as data, by a duality method[END_REF]. Should this result extend to the linear elastic system, the existence of solution would be guaranteed in any dimension for our problem. We refer to [START_REF] Mestdagh | An optimal control formulation for organ registration in augmented surgery[END_REF]Chapter 4] for additional explanations. 2. In R 3 , assume that Ω 0 and ∂Ω N are such that q 0 > 3 with the notations of Proposition 3.1. Then Problem (3) has a solution.

Characterization of optimizers

We now turn to first-order conditions satisfied by local minimizers of Problem (3). Before we state the optimality conditions themselves, let us have a look at the properties of the discrepancy functional J, defined in (2). In the following proposition, we state the differentiability of J. This result is fully based on Danskin's theorem [START_REF] Danskin | The theory of max-min and its application to weapons allocation problems[END_REF]Chapter 3] concerning the differentiability of functions defined by a minimum. For a given displacement field u ∈ C(Ω 0 ), we use the notation

P y (u) = {x ∈ S 0 | x + u(x) -y = d(y, S u )} = (Id +u) -1 Π Su (y), (4)
where Π Su (y) is the set of orthogonal projections of y onto S u . In particular, P y (u) ⊂ S 0 .

Proposition 3.2 (Proof in Appendix A A.2). Assume that Ω 0 is bounded with Lipschitz boundary, and let u ∈ C(Ω 0 ). Then the functional J has directional derivatives at u. Its derivative in the direction

v ∈ C(Ω 0 ) reads dJ(u)(v) = Γ min x∈Py(u) [v(x) • (x + u(x) -y)] dy. (5)
In addition, if P y (u) is a singleton for almost every point y ∈ Γ, then dJ(u) is a continuous linear form on C(Ω 0 ), i.e. J is differentiable in the sense of Gateaux at u.

Though it is not Gateaux-differentiable a priori, the application J always has directional derivatives. Nondifferentiability (in the sense of Gateaux) occurs when P y (u) contains several elements for too many points in Γ. There are two reasons why P y (u) can contain several elements: either y has several orthogonal projections onto S u , or y has a single projection point p Su (y) which is the image of several points x 1 , x 2 , • • • ∈ S 0 under the transformation Id +u.

Note that the directional derivative of J also reads

∀v ∈ C(Ω 0 ), dJ(u)(v) = min ∈L(u) , v , (6)
where the set of linear forms L(u) is defined by

L(u) = : C(Ω 0 ) v → Γ v(x y ) • (x y + u(x y ) -y) dy x y ∈ P y (u) .
All the linear forms involved in the definition of L(u) belong to the concave subdifferential of J, as they satisfy dJ(u)(v) (v) for all v ∈ C(Ω 0 ). In addition, according to the expression above, for a given direction v, there exists a family (x 0,y ) y∈Γ , with x 0,y ∈ P y (u) for all y ∈ Γ, such that

dJ(u)(v) = Γ v(x 0,y ) • (x 0,y + u(x 0,y ) -y) dy.
We now state necessary first-order conditions that characterize a minimizer g of Problem (3). Due to the structure of the derivative dJ, the first-order conditions consist of an equality satisfied by each element ∈ L(u g ). In particular, the linear form ∈ L(u g ) associated to the family (x y ) y∈Γ is represented in the space of controls by an adjoint state p . From a formal point of view, we aim to introduce p as the solution to the adjoint problem

   div(Aε(p)) = 0 in Ω 0 p = 0 on ∂Ω D Aε(p) • n = Γ (x y + u(x y ) -y) δ xy dy on ∂Ω N , (7)
where δ xy is the two-dimensional (pointwise) Dirac measure at x y on ∂Ω N . To rigorously define a notion of solution for such a system, we use the so-called transposition method. To this aim, we assume, as in Theorem 3.1, that Ω 0 and ∂Ω N are such that q 0 > 3 with the notations of Proposition 3.1, and we consider a fixed q ∈ (3, q 0 ].

For a given h ∈ L q (∂Ω N ), we denote by w h the solution to the PDE

   div(Aε(w)) = 0 in Ω 0 w = 0 on ∂Ω D Aε(w) • n = h on ∂Ω N . (8)
In what follows, q denotes the conjugate exponent to q. Definition 3.2. One says that p ∈ L q (∂Ω N ) solves [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF] in the sense of transposition if

∂Ω N p • h ds = Γ w h (x y ) • (x y + u(x y ) -y) dy.
for all h ∈ L q (∂Ω N ), where w h solves [START_REF] Peterlík | Modeling and Real-Time Simulation of a Vascularized Liver Tissue[END_REF].

We refer to the monography [START_REF] Ponce | Elliptic PDEs, measures and capacities[END_REF] for additional explanations about solutions to PDEs with measure right-hand sides, as well as the seminal work [START_REF] Littman | Regular points for elliptic equations with discontinuous coefficients[END_REF]. Note that, with this definition, the adjoint state p is only defined on ∂Ω N . Though the values of the adjoint state inside Ω 0 are not needed for our analysis, p can easily be extended to the whole domain Ω 0 using [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF]. Lemma 3.1 (Proof in Appendix A A.3). Let us assume that Ω 0 ∪ ∂Ω N is Gröger-regular and let q ∈ (3, q 0 ]. Then, the adjoint problem (7) has a unique solution in the sense of transposition.

The first-order optimality conditions are stated below.

Theorem 3.2 (Proof in Appendix A A.4). Assume that Ω 0 ∪ ∂Ω N is Gröger-regular and M ∈ [0, ∞].
Let g ∈ G M a local minimizer of Problem (3) and u g the associated displacement field. If ∈ L(u), denote by p the associated adjoint state defined by [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF]. For every ∈ L(u), 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.

Remark 3.2 (Practical use of optimality conditions.). The three facts mentioned in this remark justify some implementation choices made in the numerical section that slightly differ from the theoretical framework. First, it is notable that if g ∈ G M is a local minimizer of Problem (3) such that g < M almost everywhere on ∂Ω N , then Φ is Gateaux-differentiable at g with dΦ(g)(•) = 0. Indeed, if the L ∞ constraint is inactive, then λ = 0 by the so-called slackness property. In this case, the optimality condition simply reads p = 0 on ∂Ω N . Using Equation [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF], we obtain = 0, and therefore L(u g ) = {v → 0}, which means that Φ is Gateaux-differentiable at g. Furthermore, outside from a local minimizer, a descent direction for the function J can be obtained using only one projection per point y ∈ Γ. Indeed, performing a single projection per point y is equivalent to finding one linear form 0 ∈ L(u). When u is not a local minimizer, 0 = 0, and if v is chosen so that 0 , v < 0, then, using [START_REF] Charon | Shape Spaces: From Geometry to Biological Plausibility[END_REF],

dJ(u)(v) = min ∈L(u) , v 0 , v < 0,
meaning that v is also a descent direction for the whole criterion J.

Finally, even when taking into account the L ∞ constraint in the numerical algorithm, we observe that the constraint is not saturated as soon as M is chosen large enough. This suggests (but this is only a numerical conjecture) that the problem still has a solution if we remove the L ∞ constraint in the optimal control problem definition.

Based on these remarks, we chose to ignore the L ∞ constraint in our numerical examples, and to work with J as if it were Gateaux differentiable everywhere.

Numerical methods

To solve the registration problem numerically, we discretize the initial domain Ω 0 using a tetrahedral mesh, while displacement fields on Ω 0 are represented using P1 finite element functions. Due to the complex geometry of organs, using a tetrahedral mesh is very common in augmented surgery. In addition, choosing piecewise linear functions is more convenient when it comes to computing orthogonal projections onto the deformed mesh. Concerning the observation Γ, it is provided in the form of a point cloud.

We adopt a 'discretize-then-optimize' approach, which means that the entire formulation is transformed into a finite-dimensional problem, which is then solved using numerical tools for finite-dimensional optimization. Indeed, an algorithm based on the 'optimize-then-discretize' might be very sensitive to discretization errors, and might encounter difficulties to find adequate directions of descent by applying the continuous approach on a discretized problem. We have therefore chosen to adapt mutatis mutandis the tools developed in the previous section to a discrete framework.

From now on, the discrete displacement field on Ω 0 is represented by the finite-dimensional vector u = (u 1 , . . . , u n ) ∈ R 3n where n is the number of vertices and u k ∈ R 3 is the displacement of the k-th vertex. The same bold letter u is used to denote the associated finite element function. We control the vector of nodal forces b = (b 1 , . . . , b n ) at the mesh vertices, and the resulting elastic displacement u b is determined by solving the linear system

Au = b,
where A is the stiffness matrix of the liver mesh. Nodal forces b are related to the surface force distribution g by the equality b = Sg, where the matrix S represents the L 2 inner product on ∂Ω N . Controlling directly the nodal forces b instead of the surface force distribution is a good way to lighten computations, at the expense of consistency with the continuous problem. It remains consistent with the discretize-then-optimize approach, though.

Discrete objective function.

Let us first have a look at the discretized functional J. In the discretized framework, the organ surface ∂Ω 0 is a triangular mesh, and defining S 0 consists in selecting the set of triangles that is used to evaluate J. The observation Γ is provided as the point cloud Γ = {y 1 , . . . , y p } ⊂ R 3 , and the discretized functional reads

J(u) = 1 p p i=1 j i (u) where j i (u) = 1 2 d 2 (y i , S u ). (9)
Our implementation uses the nearest-neighbour search from the Trimesh Python package1 to evaluate the distance between a point y ∈ Γ and S u . This projection procedure stores triangles from S u in a spatial indexing structure [START_REF] Guttman | R-Trees: A Dynamic Index Structure for Spatial Searching[END_REF] to perform efficient nearest-neighbour queries. Using Remark 3.2 and the fact that nondifferentiability only occurs on a negligible subset of the space of displacements, we only need the nearest-neighbour search procedure to return a single projection point a = p Su (y) per observed point y ∈ Γ. Note that handling an unknown number of projections per point y would represent an additional computational cost, and would require special attention to avoid computing every p . We have therefore chosen a compromise between efficiency and complexity of the algorithm. Since we are dealing with piecewise linear 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 first consider the elementary application j(u), where we have dropped the index i compared to [START_REF] Haouchine | Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery[END_REF]. Given a displacement field u, we denote by x 0 and x the initial and current positions of the mesh vertices, respectively. In particular, x 0 and x satisfy x 0 + u = x. Assume that the projection point a = p Su (y) falls into the triangle composed of vertices k 1 , k 2 , k 3 . We use the notations

-∂ k j(u) x k x x m y -∂ k j(u) -∂ m j(u) x k x x m y a
a = θ 1 x k 1 + θ 2 x k 2 + θ 3 x k 3 and a 0 = θ 1 x 0 k 1 + θ 2 x 0 k 2 + θ 3 x 0 k 3 ,
where the (nonnegative) barycentric coefficients θ 1 , θ 2 , θ 3 satisfy θ 1 + θ 2 + θ 3 = 1. In particular, a = a 0 + u(a 0 ). For v a perturbation of u, we obtain after adapting ( 5)

dj(u), v = v(a 0 ) • (a -y) = (θ 1 v k 1 + θ 2 v k 2 + θ 3 v k 3 ) • (a -y) .
Therefore, the gradient ∇j(u) satisfies

∀i ∈ {1, 2, 3} ∂ k i j(u) = θ i (a -y) ∀k / ∈ {k 1 , k 2 , k 3 } ∂ k j(u) = 0.
Finally, the gradient ∇J(u) reads

∇J(u) = 1 p p i=1 ∇j i (u).
Here, the elasticity system is hidden in the objective function Φ and has to be interpreted as an equality constraint of the optimization problem. In particular, each evaluation of the objective value Φ(b) requires to solve the elasticity system to compute u b . Also, based on the facts described in Remark 3.2, we chose to ignore the L ∞ constraint on surface forces in the numerical examples. In this case, the feasible set B still defines on which vertices the nodal forces b j are allowed to be nonzero. Now, to solve (10) using a first-order optimization method, computing the objective gradient ∇Φ(b) is also required. The adjoint method exploits the first-order optimality conditions derived in Theorem 3.2, keeping in mind that we only consider the differentiable case here. In a similar fashion to (7), we define the adjoint state p b , solution to the adjoint elasticity system

Ap = ∇J(u b ).
We adapt the proof of Theorem 3.2 to the discrete context. If h is a perturbation of the control b, we denote by w h the associated perturbation of u b , defined by Aw h = h. Now, we note that

∇Φ(b), h = ∇J(u b ), w h = Ap b , w h = p b , Aw h = p b , h , hence ∇Φ(b) = p b .
Algorithm 1 shows the outline of the adjoint procedure. Adjoint methods are convenient as they result in modular implementations, where each part may be handled by a separate entity. In particular, it is easy to switch functionals, mechanical models or optimization solvers, as they only communicate through simple interfaces. Our implementation2 relies on the Numpy framework. The stiffness matrix is assembled using the SOFA finite element software developed by the Inria Mimesis Team [START_REF] Allard | SOFA -an Open Source Framework for Medical Simulation[END_REF], then it is factorized before the procedure starts. Concerning the optimization procedure itself, it is taken care of by off-the-shelf solvers available in the Scipy library, namely a limited-memory quasi-Newton solver [START_REF] Byrd | A Limited Memory Algorithm for Bound Constrained Optimization[END_REF]. We now present a few numerical examples involving the adjoint method. First, we illustrate the convergence properties of the optimization procedure on a toy problem. It is also the occasion to clarify what can or can't be expected from the method in terms of registration accuracy. Then we show an example in a registration scenario involving the Sparse Data Challenge dataset. In our last result, we opt for a neo-Hookean elastic model and take profit of a simple registration scenario to illustrate the feasibility of choosing a nonlinear model. In all the following examples, the deformable object Ω 0 is represented by a tetrahedral mesh, while the observed data Γ is represented by a point cloud.

A toy problem

Let us begin with an unpleasant remark: though registration error might be reduced by taking physical considerations into account, our registration method comes with no guarantee in terms of displacement accuracy. To illustrate this disclaimer, we evaluate the displacement error of the procedure on a toy problem.

We create a truncated sphere mesh with radius 1, where the distance from the sphere centre to the truncating plane is 1/2. The volume mesh contains 10,385 nodes, 53,992 tetrahedral elements and 6,702 triangular faces. The flat surface of the mesh is subject to a Dirichlet condition, while synthetic forces are applied on the round surface to generate a ground truth displacement field u true . A linear elastic model is used, with E = 1 and ν = 0.49. After the elastic deformation is applied to the mesh, the round surface is sampled to create a point cloud of 10,000 points. Figure 3a shows the initial mesh and the synthetic deformation that we generated.

We create very favourable conditions for the registration procedure. First, the point cloud, very dense and noise-free, 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. For the reconstruction, we maintain the Dirichlet boundary condition on the flat surface and allow reconstructed forces to be supported by the whole round surface ∂Ω N . The round part of the mesh surface is also chosen as the matching surface S 0 (i.e. S 0 = ∂Ω N ).

Figure 3b shows the objective value and gradient norm along iterations of the optimization solver. After 176 iterations (80 seconds on Apple M1 Pro with 32 GB RAM), the gradient norm has decreased by 5 orders of magnitude. The objective function keeps decreasing and evaluates to 2 • 10 -7 as the procedure ends, meaning that the quadratic mean of the distance between data points and the deformed surface is approximately 6 • 10 -4 . Figure 4a shows the point cloud, along with a superposition of the true (green) and reconstructed (yellow) surfaces, confirming the good matching between surfaces.

We define the displacement error for each node (including nodes inside the mesh) as the distance between the node position in the reconstructed configuration and its position in the generated ground truth configuration. Despite a tight surface matching, the displacement error (Figure 3c) settles around 10 -1 , with a maximum error of 0.6 (the sphere radius is 1). In Figure 4b, we illustrated the displacement discrepancy between the true (green) and reconstructed (yellow) surfaces, by applying a pattern on the initial mesh surface. The displacement error is plotted in Figure 4c. Without surprise, the displacement error increases as we move away from the Dirichlet boundary condition. The last graph in Figure 3c shows the error on nodal forces for the 2,522 vertices of the Neumann surface. The nodal force errors settle to the same order of magnitude as the nodal forces themselves, suggesting an average 100% error. In general, it should be difficult to give a meaning to the reconstructed distribution, especially when it contains many degrees of freedom.

Should we give up every hope for all that? Actually, the sphere test case is not as easy as it seems. Due to the sphere's symmetry, the registration method has no clue about the best way to match the point cloud. As a consequence, an inaccurate displacement is not penalized as long as the deformed surface matches the point cloud. In the case of an organ, things are different, as the point cloud reflects the shape of the surface part it is supposed to match. In addition, in the application case, some landmarks, detected by the laparoscopic camera [20, Figure 3], may be used to improve registration accuracy.

Sparse Data Challenge dataset

In this section, we present an organ registration result in the context of augmented surgery. This example, developed in our conference paper [29, Section 3.1], involves the Sparse Data Challenge3 dataset. The dataset contains one tetrahedral mesh representing a silicone liver phantom in its initial configuration and 112 point clouds acquired from 112 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-tophysical liver registration sparse data challenge[END_REF]. To generate deformations, the challenge organizers laid the silicone liver on irregular supports on its posterior face [24, Figure 2]. They produced the point clouds by acquiring pictures of the anterior face with a camera. The silicone liver also contains 159 targets in its volume whose position is measured in reference and deformed configurations to establish ground truth data. After we upload our reconstructed mesh positions to the challenge website, the organisers evaluate the target registration error by comparing the positions of targets according to our reconstruction with their positions in the ground truth data. By keeping the ground truth hidden from participants, the challenge organizers eliminate bias associated with knowledge of the true deformations.

To begin with, we perform a rigid alignment between the initial mesh and the point cloud using the Iterative closest point algorithm [START_REF] Besl | Method for registration of 3-D shapes[END_REF], and afterwards we set a homogeneous Dirichlet boundary condition on a small zone in the liver posterior face to enforce the existence and uniqueness of a solution to the elasticity system. The Dirichlet boundary is chosen arbitrarily, as available information does not suggest an area rather than another. Here, the deformation is caused by the irregular support applying contact forces on the posterior surface, and for this reason we label the posterior surface as the loaded surface, where nodal forces are allowed to take nonzero values. On the other hand, the point cloud represents the anterior surface, and in this case, the anterior surface is selected as the matching surface S 0 , which is taken into account when evaluating J. Figure 5a shows the initial mesh with the two labelled surfaces. We run the adjoint method, using a linear elastic model with E = 1 and ν = 0.4, and we stop the procedure after 200 iterations. The resulting deformations in four cases are shown in Figure 5b.

Displacement errors for the whole dataset are gathered in Table 5c. In Figure 5d, they are compared with other submissions displayed on the challenge website. We obtain the second-best result among all submissions, close to the best-performing method [START_REF] Heiselman | Intraoperative Correction of Liver Deformation Using Sparse Surface and Vascular Features via Linearized Iterative Boundary Reconstruction[END_REF], also based on an inverse biomechanical simulation. In the method from VEIDA [START_REF] Jia | Improving Intraoperative Liver Registration in Image-Guided Surgery with Learning-Based Reconstruction[END_REF], a neural network predicts a shape from the point cloud in the form level-set function, and then the liver deformation is estimated by fitting this shape. In the V2SNet method [START_REF] Pfeiffer | Non-Rigid Volume to Surface Registration Using a Data-Driven Biomechanical Model[END_REF], a neural network is trained from synthetic simulations to predict a displacement field in one step from the point cloud. Finally, the BML2 method is not specified by the challenge organizers. A more in-depth analysis of the challenge results is proposed in [START_REF] Heiselman | Comparison study of sparse data-driven soft tissue registration: preliminary results from the image-to-physical liver registration sparse data challenge[END_REF], along with a more detailed description of the compared methods.

Despite a greater sensibility to surface coverage than the leading team, our average registration error remains well below the 5 mm limit required by our clinical partners. This results highlight the interest of methods based on an inverse problem in intra-operative organ registration.

Toward nonlinear mechanics

All theoretical and numerical results presented until now rely on the linear elastic model. However, while linear elasticity is only relevant for small deformations, organs are made of soft living tissues undergoing large deformations, and a successful registration procedure should allow hyperelastic deformations. In this section, we illustrate the extension of our results to nonlinear elasticity on a simple example involving the neo-Hookean model [START_REF] Bonet | Nonlinear continuum mechanics for finite element analysis[END_REF]Section 5.4.3].

The Neo-Hookean model is defined by the local energy function

ψ(∇u) = µ tr(e(u)) -µ ln(det F ) + λ 2 ln(det F ) 2 ,
where F = I + ∇u is the deformation gradient and e(u) = 1 2 F T F -I is the Green-Saint Venant strain tensor. The elastic displacement u g is solution to the partial differential equation

   -div(ψ (∇u)) = 0 in Ω 0 u = 0 on ∂Ω D ψ (∇u) • n = g on ∂Ω N ,
where the tensor ψ (∇u) is also known as the first Piola-Kirchhoff stress tensor.

We consider a liver mesh with 3,046 vertices, embedded with a Neo-Hookean model. Dirichlet boundary conditions are applied at the hepatic vein entry and in the falciform region. We create a synthetic deformation by applying a local force on one lobe of the liver, and we sample the deformed surface to create a point cloud with 500 points. We perform the registration using the same mesh and elastic model. The reconstructed force distributions is restricted to a zone slightly larger than the zone used to generated the deformation.

A few modifications are required in the adjoint pipeline to use hyperelastic models. In particular, the direct problem is now a nonlinear variational equation, denoted It is usually solved using a Newton [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] or quasi-Newton [START_REF] Gelin | Use of quasi-Newton methods for large strain elastic-plastic finite element computations[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] method. In the Neo-Hookean case, though, ensuring stability of the iterative method is critical to avoid generating displacement fields where the elastic residual F(u) cannot be evaluated. For this reason, we use a trust-region Newton method, for which convergence is guaranteed. Note that the adjoint problem is the linear system

F(u) = b.
F (u b )p = ∇J(u b ),
where the residual Jacobian F (u b ) is evaluated at the solution of the direct problem.

Figure 6 show registration examples involving the Neo-Hookean model. is noteworthy that the number of iterations if the optimization solver is very different between cases. Actually, the registration procedure did not stop because the optimality tolerance was met, but because the iterative direct solver accuracy was not sufficient for the optimization solver to keep generating descent directions. In each case, the original deformation is shown side by side with the reconstructed deformation. The arrows represent nodal forces. Note that the reconstructed force distribution is much noisier and messier than the true one, but it somehow exhibits a trend in the right direction.

These results illustrate the feasibility of the optimal control approach for hyperelastic organ registration. Subsequent work should seek to produce a more mature and efficient procedure involving the nonlinear model.

Discussion and Conclusion

The first numerical results from the implementation of the method described in this paper are very encouraging, even when considering nonlinear elasticity models. We conclude by mentioning below some avenues of improvement that are the subject of current and future work.

Let us first comment the choice to consider static elasticity problems involving a Dirichlet boundary condition. The main role of the Dirichlet boundary condition here is to guarantee the well-posed character of the direct and optimal control problems, in particular the existence and uniqueness of a solution displacement u given a force distribution g. Its role is also, to a lesser extent, to take into account anatomical factors that limit the organ 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). However, we believe that the whole procedure could benefit from a finer modelling of those movement restrictions. We plan to no longer consider a Dirichlet hard-constraint, but instead account for the presence of the hepatic vein entry using penalty terms in the cost functional.

On another level, even though we made choices in order to make the algorithm as light as possible, the question of computation time still needs to be addressed. The adjoint method shares processor time with other parts of the pipeline, including point cloud generation from laparoscopic images and displaying the augmented view. In this article, we do not clarify whether 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. For this reason, we are developing an implementation where the iterative nonlinear elastic solver is replaced with an artificial neural network. Preliminary results [START_REF] Odot | Real-time elastic partial shape matching using a neural network-based adjoint method[END_REF] show a significant speed-up of the procedure due to the efficiency of neural networks.

Finally, we are currently thinking about modeling issues, around the choice of cost functionals, in order to improve the robustness against possible measurement errors on the Γ points. Indeed, we have in practice a cloud of points allowing to reconstruct Γ. We would like to introduce a worst-case functional to make the result as insensitive as possible to possible measurement errors on the data.

A Proofs of Theorem 3.1, Proposition 3.2 and Theorem 3.2

In this section we provide details about the results stated in Section 3. 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 ), towards a limit that we denote g. As G M is closed for the weak-topology, g is an element of G M . We denote respectively by u j = u g j the state associated to g j and by u = u g the state associated to g. The proof falls into two parts. We first prove that the sequence (u j ) converges towards u, and then we finish the proof by showing that the convergence is uniform.

The state equations for u j and u yield

∀v ∈ H 1 D (Ω) Ω 0 Aε(u j ) : ε(v) dx = ∂Ω N g j • v ds ∀v ∈ H 1 D (Ω) Ω 0 Aε(u) : ε(v) dx = ∂Ω N g • v ds.
Taking into account the weak-convergence of (g j ) towards g leads to

∀v ∈ H 1 D (Ω) Ω 0 Aε(u j ) : ε(v) dx → Ω 0 Aε(u) : ε(v) dx. (11)
As ∂Ω D has positive measure in ∂Ω, the linear elasticity bilinear form is an inner product on H 1 D (Ω 0 ) [70, Theorems 6.3-3 and 6.3-4], and ( 11) is equivalent to u j u in H 1 D (Ω 0 ). In particular, u j u in L 2 (Ω 0 ).

We now check that (u j ) converges uniformly. Let q the exponent from Proposition 3.1 and denote by q = q/(q -1) its conjugate exponent. Remember that we assumed q > 3. Note that the trace

v| ∂Ω N of a function v ∈ W 1,q D (Ω 0 ) is in L q (∂Ω N ) [71, §5.5, Theorem 1]
. This confirms that an element of L ∞ (∂Ω N ) ⊂ L q (∂Ω N ) represents an element of W -1,q D (Ω 0 ) (with continuous injection). Because elements of (g j ) are in G M , the sequence is bounded in W -1,q D (Ω 0 ). Then, due to Proposition 3.1, the sequence (u j ) is uniformly bounded in W 1,q D (Ω 0 ), and, as a consequence of the Rellich-Kondrachov theorem [72, Theorem 6.3], it converges uniformly, up to a subsequence. We check that the uniform limit is u by noting that uniform convergence implies weak convergence in L 2 .

Finally, as J is continuous on C(Ω 0 ), we obtain that J(u) = lim J(u j ) = min G M Φ, and g is a solution to (3).

A.2 Proof of Proposition 3.2

The continuous version of J reads J(u) = Γ j y (u) dy where j y (u) = 1 2 d 2 (y, S u ).

We begin by studying the properties of the elementary functional j y in a first lemma, and then we finish the proof of Proposition 3.2 using the dominated convergence theorem. Remember that the definition of P y (u) is given in (4). Proof. Recalling that S u = (Id +u)(S 0 ), we write j y (u) = min x∈S 0 1 2 x + u(x) -y 2 .

First, we check Lipschitz continuity. Let u 1 , u 2 ∈ C(Ω 0 ). As S 0 is compact, there is a x 1 ∈ S 0 with j y (u 1 ) = 1 2 x 1 -u 1 (x 1 ) -y 2 . Then, 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, we obtain the Lipschitz constant

L = max y∈Γ max x∈S 0 x -y + 1 2 u 2 + u 1 L ∞ (Ω 0 ) .
To compute directional derivatives, we define the application f : R + × S 0 → R (t, x) →

∂ t f (0, x),
which is the formula we wanted to end up with.

We now use the dominated convergence theorem to finish the proof. Denote by B = B(0, r) a small ball of radius r > 0. As j y is Lipschitz continuous on u + B, there is a constant L independent from y such that ∀v ∈ B ∀t ∈ (0, 1) j y (u + tv) -j y (u) t L v L ∞ (Ω 0 ) Lr.

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 --→ A.3 Proof of Lemma 3.1

Let us construct the solution with the help of an approximation process. We first consider a single point y ∈ Γ, and we denote z = x y + u(x y ) -y ∈ R 3 . We look for an elementary adjoint state p ,y such that ∀h ∈ L q (∂Ω N )

∂Ω N p ,y • h ds = w h (x y ) • z. ( 12
)
We now characterize the cone of admissible directions at g. Let us introduce the active set A g = {x ∈ ∂Ω N | g(x) = M }. As the pointwise norm g(x) should not exceed M on ∂Ω N , admissible directions h ∈ L ∞ (∂Ω N ) are characterized (see for instance [START_REF] Bednarczuk | Tangent Sets in Some Functional Spaces[END_REF]) by h(x) • g(x) 0 almost everywhere in A g .

An alternative formulation for 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 + ) Ag λ(h • g) ds 0 .
Let us briefly check this expression. First, it is clear that a direction h which satisfies h • g 0 almost everywhere in A g is in K. On the other hand, if there is a subset S ⊂ A g with positive Lebesgue measure where h • g > 0, then, by choosing λ = 1 S , we obtain h / ∈ K. We notice that this new expression for K is the expression of a polar cone, namely

K = Q • = h ∈ L 2 (∂Ω N ) ∀q ∈ Q h, q 0 ,
where the closed cone Q is defined by

Q = 1 Ag λg, λ ∈ L 2 (∂Ω N , R + ) .
We denote by -p = p K + p Q the so-called Moreau decomposition of -p [START_REF] Moreau | Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires[END_REF] with p K ∈ K, p Q ∈ Q and p K , p Q = 0. Now coming back to the optimality condition, as g is a local minimizer of Φ, the Euler inequation stipulates that dΦ(g)(h) 0 for every element h of the cone of admissible directions. In particular, we consider the admissible direction

h = p K max(1, p K ) ∈ L ∞ (∂Ω N ) ∩ K.
It results from the Euler inequation that 0 dΦ(g)(h)

∂Ω N p • h ds = - ∂Ω N min p K , p K 2 ds.
Finally, we obtain p K (•) = 0, which yields the first-order optimality condition p + p Q = 0. Observing that this condition rewrites as in the statement of Theorem 3.2, we are done.

Figure 1 :

 1 Figure 1: Problem geometry. The thick dashed line represents the deformed matching surface S u . The distance between a point y ∈ Γ and its orthogonal projection onto S u is denoted by d(y, S u ).

1 2

 1 ∇u + ∇u T and Aε = 2µε + λ tr(ε)I, respectively. The Lamé parameters µ and λ describe the material resistance to deformation and volume change. They are defined from the Young modulus E and the Poisson ratio ν by µ = E/(2(1 + ν)) and λ = 2νµ/(1 -2ν).

Theorem 3 . 1 (

 31 Proof in Appendix A A.1). Assume that Ω 0 ∪ ∂Ω N is Gröger-regular and let M > 0.1. Let us consider the equivalent of Problem (3) in R 2 . Then, Problem (3) has at least one solution.

Figure 2 :

 2 Figure 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.

Figure 2

 2 Figure2illustrates the components of ∇j(u) in a two-dimensional configuration.

Algorithm 1 :

 1 Computation of the objective gradient using an adjoint method. Data: Current iterate b Compute the displacement u b by solving Au = b Evaluate J(u b ) and ∇J(u b ) Compute the adjoint state p b by solving Ap = ∇J(u b ) Result: ∇Φ(b) = p b 4 Numerical results

  (a) Initial sphere mesh (left) and synthetic deformation (right). The matching surface S 0 is the round surface of the initial mesh. It is transformed into S utrue in the right figure.

  Objective value Φ(b) and gradient norm ∇Φ(b) along iterations.

  Reconstruction error statistics for vertex displacements in the whole volume (left) and nodal forces on the mesh surface (right).

Figure 3 :

 3 Figure 3: Synthetic deformation and convergence statistics for the toy problem.

  (a) Left: Point cloud Γ sampled from the true deformed surface S utrue . Right: superposition of the true surface and the reconstructed surface. (b) Discrepancy between the true (left) and reconstructed (right) surfaces. Despite both surfaces having the same shape, the true displacement field was not recovered, as pattern elements were not transported at the same place on the left and on the right. (c) Node position error in the mesh volume between the ground truth and the reconstructed displacement field.

Figure 4 :

 4 Figure 4: Registration results for the toy problem. The yellow surface corresponds to the reconstructed displacement field, while the green surface corresponds to the true generated displacement field.

  (a) Initial liver mesh. The bottom surface (in blue) is the loaded surface , where forces are applied, while the top surface (in magenta) is the matching surface, that should match with the point cloud.(b) Observed point cloud (magenta) and deformed mesh computed by our algorithm for four cases.

  (d) Comparison of our registration results with other participants (our column is the second one).

Figure 5 :

 5 Figure 5: Results for the Sparse Data Challenge dataset [29].

Figure 6 :

 6 Figure 6: Convergence statistics (top), along with the true and reconstructed deformations (bottom) for two test cases with the nonlinear model.

Lemma A. 1 .

 1 Let y a fixed point in Γ. The application j y is locally Lipschitz continuous on C(Ω 0 ). The Lipschitz constant can be chosen independent from y but still depends on Γ. In addition, for u ∈ C(Ω 0 ), j y has directional derivatives. The directional derivative in the direction v ∈ C(Ω 0 ) reads dj y (u)(v) = min x∈Py(u) v(x) • (x + u(x) -y).

  dj y (u)(v) dy to obtain the derivative. Now, if P y (u) is a singleton {x y } for almost every y ∈ Γ, we obtaindJ(u), v = Γ v(x y ) • (x y + u(x y ) -y) dy |Γ| max y∈Γ d(y, S u ) v L ∞ (Ω 0 ) ,i.e. dJ(u) is a continuous linear form.

  1 2 x + u(x) + tv(x) -y 2 , and we note that j y (u + tv) = min x∈Py(u+tv) f (t, x). The partial derivative ∂ t f exists and is continous, and for t > 0 the set P y (u + tv) is nonempty and compact. All conditions are gathered to apply Danskin's theorem [53, Chapter 3, Theorem I], which stipulates that j y is differentiable at u in the direction v, with dj

y (u)(v) = lim t 0 j y (u + tv) -j y (u) t = min x∈Py(u)

Website: https://trimsh.org

Our code is provided as supplementary material.

Challenge website: sparsedatachallenge.org. Accessed September 5th, 2023.

The regularity assumption on ∂Ω N yields to the existence of a family of bounded open subsets ω δ of ∂Ω N , and such that (ω δ ) δ>0 shrinks to {x y } as δ 0. Setting f δ = 1 ω δ /H 2 (ω δ ), where H 2 (ω δ ) is the two-dimensional Hausdorff measure of ω δ , we introduce p δ as the solution of the PDE

According to Proposition 3.1, one has p δ ∈ W 1,q D (Ω 0 ), and using the Green formula, one gets

Combining this identity to Proposition 3.1 yields

where the constants C and C are independent from δ (and from y). This means that (p δ ) δ>0 is bounded as a family in L q (∂Ω N ). One can thus extract from (p δ | ∂Ω N ) δ>0 a subsequence (p k ) k∈N converging to some p * ∈ L q (∂Ω N ) for the weak-topology of L q (∂Ω N ). Since

we conclude that p ,y = p * is a solution to [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. Now, if we define p = Γ p ,y dy, we can use the Tonelli and Fubini theorems [49, Theorems IV.4 and IV.5] to check that ∀h ∈ L q (∂Ω N )

which makes p a solution to the considered PDE. Finally, regarding the uniqueness, one easily gets that if p 1 and p 2 are two solutions of the PDE, then

for all h ∈ L q (∂Ω N ), hence the result.

A.4 Proof of Theorem 3.2

We first give a justification for the structure of dΦ(g). For a direction h ∈ L ∞ (∂Ω N ), we consider the associated displacement field perturbation w h , which satisfies

By exploiting the notion of solution in the sense of transposition for Problem [START_REF] Marchesseau | Nonlinear Biomechanical Model of the Liver[END_REF], we obtain for ∈ L(u)

, w h =