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

Stéphane Cotin, Guillaume Mestdagh, Yannick Privat

► To cite this version:

Stéphane Cotin, Guillaume Mestdagh, Yannick Privat. Organ registration from partial surface data in augmented surgery from an optimal control perspective. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2023, 480, pp.20230197. 10.1098/rspa.2023.0197. hal-04043695v2

HAL Id: hal-04043695

https://hal.science/hal-04043695v2

Submitted on 11 Jan 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Organ registration from partial surface data in augmented surgery from an optimal control perspective

Stéphane Cotin, Guillaume Mestdagh, Yannick Privat

► To cite this version:

Stéphane Cotin, Guillaume Mestdagh, Yannick Privat. Organ registration from partial surface data in augmented surgery from an optimal control perspective. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2024, 480, pp.20230197. 10.1098/rspa.2023.0197. hal-04043695v2

HAL Id: hal-04043695

https://hal.science/hal-04043695v2

Submitted on 11 Jan 2024

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

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

<figure><img src="image_1.png" /></figure>
# Organ registration from partial surface data in augmented surgery from an optimal control perspective

Stéphane Cotin*¹, Guillaume Mestdagh†¹,², and Yannick Privat‡³,⁴

¹Inria, Strasbourg, France.

²Inria, Lyon, France.

³Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France.

⁴Institut Universitaire de France (IUF).

January 11, 2024

## Abstract

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. Physics-based 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.

**Keywords:** Shape registration, Augmented Surgery, Optimal Control

## 1 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[1]. 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

<footnote>* stephane.cotin@inria.fr</footnote>

<footnote>† guillaume.mestdagh@inria.fr</footnote>

<footnote>‡ yannick.privat@unistra.fr</footnote>
# Organ registration from partial surface data in augmented surgery from an optimal control perspective

Stéphane Cotin*¹, Guillaume Mestdagh†¹,², and Yannick Privat‡³,⁴

¹Inria, Strasbourg, France.

²Inria, Lyon, France.

³Université de Lorraine, CNRS, Inria, IECL, F-54000 Nancy, France.

⁴Institut Universitaire de France (IUF).

January 11, 2024

## Abstract

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. Physics-based 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.

**Keywords:** Shape registration, Augmented Surgery, Optimal Control

## 1 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[1]. 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

<footnote>* stephane.cotin@inria.fr</footnote>

<footnote>† guillaume.mestdagh@inria.fr</footnote>

<footnote>‡ yannick.privat@unistra.fr</footnote>
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 [2, 3, 4, 5, 6]. 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 [5], 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 [7], additional stiffness due to blood vessels [8, 9], and boundary conditions [10, 11].

A whole category of existing methods for intra-operative registration inherits from so-called elastic image registration methods [12, 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 [13, 14, 15, 16, 17]. These methods simulate fictitious forces that attract the liver toward the observed surface, in the same fashion as the Iterative Closest Point algorithm [18]. Artificial forces applied to the organ surface include electrostatic forces [19], linear and nonlinear springs [9, 20], or Lagrange multipliers that enforce a so-called sliding constraint between the deformed surface and the observed point cloud [21, 22, 23].

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 [24]. 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 [25, 26] 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 [27], the initial liver pose is estimated by taking into account the effects of gravity and gas insufflation in the patient's abdomen. In [28], 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
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 [2, 3, 4, 5, 6]. 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 [5], 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 [7], additional stiffness due to blood vessels [8, 9], and boundary conditions [10, 11].

A whole category of existing methods for intra-operative registration inherits from so-called elastic image registration methods [12, 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 [13, 14, 15, 16, 17]. These methods simulate fictitious forces that attract the liver toward the observed surface, in the same fashion as the Iterative Closest Point algorithm [18]. Artificial forces applied to the organ surface include electrostatic forces [19], linear and nonlinear springs [9, 20], or Lagrange multipliers that enforce a so-called sliding constraint between the deformed surface and the observed point cloud [21, 22, 23].

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 [24]. 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 [25, 26] 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 [27], the initial liver pose is estimated by taking into account the effects of gravity and gas insufflation in the patient's abdomen. In [28], 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 [29], 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 [30], the authors control the accretion process leading to the shape of a horn. In [31], the forces driving the motion of a viscoplastic material are estimated. The authors in [32] 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 [33, 34]. 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 [35] 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 [36]. 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.

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

### 2.1 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 $\Omega_0 \subset \mathbb{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 $\partial\Omega_0$ falls into two parts, $\partial\Omega_D$ and $\partial\Omega_N$. On $\partial\Omega_D$, a homogeneous Dirichlet boundary condition applies, whereas the Neumann boundary $\partial\Omega_N$ is subject to a surface force distribution $g \in L^2(\partial\Omega_N, \mathbb{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_D^1(\Omega_0) = H_D^1(\Omega_0, \mathbb{R}^3) = \{u \in H^1(\Omega_0, \mathbb{R}^3) \mid u = 0 \text{ on } \partial\Omega_D\},$$

and we denote by $u_g \in H_D^1(\Omega_0)$ the elastic displacement field generated by a given surface force distribution $g$ on $\partial\Omega_N$.

Elastic stress-strain laws used in the literature to describe the behaviour of living tissues range from the linear elastic model [25] to hyperelastic models such as Neo-Hookean [37] or Ogden [11] models (see also [7] and references therein). When computation time is a constraint, the linear co-rotational model is sometimes preferred [38, 20], 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
paper [29], 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 [30], the authors control the accretion process leading to the shape of a horn. In [31], the forces driving the motion of a viscoplastic material are estimated. The authors in [32] 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 [33, 34]. 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 [35] 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 [36]. 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.

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

## 2.1 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 $\Omega_0 \subset \mathbb{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 $\partial\Omega_0$ falls into two parts, $\partial\Omega_D$ and $\partial\Omega_N$. On $\partial\Omega_D$, a homogeneous Dirichlet boundary condition applies, whereas the Neumann boundary $\partial\Omega_N$ is subject to a surface force distribution $g \in L^2(\partial\Omega_N, \mathbb{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_D^1(\Omega_0) = H_D^1(\Omega_0, \mathbb{R}^3) = \{u \in H^1(\Omega_0, \mathbb{R}^3) \mid u = 0 \text{ on } \partial\Omega_D\},$$

and we denote by $u_g \in H_D^1(\Omega_0)$ the elastic displacement field generated by a given surface force distribution $g$ on $\partial\Omega_N$.

Elastic stress-strain laws used in the literature to describe the behaviour of living tissues range from the linear elastic model [25] to hyperelastic models such as Neo-Hookean [37] or Ogden [11] models (see also [7] and references therein). When computation time is a constraint, the linear co-rotational model is sometimes preferred [38, 20], 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
<figure><img src="image_2.png" /><figcaption>Figure 1: Problem geometry. The thick dashed line represents the deformed matching surface $S_u$. The distance between a point $y \in \Gamma$ and its orthogonal projection onto $S_u$ is denoted by $d(y, S_u)$.</figcaption></figure>

equation (PDE)

$$ \left\{ \begin{array}{ll} -\mathrm{div}(A\varepsilon(u)) = 0 & \text{in } \Omega_0 \\ u = 0 & \text{on } \partial\Omega_D \\ A\varepsilon(u) \cdot n = g & \text{on } \partial\Omega_N, \end{array} \right. \qquad (1) $$

where the linearized strain tensor and the Hooke tensor are defined by

$$ \varepsilon(u) = \frac{1}{2}(\nabla u + \nabla u^T) \quad \text{and} \quad A\varepsilon = 2\mu\varepsilon + \lambda \operatorname{tr}(\varepsilon)I, $$

respectively. The Lamé parameters $\mu$ and $\lambda$ describe the material resistance to deformation and volume change. They are defined from the Young modulus $E$ and the Poisson ratio $\nu$ by $\mu = E/(2(1+\nu))$ and $\lambda = 2\nu\mu/(1-2\nu)$.

We denote by $\Omega_u = (\mathrm{Id}+u)(\Omega_0)$ the volume occupied by the deformed organ. Note that we do not assume $(\mathrm{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.

## 2.2 Intra-operative data and surface-matching functional

While the reference configuration $\Omega_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 $\Gamma \subset \mathbb{R}^3$. We denote by $S_0 \subset \partial\Omega_0$ the matching surface, i.e. a part of the initial organ boundary that is supposed to partially match with $\Gamma$ as the registration ends. We assume that $S_0$ is a compact subset of $\partial\Omega_0$ known by the user, and we will use the notation

$$ S_u = (\mathrm{Id} + u)(S_0) \subset (\mathrm{Id} + u)(\partial\Omega_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 $\Gamma$. If no information is available, it is still possible to set $S_0 = \partial\Omega_0$.

To measure the discrepancy between a given displacement and the observed data, we introduce the least-squares function $J: C(\bar{\Omega}_0) \to \mathbb{R}$, defined by

$$ J(u) = \frac{1}{2} \int_{\Gamma} d^2(y, S_u) \, dy, \qquad (2) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Problem geometry. The thick dashed line represents the deformed matching surface $S_u$. The distance between a point $y \in \Gamma$ and its orthogonal projection onto $S_u$ is denoted by $d(y, S_u)$.</figcaption></figure>

equation (PDE)

$$ \left\{ \begin{array}{ll} -\mathrm{div}(A\varepsilon(u)) = 0 & \text{in } \Omega_0 \\ u = 0 & \text{on } \partial\Omega_D \\ A\varepsilon(u) \cdot n = g & \text{on } \partial\Omega_N, \end{array} \right. \qquad (1) $$

where the linearized strain tensor and the Hooke tensor are defined by

$$ \varepsilon(u) = \frac{1}{2}(\nabla u + \nabla u^T) \quad \text{and} \quad A\varepsilon = 2\mu\varepsilon + \lambda \operatorname{tr}(\varepsilon)I, $$

respectively. The Lamé parameters $\mu$ and $\lambda$ describe the material resistance to deformation and volume change. They are defined from the Young modulus $E$ and the Poisson ratio $\nu$ by $\mu = E/(2(1+\nu))$ and $\lambda = 2\nu\mu/(1-2\nu)$.

We denote by $\Omega_u = (\mathrm{Id}+u)(\Omega_0)$ the volume occupied by the deformed organ. Note that we do not assume $(\mathrm{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.

## 2.2 Intra-operative data and surface-matching functional

While the reference configuration $\Omega_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 $\Gamma \subset \mathbb{R}^3$. We denote by $S_0 \subset \partial\Omega_0$ the matching surface, i.e. a part of the initial organ boundary that is supposed to partially match with $\Gamma$ as the registration ends. We assume that $S_0$ is a compact subset of $\partial\Omega_0$ known by the user, and we will use the notation

$$ S_u = (\mathrm{Id} + u)(S_0) \subset (\mathrm{Id} + u)(\partial\Omega_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 $\Gamma$. If no information is available, it is still possible to set $S_0 = \partial\Omega_0$.

To measure the discrepancy between a given displacement and the observed data, we introduce the least-squares function $J: C(\bar{\Omega}_0) \to \mathbb{R}$, defined by

$$ J(u) = \frac{1}{2} \int_{\Gamma} d^2(y, S_u) \, dy, \qquad (2) $$
where $d(y, S_u)^2 = \min_{x \in S_u} \|x - y\|^2$ denote the Euclidean distance between $y \in \Gamma$ and $S_u$. Therefore, $J(u)$ vanishes whenever $\Gamma \subset S_u$, up to a zero Lebesgue measure set. There is a whole zoo of cost functions for partial shape-matching (see the reviews [39, 40] and reference therein), sometimes involving frameworks such as varifolds [41, 5] or functional maps [42]. 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 (2), which is a continuous version of functions previously used in augmented surgery [24, 21]. Although it is not Gateau-differentiable and requires a special treatment when computing optimality conditions, this function can be quickly evaluated when $\Gamma$ is a point cloud, and therefore is more relevant in the context of augmented surgery.

We consider the constrained optimal control problem

$$ \min_{g \in \mathcal{G}_M} \Phi(g) \quad \text{where} \quad \Phi(g) = J(u_g), \qquad (3) $$

where $g$ and $u_g$ play the role of the control and the state, respectively. The feasible set reads

$$ \mathcal{G}_M = \{ g \in L^\infty(\partial\Omega_N, \mathbb{R}^3) \mid \|g\|_{L^\infty(\partial\Omega_N)} < M \}. $$

Here, the definition of $\mathcal{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 $\Gamma$. In the medical application, $\Gamma$ 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 $\Gamma$ 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 $\Gamma$.

# 3 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_D^{1,q}(\Omega_0, \mathbb{R}^3)$ to denote the Sobolev space

$$ W_D^{1,q}(\Omega_0, \mathbb{R}^3) = \{ u \in W^{1,q}(\Omega_0, \mathbb{R}^3) \mid u = 0 \text{ on } \partial\Omega_D \}, $$

and $q' > 0$ is the conjugate exponent to $q$, given by $1/q + 1/q' = 1$. Then, the topological dual of $W_D^{1,q'}(\Omega_0, \mathbb{R}^3)$ with respect to the pivot space $L^2(\Omega_0, \mathbb{R}^3)$ is denoted by $W_D^{-1,q'}(\Omega_0, \mathbb{R}^3)$.

## 3.1 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 expected for numerical solutions.
where $d(y, S_u)^2 = \min_{x \in S_u} \|x - y\|^2$ denote the Euclidean distance between $y \in \Gamma$ and $S_u$. Therefore, $J(u)$ vanishes whenever $\Gamma \subset S_u$, up to a zero Lebesgue measure set. There is a whole zoo of cost functions for partial shape-matching (see the reviews [39, 40] and reference therein), sometimes involving frameworks such as varifolds [41, 5] or functional maps [42]. 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 (2), which is a continuous version of functions previously used in augmented surgery [24, 21]. Although it is not Gateau-differentiable and requires a special treatment when computing optimality conditions, this function can be quickly evaluated when $\Gamma$ is a point cloud, and therefore is more relevant in the context of augmented surgery.

We consider the constrained optimal control problem

$$ \min_{g \in \mathcal{G}_M} \Phi(g) \quad \text{where} \quad \Phi(g) = J(u_g), \qquad (3) $$

where $g$ and $u_g$ play the role of the control and the state, respectively. The feasible set reads

$$ \mathcal{G}_M = \{ g \in L^\infty(\partial\Omega_N, \mathbb{R}^3) \mid \|g\|_{L^\infty(\partial\Omega_N)} < M \}. $$

Here, the definition of $\mathcal{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 $\Gamma$. In the medical application, $\Gamma$ 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 $\Gamma$ 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 $\Gamma$.

# 3 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_D^{1,q}(\Omega_0, \mathbb{R}^3)$ to denote the Sobolev space

$$ W_D^{1,q}(\Omega_0, \mathbb{R}^3) = \{ u \in W^{1,q}(\Omega_0, \mathbb{R}^3) \mid u = 0 \text{ on } \partial\Omega_D \}, $$

and $q' > 0$ is the conjugate exponent to $q$, given by $1/q + 1/q' = 1$. Then, the topological dual of $W_D^{1,q'}(\Omega_0, \mathbb{R}^3)$ with respect to the pivot space $L^2(\Omega_0, \mathbb{R}^3)$ is denoted by $W_D^{-1,q'}(\Omega_0, \mathbb{R}^3)$.

## 3.1 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 expected for numerical solutions.
The following results are based on the work by Gröger [43] 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 $\mathbb{R}^3$ (with the notation $x = (x_1, x_2, x_3)$):

$$B = \{x \in \mathbb{R}^3 \mid \|x\| < 1\}, \quad B_+ = \{x \in B \mid x_3 > 0\}, \quad D = \{x \in B \mid x_3 = 0\}, \quad D_0 = \{x \in D \mid x_1 < 0\}.$$

**Definition 3.1.** Let $\Omega_0$ be a bounded open subset of $\mathbb{R}^3$ and let $\partial\Omega_N$ be a relatively open part of its boundary. The set $G = \Omega_0 \cup \partial\Omega_N$ is said to be regular in the sense of Gröger if, for each $x \in \partial\Omega_0$, there is a neighbourhood $U$ of $x$ and a Lipschitz diffeomorphism $\Psi: U \to B$ such that $\Psi(U \cap G)$ either coincides with $B_+$, or $B_+ \cup D$, or $B_+ \cup 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 $\Omega_0$ is a bounded domain of class $W^{1,\infty}$, $\partial\Omega_D = \partial\Omega_0 \setminus \partial\Omega_N$ is closed, and the boundary $\partial(\partial\Omega_N)$, seen as a subset of $\partial\Omega_0$, is $W^{1,\infty}$.

Let us mention [43, 45, 46, 44], investigating the $W^{1,q}$ or $C^{0,\alpha}$ 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 [47], $W^{1,q}$ regularity is obtained for solutions to the linear elastic system on a domain of class $C^1$. In [48], 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 $\varepsilon(u)$. We state their result in the specific case of the linear elasticity system.

**Proposition 3.1.** [48, Theorem 1.1] Assume that $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular, and that the Lamé coefficients $\mu, \lambda$ satisfy $\mu > 0$ and $2\mu + 3\lambda > 0$. Then there exists a $q_0 > 2$ such that, for all $q \in [2, q_0]$, the solution to (1) is in $W_D^{1,q}(\Omega_0, \mathbb{R}^3)$ provided that $g$ defines an element of $W_D^{-1,q}(\Omega_0, \mathbb{R}^3)$. In addition, there exists $C > 0$, depending only on $\Omega_0$ and $q$, such that

$$\|u\|_{W^{1,q}(\Omega_0)} \le C \|g\|_{W^{-1,q}(\Omega_0)}.$$

Our existence result relies on the compact embedding $W^{1,q}(\Omega_0) \hookrightarrow C(\bar{\Omega}_0)$ for $q > \dim \mathbb{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 [50, Appendix] and [51] 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 $\mathbb{R}^2$, such an assumption is not necessary as $q_0 > \dim \mathbb{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,\alpha}$ result from [46]. 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 [52, Chapter 4] for additional explanations.

**Theorem 3.1** (Proof in Appendix A A.1). Assume that $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular and let $M > 0$.

1. Let us consider the equivalent of Problem (3) in $\mathbb{R}^2$. Then, Problem (3) has at least one solution.

2. In $\mathbb{R}^3$, assume that $\Omega_0$ and $\partial\Omega_N$ are such that $q_0 > 3$ with the notations of Proposition 3.1. Then Problem (3) has a solution.

## 3.2 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
The following results are based on the work by Gröger [43] 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 $\mathbb{R}^3$ (with the notation $x = (x_1, x_2, x_3)$):

$$B = \{x \in \mathbb{R}^3 \mid \|x\| < 1\}, \quad B_+ = \{x \in B \mid x_3 > 0\}, \quad D = \{x \in B \mid x_3 = 0\}, \quad D_0 = \{x \in D \mid x_1 < 0\}.$$

**Definition 3.1.** Let $\Omega_0$ be a bounded open subset of $\mathbb{R}^3$ and let $\partial\Omega_N$ be a relatively open part of its boundary. The set $G = \Omega_0 \cup \partial\Omega_N$ is said to be regular in the sense of Gröger if, for each $x \in \partial\Omega_0$, there is a neighbourhood $U$ of $x$ and a Lipschitz diffeomorphism $\Psi: U \to B$ such that $\Psi(U \cap G)$ either coincides with $B_+$, or $B_+ \cup D$, or $B_+ \cup 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 $\Omega_0$ is a bounded domain of class $W^{1,\infty}$, $\partial\Omega_D = \partial\Omega_0 \setminus \partial\Omega_N$ is closed, and the boundary $\partial(\partial\Omega_N)$, seen as a subset of $\partial\Omega_0$, is $W^{1,\infty}$.

Let us mention [43, 45, 46, 44], investigating the $W^{1,q}$ or $C^{0,\alpha}$ 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 [47], $W^{1,q}$ regularity is obtained for solutions to the linear elastic system on a domain of class $C^1$. In [48], 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 $\varepsilon(u)$. We state their result in the specific case of the linear elasticity system.

**Proposition 3.1.** [48, Theorem 1.1] Assume that $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular, and that the Lamé coefficients $\mu, \lambda$ satisfy $\mu > 0$ and $2\mu + 3\lambda > 0$. Then there exists a $q_0 > 2$ such that, for all $q \in [2, q_0]$, the solution to (1) is in $W_D^{1,q}(\Omega_0, \mathbb{R}^3)$ provided that $g$ defines an element of $W_D^{-1,q}(\Omega_0, \mathbb{R}^3)$. In addition, there exists $C > 0$, depending only on $\Omega_0$ and $q$, such that

$$\|u\|_{W^{1,q}(\Omega_0)} \le C \|g\|_{W^{-1,q}(\Omega_0)}.$$

Our existence result relies on the compact embedding $W^{1,q}(\Omega_0) \hookrightarrow C(\bar{\Omega}_0)$ for $q > \dim \mathbb{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 [50, Appendix] and [51] 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 $\mathbb{R}^2$, such an assumption is not necessary as $q_0 > \dim \mathbb{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,\alpha}$ result from [46]. 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 [52, Chapter 4] for additional explanations.

**Theorem 3.1** (Proof in Appendix A A.1). Assume that $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular and let $M > 0$.

1. Let us consider the equivalent of Problem (3) in $\mathbb{R}^2$. Then, Problem (3) has at least one solution.

2. In $\mathbb{R}^3$, assume that $\Omega_0$ and $\partial\Omega_N$ are such that $q_0 > 3$ with the notations of Proposition 3.1. Then Problem (3) has a solution.

## 3.2 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 [53, Chapter 3] concerning the differentiability of functions defined by a minimum. For a given displacement field $u \in C(\bar{\Omega}_0)$, we use the notation

$$P_y(u) = \{x \in S_0 \mid \|x + u(x) - y\| = d(y, S_u)\} = (\mathrm{Id} + u)^{-1} \Pi_{S_u}(y), \quad (4)$$

where $\Pi_{S_u}(y)$ is the set of orthogonal projections of $y$ onto $S_u$. In particular, $P_y(u) \subset S_0$.

**Proposition 3.2 (Proof in Appendix A A.2).** Assume that $\Omega_0$ is bounded with Lipschitz boundary, and let $u \in C(\bar{\Omega}_0)$. Then the functional J has directional derivatives at u. Its derivative in the direction $v \in C(\bar{\Omega}_0)$ reads

$$dJ(u)(v) = \int_{\Gamma} \min_{x \in P_y(u)} [v(x) \cdot (x + u(x) - y)] dy. \quad (5)$$

In addition, if $P_y(u)$ is a singleton for almost every point $y \in \Gamma$, then $dJ(u)$ is a continuous linear form on $C(\bar{\Omega}_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 $\Gamma$. 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_{S_u}(y)$ which is the image of several points $x_1, x_2, \dots \in S_0$ under the transformation $\mathrm{Id} + u$.

Note that the directional derivative of J also reads

$$\forall v \in C(\bar{\Omega}_0), \quad dJ(u)(v) = \min_{\ell \in L(u)} \langle \ell, v \rangle, \quad (6)$$

where the set of linear forms $L(u)$ is defined by

$$L(u) = \left\{ \ell : C(\bar{\Omega}_0) \ni v \mapsto \int_{\Gamma} v(x_y) \cdot (x_y + u(x_y) - y) dy \mid x_y \in P_y(u) \right\}.$$

All the linear forms $\ell$ involved in the definition of $L(u)$ belong to the concave subdifferential of J, as they satisfy $dJ(u)(v) \le \ell(v)$ for all $v \in C(\bar{\Omega}_0)$. In addition, according to the expression above, for a given direction $v$, there exists a family $(x_{0,y})_{y \in \Gamma}$, with $x_{0,y} \in P_y(u)$ for all $y \in \Gamma$, such that

$$dJ(u)(v) = \int_{\Gamma} v(x_{0,y}) \cdot (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 $\ell \in L(u_g)$. In particular, the linear form $\ell \in L(u_g)$ associated to the family $(x_y)_{y \in \Gamma}$ is represented in the space of controls by an adjoint state $p_\ell$. From a formal point of view, we aim to introduce $p_\ell$ as the solution to the adjoint problem

$$\left\{ \begin{array}{ll} \mathrm{div}(A\varepsilon(p)) = 0 & \text{in } \Omega_0 \\ p = 0 & \text{on } \partial\Omega_D \\ A\varepsilon(p) \cdot n = \int_{\Gamma}(x_y + u(x_y) - y) \delta_{x_y} dy & \text{on } \partial\Omega_N, \end{array} \right. \quad (7)$$

where $\delta_{x_y}$ is the two-dimensional (pointwise) Dirac measure at $x_y$ on $\partial\Omega_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 $\Omega_0$ and $\partial\Omega_N$ are such that $q_0 > 3$ with the notations of Proposition 3.1, and we consider a fixed $q \in (3, q_0]$.

For a given $h \in L^q(\partial\Omega_N)$, we denote by $w_h$ the solution to the PDE

$$\left\{ \begin{array}{ll} \mathrm{div}(A\varepsilon(w)) = 0 & \text{in } \Omega_0 \\ w = 0 & \text{on } \partial\Omega_D \\ A\varepsilon(w) \cdot n = h & \text{on } \partial\Omega_N. \end{array} \right. \quad (8)$$

In what follows, $q'$ denotes the conjugate exponent to $q$.
J. This result is fully based on Danskin's theorem [53, Chapter 3] concerning the differentiability of functions defined by a minimum. For a given displacement field $u \in C(\bar{\Omega}_0)$, we use the notation

$$P_y(u) = \{x \in S_0 \mid \|x + u(x) - y\| = d(y, S_u)\} = (\mathrm{Id} + u)^{-1} \Pi_{S_u}(y), \quad (4)$$

where $\Pi_{S_u}(y)$ is the set of orthogonal projections of $y$ onto $S_u$. In particular, $P_y(u) \subset S_0$.

**Proposition 3.2 (Proof in Appendix A A.2).** Assume that $\Omega_0$ is bounded with Lipschitz boundary, and let $u \in C(\bar{\Omega}_0)$. Then the functional J has directional derivatives at $u$. Its derivative in the direction $v \in C(\bar{\Omega}_0)$ reads

$$dJ(u)(v) = \min_{\Gamma, x \in P_y(u)} [v(x) \cdot (x + u(x) - y)] dy. \quad (5)$$

In addition, if $P_y(u)$ is a singleton for almost every point $y \in \Gamma$, then $dJ(u)$ is a continuous linear form on $C(\bar{\Omega}_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 $\Gamma$. 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_{S_u}(y)$ which is the image of several points $x_1, x_2, \dots \in S_0$ under the transformation $\mathrm{Id} + u$.

Note that the directional derivative of J also reads

$$\forall v \in C(\bar{\Omega}_0), \quad dJ(u)(v) = \min_{\ell \in L(u)} \langle \ell, v \rangle, \quad (6)$$

where the set of linear forms $L(u)$ is defined by

$$L(u) = \left\{ \ell : C(\bar{\Omega}_0) \ni v \mapsto \int_{\Gamma} v(x_y) \cdot (x_y + u(x_y) - y) dy \mid x_y \in P_y(u) \right\}.$$

All the linear forms $\ell$ involved in the definition of $L(u)$ belong to the concave subdifferential of J, as they satisfy $dJ(u)(v) \le \ell(v)$ for all $v \in C(\bar{\Omega}_0)$. In addition, according to the expression above, for a given direction $v$, there exists a family $(x_{0,y})_{y \in \Gamma}$, with $x_{0,y} \in P_y(u)$ for all $y \in \Gamma$, such that

$$dJ(u)(v) = \int_{\Gamma} v(x_{0,y}) \cdot (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 $\ell \in L(u_g)$. In particular, the linear form $\ell \in L(u_g)$ associated to the family $(x_y)_{y \in \Gamma}$ is represented in the space of controls by an adjoint state $p_\ell$. From a formal point of view, we aim to introduce $p_\ell$ as the solution to the adjoint problem

$$\left\{ \begin{array}{ll} \mathrm{div}(A\varepsilon(p)) = 0 & \text{in } \Omega_0 \\ p = 0 & \text{on } \partial\Omega_D \\ A\varepsilon(p) \cdot n = \int_{\Gamma}(x_y + u(x_y) - y) \delta_{x_y} dy & \text{on } \partial\Omega_N, \end{array} \right. \quad (7)$$

where $\delta_{x_y}$ is the two-dimensional (pointwise) Dirac measure at $x_y$ on $\partial\Omega_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 $\Omega_0$ and $\partial\Omega_N$ are such that $q_0 > 3$ with the notations of Proposition 3.1, and we consider a fixed $q \in (3, q_0]$.

For a given $h \in L^q(\partial\Omega_N)$, we denote by $w_h$ the solution to the PDE

$$\left\{ \begin{array}{ll} \mathrm{div}(A\varepsilon(w)) = 0 & \text{in } \Omega_0 \\ w = 0 & \text{on } \partial\Omega_D \\ A\varepsilon(w) \cdot n = h & \text{on } \partial\Omega_N. \end{array} \right. \quad (8)$$

In what follows, $q'$ denotes the conjugate exponent to $q$.
**Definition 3.2.** One says that $p_{\ell} \in L^q(\partial\Omega_N)$ solves (7) in the sense of transposition if

$$ \int_{\partial\Omega_N} p_{\ell} \cdot h \, ds = \int_{\Gamma} w_h(x_y) \cdot (x_y + u(x_y) - y) \, dy. $$

for all $h \in L^q(\partial\Omega_N)$, where $w_h$ solves (8).

We refer to the monography [54] for additional explanations about solutions to PDEs with measure right-hand sides, as well as the seminal work [55]. Note that, with this definition, the adjoint state $p_{\ell}$ is only defined on $\partial\Omega_N$. Though the values of the adjoint state inside $\Omega_0$ are not needed for our analysis, $p_{\ell}$ can easily be extended to the whole domain $\Omega_0$ using (7).

**Lemma 3.1** (Proof in Appendix A A.3). Let us assume that $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular and let $q \in (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 $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular and $M \in [0, \infty]$. Let $g \in \mathcal{G}_M$ a local minimizer of Problem (3) and $u_g$ the associated displacement field. If $\ell \in L(u)$, denote by $p_{\ell}$ the associated adjoint state defined by (7). For every $\ell \in L(u)$, there exists a Lagrange multiplier $\lambda_{\ell} \in L^2(\partial\Omega_N, \mathbb{R})$, with

$$ \text{for a.e. } x \in \partial\Omega_N \quad \begin{cases} \lambda_{\ell}(x) = 0 & \text{if } \|g(x)\| < M \\ \lambda_{\ell}(x) \ge 0 & \text{if } \|g(x)\| = M, \end{cases} $$

such that $g$ satisfies the first-order optimality condition

$$ \text{for a.e. } x \in \partial\Omega_N \quad p_{\ell}(x) + \lambda_{\ell}(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 \in \mathcal{G}_M$ is a local minimizer of Problem (3) such that $\|g\| < M$ almost everywhere on $\partial\Omega_N$, then $\Phi$ is Gateaux-differentiable at $g$ with $d\Phi(g)(\cdot) = 0$. Indeed, if the $L^\infty$ constraint is inactive, then $\lambda_{\ell} = 0$ by the so-called slackness property. In this case, the optimality condition simply reads $p_{\ell} = 0$ on $\partial\Omega_N$. Using Equation (7), we obtain $\ell = 0$, and therefore $L(u_g) = \{v \mapsto 0\}$, which means that $\Phi$ 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 \in \Gamma$. Indeed, performing a single projection per point $y$ is equivalent to finding one linear form $\ell_0 \in L(u)$. When $u$ is not a local minimizer, $\ell_0 \ne 0$, and if $v$ is chosen so that $\langle \ell_0, v \rangle < 0$, then, using (6),

$$ dJ(u)(v) = \min_{\ell \in L(u)} \langle \ell, v \rangle \le \langle \ell_0, v \rangle < 0, $$

meaning that $v$ is also a descent direction for the whole criterion $J$.

Finally, even when taking into account the $L^\infty$ 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^\infty$ constraint in the optimal control problem definition.

Based on these remarks, we chose to ignore the $L^\infty$ constraint in our numerical examples, and to work with $J$ as if it were Gateaux differentiable everywhere.
**Definition 3.2.** One says that $p_{\ell} \in L^q(\partial\Omega_N)$ solves (7) in the sense of transposition if

$$ \int_{\partial\Omega_N} p_{\ell} \cdot h \, ds = \int_{\Gamma} w_h(x_y) \cdot (x_y + u(x_y) - y) \, dy. $$

for all $h \in L^q(\partial\Omega_N)$, where $w_h$ solves (8).

We refer to the monography [54] for additional explanations about solutions to PDEs with measure right-hand sides, as well as the seminal work [55]. Note that, with this definition, the adjoint state $p_{\ell}$ is only defined on $\partial\Omega_N$. Though the values of the adjoint state inside $\Omega_0$ are not needed for our analysis, $p_{\ell}$ can easily be extended to the whole domain $\Omega_0$ using (7).

**Lemma 3.1** (Proof in Appendix A A.3). Let us assume that $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular and let $q \in (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 $\Omega_0 \cup \partial\Omega_N$ is Gröger-regular and $M \in [0, \infty]$. Let $g \in \mathcal{G}_M$ a local minimizer of Problem (3) and $u_g$ the associated displacement field. If $\ell \in L(u)$, denote by $p_{\ell}$ the associated adjoint state defined by (7). For every $\ell \in L(u)$, there exists a Lagrange multiplier $\lambda_{\ell} \in L^2(\partial\Omega_N, \mathbb{R})$, with

$$ \text{for a.e. } x \in \partial\Omega_N \quad \begin{cases} \lambda_{\ell}(x) = 0 & \text{if } \|g(x)\| < M \\ \lambda_{\ell}(x) \ge 0 & \text{if } \|g(x)\| = M, \end{cases} $$

such that $g$ satisfies the first-order optimality condition

$$ \text{for a.e. } x \in \partial\Omega_N \quad p_{\ell}(x) + \lambda_{\ell}(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 \in \mathcal{G}_M$ is a local minimizer of Problem (3) such that $\|g\| < M$ almost everywhere on $\partial\Omega_N$, then $\Phi$ is Gateaux-differentiable at $g$ with $d\Phi(g)(\cdot) = 0$. Indeed, if the $L^\infty$ constraint is inactive, then $\lambda_{\ell} = 0$ by the so-called slackness property. In this case, the optimality condition simply reads $p_{\ell} = 0$ on $\partial\Omega_N$. Using Equation (7), we obtain $\ell = 0$, and therefore $L(u_g) = \{v \mapsto 0\}$, which means that $\Phi$ 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 \in \Gamma$. Indeed, performing a single projection per point $y$ is equivalent to finding one linear form $\ell_0 \in L(u)$. When $u$ is not a local minimizer, $\ell_0 \ne 0$, and if $v$ is chosen so that $\langle \ell_0, v \rangle < 0$, then, using (6),

$$ dJ(u)(v) = \min_{\ell \in L(u)} \langle \ell, v \rangle \le \langle \ell_0, v \rangle < 0, $$

meaning that $v$ is also a descent direction for the whole criterion $J$.

Finally, even when taking into account the $L^\infty$ 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^\infty$ constraint in the optimal control problem definition.

Based on these remarks, we chose to ignore the $L^\infty$ constraint in our numerical examples, and to work with $J$ as if it were Gateaux differentiable everywhere.
## 3.3 Numerical methods

To solve the registration problem numerically, we discretize the initial domain $\Omega_0$ using a tetrahedral mesh, while displacement fields on $\Omega_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 $\Gamma$, 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 $\Omega_0$ is represented by the finite-dimensional vector $\mathbf{u} = (u_1, \dots, u_n) \in \mathbb{R}^{3n}$ where $n$ is the number of vertices and $u_k \in \mathbb{R}^3$ is the displacement of the $k$-th vertex. The same bold letter $\mathbf{u}$ is used to denote the associated finite element function. We control the vector of nodal forces $\mathbf{b} = (b_1, \dots, b_n)$ at the mesh vertices, and the resulting elastic displacement $\mathbf{u}_b$ is determined by solving the linear system

$$ \mathbf{A}\mathbf{u} = \mathbf{b}, $$

where $\mathbf{A}$ is the stiffness matrix of the liver mesh. Nodal forces $\mathbf{b}$ are related to the surface force distribution $\mathbf{g}$ by the equality $\mathbf{b} = \mathbf{S}\mathbf{g}$, where the matrix $\mathbf{S}$ represents the $L^2$ inner product on $\partial\Omega_N$. Controlling directly the nodal forces $\mathbf{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 $\partial\Omega_0$ is a triangular mesh, and defining $S_0$ consists in selecting the set of triangles that is used to evaluate $J$. The observation $\Gamma$ is provided as the point cloud $\Gamma = \{y_1, \dots, y_p\} \subset \mathbb{R}^3$, and the discretized functional reads

$$ J(\mathbf{u}) = \frac{1}{p} \sum_{i=1}^{p} j_i(\mathbf{u}) \quad \text{where} \quad j_i(\mathbf{u}) = \frac{1}{2} d^2(y_i, S_{\mathbf{u}}). \qquad (9) $$

Our implementation uses the nearest-neighbour search from the Trimesh Python package<fnref n="1" /> to evaluate the distance between a point $y \in \Gamma$ and $S_{\mathbf{u}}$. This projection procedure stores triangles from $S_{\mathbf{u}}$ in a spatial indexing structure [56] 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_{S_{\mathbf{u}}}(y)$ per observed point $y \in \Gamma$. 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_\ell$. 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 $\mathbf{u}$ means differentiating $J$ with respect to the displacements $(u_1, \dots, u_n)$ at the mesh vertices. Thus, the gradient $\nabla J(\mathbf{u})$ is defined by its vertex-wise components $(\partial_1 J(\mathbf{u}), \dots, \partial_n J(\mathbf{u}))$. We first consider the elementary application $j(\mathbf{u})$, where we have dropped the index $i$ compared to (9). Given a displacement field $\mathbf{u}$, we denote by $\mathbf{x}^0$ and $\mathbf{x}$ the initial and current positions of the mesh vertices, respectively. In particular, $\mathbf{x}^0$ and $\mathbf{x}$ satisfy $\mathbf{x}^0 + \mathbf{u} = \mathbf{x}$.

<footnote n="1">¹Website: https://trimsh.org</footnote>
## 3.3 Numerical methods

To solve the registration problem numerically, we discretize the initial domain $\Omega_0$ using a tetrahedral mesh, while displacement fields on $\Omega_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 $\Gamma$, 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 $\Omega_0$ is represented by the finite-dimensional vector $\mathbf{u} = (u_1, \dots, u_n) \in \mathbb{R}^{3n}$ where $n$ is the number of vertices and $u_k \in \mathbb{R}^3$ is the displacement of the $k$-th vertex. The same bold letter $\mathbf{u}$ is used to denote the associated finite element function. We control the vector of nodal forces $\mathbf{b} = (b_1, \dots, b_n)$ at the mesh vertices, and the resulting elastic displacement $\mathbf{u}_b$ is determined by solving the linear system

$$ \mathbf{A}\mathbf{u} = \mathbf{b}, $$

where $\mathbf{A}$ is the stiffness matrix of the liver mesh. Nodal forces $\mathbf{b}$ are related to the surface force distribution $\mathbf{g}$ by the equality $\mathbf{b} = \mathbf{S}\mathbf{g}$, where the matrix $\mathbf{S}$ represents the $L^2$ inner product on $\partial\Omega_N$. Controlling directly the nodal forces $\mathbf{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 $\partial\Omega_0$ is a triangular mesh, and defining $S_0$ consists in selecting the set of triangles that is used to evaluate $J$. The observation $\Gamma$ is provided as the point cloud $\Gamma = \{y_1, \dots, y_p\} \subset \mathbb{R}^3$, and the discretized functional reads

$$ J(\mathbf{u}) = \frac{1}{p} \sum_{i=1}^{p} j_i(\mathbf{u}) \quad \text{where} \quad j_i(\mathbf{u}) = \frac{1}{2} d^2(y_i, S_{\mathbf{u}}). \qquad (9) $$

Our implementation uses the nearest-neighbour search from the Trimesh Python package<fnref n="1" /> to evaluate the distance between a point $y \in \Gamma$ and $S_{\mathbf{u}}$. This projection procedure stores triangles from $S_{\mathbf{u}}$ in a spatial indexing structure [56] 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_{S_{\mathbf{u}}}(y)$ per observed point $y \in \Gamma$. 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_\ell$. 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 $\mathbf{u}$ means differentiating $J$ with respect to the displacements $(u_1, \dots, u_n)$ at the mesh vertices. Thus, the gradient $\nabla J(\mathbf{u})$ is defined by its vertex-wise components $(\partial_1 J(\mathbf{u}), \dots, \partial_n J(\mathbf{u}))$. We first consider the elementary application $j(\mathbf{u})$, where we have dropped the index $i$ compared to (9). Given a displacement field $\mathbf{u}$, we denote by $\mathbf{x}^0$ and $\mathbf{x}$ the initial and current positions of the mesh vertices, respectively. In particular, $\mathbf{x}^0$ and $\mathbf{x}$ satisfy $\mathbf{x}^0 + \mathbf{u} = \mathbf{x}$.

<footnote n="1">¹Website: https://trimsh.org</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of $\nabla j(\mathbf{u})$ in two configurations. Points $x_\ell$, $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 $-\nabla j(\mathbf{u})$. On the right, the equality $y - a = -\partial_k j(\mathbf{u}) - \partial_m j(\mathbf{u})$ holds.</figcaption></figure>

Assume that the projection point $a = p_{S_u}(y)$ falls into the triangle composed of vertices $k_1, k_2, k_3$. We use the notations

$$a = \theta_1 x_{k_1} + \theta_2 x_{k_2} + \theta_3 x_{k_3} \quad \text{and} \quad a^0 = \theta_1 x_{k_1}^0 + \theta_2 x_{k_2}^0 + \theta_3 x_{k_3}^0,$$

where the (nonnegative) barycentric coefficients $\theta_1, \theta_2, \theta_3$ satisfy $\theta_1 + \theta_2 + \theta_3 = 1$. In particular, $a = a^0 + \mathbf{u}(a^0)$. For $\mathbf{v}$ a perturbation of $\mathbf{u}$, we obtain after adapting (5)

$$\langle dj(\mathbf{u}), \mathbf{v} \rangle = \mathbf{v}(a^0) \cdot (a - y) = (\theta_1 v_{k_1} + \theta_2 v_{k_2} + \theta_3 v_{k_3}) \cdot (a - y).$$

Therefore, the gradient $\nabla j(\mathbf{u})$ satisfies

$$\forall i \in \{1, 2, 3\} \quad \partial_{k_i} j(\mathbf{u}) = \theta_i (a - y) \quad \forall k_i \notin \{k_1, k_2, k_3\} \quad \partial_k j(\mathbf{u}) = 0.$$

Finally, the gradient $\nabla J(\mathbf{u})$ reads

$$\nabla J(\mathbf{u}) = \frac{1}{p} \sum_{i=1}^{p} \nabla j_i(\mathbf{u}).$$

Figure 2 illustrates the components of $\nabla j(\mathbf{u})$ in a two-dimensional configuration.

**Adjoint method.** As suggested by the problem formulation in (3), we solve the optimization problem numerically using an adjoint method, where the only variable controlled by the optimization solver is the nodal force distribution $\mathbf{b}$. Keeping the same notations as in (3), the discrete optimization problem reads

$$\min_{\mathbf{b} \in \mathcal{B}} \Phi(\mathbf{b}) \quad \text{where} \quad \Phi(\mathbf{b}) = J(\mathbf{u}_{\mathbf{b}}). \qquad (10)$$

Here, the elasticity system is hidden in the objective function $\Phi$ and has to be interpreted as an equality constraint of the optimization problem. In particular, each evaluation of the objective value $\Phi(\mathbf{b})$ requires to solve the elasticity system to compute $\mathbf{u}_{\mathbf{b}}$.

Also, based on the facts described in Remark 3.2, we chose to ignore the $L^\infty$ constraint on surface forces in the numerical examples. In this case, the feasible set $\mathcal{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 $\nabla\Phi(\mathbf{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 $\mathbf{p}_{\mathbf{b}}$, solution to the adjoint elasticity system

$$\mathbf{A}\mathbf{p} = \nabla J(\mathbf{u}_{\mathbf{b}}).$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of $\nabla j(\mathbf{u})$ in two configurations. Points $x_\ell$, $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 $-\nabla j(\mathbf{u})$. On the right, the equality $y - a = -\partial_k j(\mathbf{u}) - \partial_m j(\mathbf{u})$ holds.</figcaption></figure>

Assume that the projection point $a = p_{S_u}(y)$ falls into the triangle composed of vertices $k_1, k_2, k_3$. We use the notations

$$a = \theta_1 x_{k_1} + \theta_2 x_{k_2} + \theta_3 x_{k_3} \quad \text{and} \quad a^0 = \theta_1 x_{k_1}^0 + \theta_2 x_{k_2}^0 + \theta_3 x_{k_3}^0,$$

where the (nonnegative) barycentric coefficients $\theta_1, \theta_2, \theta_3$ satisfy $\theta_1 + \theta_2 + \theta_3 = 1$. In particular, $a = a^0 + \mathbf{u}(a^0)$. For $\mathbf{v}$ a perturbation of $\mathbf{u}$, we obtain after adapting (5)

$$\langle dj(\mathbf{u}), \mathbf{v} \rangle = \mathbf{v}(a^0) \cdot (a - y) = (\theta_1 v_{k_1} + \theta_2 v_{k_2} + \theta_3 v_{k_3}) \cdot (a - y).$$

Therefore, the gradient $\nabla j(\mathbf{u})$ satisfies

$$\forall i \in \{1, 2, 3\} \quad \partial_{k_i} j(\mathbf{u}) = \theta_i (a - y) \quad \forall k \notin \{k_1, k_2, k_3\} \quad \partial_k j(\mathbf{u}) = 0.$$

Finally, the gradient $\nabla J(\mathbf{u})$ reads

$$\nabla J(\mathbf{u}) = \frac{1}{p} \sum_{i=1}^{p} \nabla j_i(\mathbf{u}).$$

Figure 2 illustrates the components of $\nabla j(\mathbf{u})$ in a two-dimensional configuration.

**Adjoint method.** As suggested by the problem formulation in (3), we solve the optimization problem numerically using an adjoint method, where the only variable controlled by the optimization solver is the nodal force distribution $\mathbf{b}$. Keeping the same notations as in (3), the discrete optimization problem reads

$$\min_{\mathbf{b} \in \mathcal{B}} \Phi(\mathbf{b}) \quad \text{where} \quad \Phi(\mathbf{b}) = J(\mathbf{u}_{\mathbf{b}}). \qquad (10)$$

Here, the elasticity system is hidden in the objective function $\Phi$ and has to be interpreted as an equality constraint of the optimization problem. In particular, each evaluation of the objective value $\Phi(\mathbf{b})$ requires to solve the elasticity system to compute $\mathbf{u}_{\mathbf{b}}$.

Also, based on the facts described in Remark 3.2, we chose to ignore the $L^\infty$ constraint on surface forces in the numerical examples. In this case, the feasible set $\mathcal{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 $\nabla \Phi(\mathbf{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 $\mathbf{p}_{\mathbf{b}}$, solution to the adjoint elasticity system

$$\mathbf{A}\mathbf{p} = \nabla J(\mathbf{u}_{\mathbf{b}}).$$
We adapt the proof of Theorem 3.2 to the discrete context. If $\mathbf{h}$ is a perturbation of the control $\mathbf{b}$, we denote by $\mathbf{w}_h$ the associated perturbation of $\mathbf{u}_b$, defined by $\mathbf{A}\mathbf{w}_h = \mathbf{h}$. Now, we note that

$$ \langle \nabla \Phi(\mathbf{b}), \mathbf{h} \rangle = \langle \nabla J(\mathbf{u}_b), \mathbf{w}_h \rangle = \langle \mathbf{A}\mathbf{p}_b, \mathbf{w}_h \rangle = \langle \mathbf{p}_b, \mathbf{A}\mathbf{w}_h \rangle = \langle \mathbf{p}_b, \mathbf{h} \rangle, $$

hence

$$ \nabla \Phi(\mathbf{b}) = \mathbf{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 implementation<fnref n="2" /> relies on the Numpy framework. The stiffness matrix is assembled using the SOFA finite element software developed by the Inria Mimesis Team [57], 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 [58].

<caption>**Algorithm 1:** Computation of the objective gradient using an adjoint method.</caption>

**Data:** Current iterate $\mathbf{b}$
Compute the displacement $\mathbf{u}_b$ by solving $\mathbf{A}\mathbf{u} = \mathbf{b}$
Evaluate $J(\mathbf{u}_b)$ and $\nabla J(\mathbf{u}_b)$
Compute the adjoint state $\mathbf{p}_b$ by solving $\mathbf{A}\mathbf{p} = \nabla J(\mathbf{u}_b)$
**Result:** $\nabla \Phi(\mathbf{b}) = \mathbf{p}_b$

# 4 Numerical results

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 $\Omega_0$ is represented by a tetrahedral mesh, while the observed data $\Gamma$ is represented by a point cloud.

## 4.1 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_{\text{true}}$. A linear elastic model is used, with $E = 1$ and $\nu = 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.

<footnote n="2">²Our code is provided as supplementary material.</footnote>
We adapt the proof of Theorem 3.2 to the discrete context. If $\mathbf{h}$ is a perturbation of the control $\mathbf{b}$, we denote by $\mathbf{w}_h$ the associated perturbation of $\mathbf{u}_b$, defined by $\mathbf{A}\mathbf{w}_h = \mathbf{h}$. Now, we note that

$$ \langle \nabla \Phi(\mathbf{b}), \mathbf{h} \rangle = \langle \nabla J(\mathbf{u}_b), \mathbf{w}_h \rangle = \langle \mathbf{A}\mathbf{p}_b, \mathbf{w}_h \rangle = \langle \mathbf{p}_b, \mathbf{A}\mathbf{w}_h \rangle = \langle \mathbf{p}_b, \mathbf{h} \rangle, $$

hence

$$ \nabla \Phi(\mathbf{b}) = \mathbf{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 implementation<fnref n="2" /> relies on the Numpy framework. The stiffness matrix is assembled using the SOFA finite element software developed by the Inria Mimesis Team [57], 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 [58].

<caption>**Algorithm 1:** Computation of the objective gradient using an adjoint method.</caption>

**Data:** Current iterate $\mathbf{b}$
Compute the displacement $\mathbf{u}_b$ by solving $\mathbf{A}\mathbf{u} = \mathbf{b}$
Evaluate $J(\mathbf{u}_b)$ and $\nabla J(\mathbf{u}_b)$
Compute the adjoint state $\mathbf{p}_b$ by solving $\mathbf{A}\mathbf{p} = \nabla J(\mathbf{u}_b)$
**Result:** $\nabla \Phi(\mathbf{b}) = \mathbf{p}_b$

# 4 Numerical results

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 $\Omega_0$ is represented by a tetrahedral mesh, while the observed data $\Gamma$ is represented by a point cloud.

## 4.1 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_{\text{true}}$. A linear elastic model is used, with $E = 1$ and $\nu = 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.

<footnote n="2">²Our code is provided as supplementary material.</footnote>
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
∂Ω<sub>N</sub>. The round part of the mesh surface is also chosen as the matching surface S<sub>0</sub> (i.e. S<sub>0</sub> = ∂Ω<sub>N</sub>).

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<sup>-7</sup> as the procedure ends, meaning that the quadratic mean of the distance between data points and the deformed surface is approximately 6 · 10<sup>-4</sup>. 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<sup>-1</sup>, 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.

## 4.2 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 Challenge*<fnref n="3" />
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
[59, 60]. 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 [18], 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

<footnote n="3">³Challenge website: sparsedatachallenge.org. Accessed September 5th, 2023.</footnote>
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
∂Ω<sub>N</sub>. The round part of the mesh surface is also chosen as the matching surface S<sub>0</sub> (i.e. S<sub>0</sub> = ∂Ω<sub>N</sub>).

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<sup>-7</sup> as the procedure ends, meaning that the quadratic mean of the distance between data points and the deformed surface is approximately 6 · 10<sup>-4</sup>. 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<sup>-1</sup>, 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.

## 4.2 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 Challenge*<fnref n="3" /> 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 [59, 60]. 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 [18], 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

<footnote n="3">³Challenge website: sparsedatachallenge.org. Accessed September 5th, 2023.</footnote>
<figure><img src="image_4.png" /><figcaption>(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_{u_{\text{true}}}$ in the right figure.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Reconstruction error statistics for vertex displacements in the whole volume (left) and nodal forces on the mesh surface (right).</figcaption></figure>

<caption>Figure 3: Synthetic deformation and convergence statistics for the toy problem.</caption>
<figure><img src="image_4.png" /><figcaption>(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_{u_{\text{true}}}$ in the right figure.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Reconstruction error statistics for vertex displacements in the whole volume (left) and nodal forces on the mesh surface (right).</figcaption></figure>

<caption>Figure 3: Synthetic deformation and convergence statistics for the toy problem.</caption>
<figure><img src="image_6.png" /><figcaption>(a) Left: Point cloud $\Gamma$ sampled from the true deformed surface $S_{u,\text{true}}$. Right: superposition of the true surface and the reconstructed surface.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Node position error in the mesh volume between the ground truth and the reconstructed displacement field.</figcaption></figure>

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.
<figure><img src="image_6.png" /><figcaption>(a) Left: Point cloud $\Gamma$ sampled from the true deformed surface $S_{u,\text{true}}$. Right: superposition of the true surface and the reconstructed surface.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Node position error in the mesh volume between the ground truth and the reconstructed displacement field.</figcaption></figure>

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.
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 $\nu = 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 [26], also based on an inverse biomechanical simulation. In the method from VEIDA [61], 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 [62], 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 [63], 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.

## 4.3 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 [64, Section 5.4.3].

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

$$ \psi(\nabla u) = \mu \operatorname{tr}(e(u)) - \mu \ln(\det F) + \frac{\lambda}{2} \ln(\det F)^2, $$

where $F = I + \nabla u$ is the deformation gradient and $e(u) = \frac{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

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\psi'(\nabla u)) & = 0 \quad \text{in } \Omega_0 \\ u & = 0 \quad \text{on } \partial\Omega_D \\ \psi'(\nabla u) \cdot n & = g \quad \text{on } \partial\Omega_N, \end{array} \right. $$

where the tensor $\psi'(\nabla 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

$$ F(\mathbf{u}) = \mathbf{b}. $$
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 $\nu = 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 [26], also based on an inverse biomechanical simulation. In the method from VEIDA [61], 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 [62], 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 [63], 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.

## 4.3 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 [64, Section 5.4.3].

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

$$ \psi(\nabla u) = \mu \operatorname{tr}(e(u)) - \mu \ln(\det F) + \frac{\lambda}{2} \ln(\det F)^2, $$

where $F = I + \nabla u$ is the deformation gradient and $e(u) = \frac{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

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\psi'(\nabla u)) & = 0 \quad \text{in } \Omega_0 \\ u & = 0 \quad \text{on } \partial\Omega_D \\ \psi'(\nabla u) \cdot n & = g \quad \text{on } \partial\Omega_N, \end{array} \right. $$

where the tensor $\psi'(\nabla 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

$$ F(\mathbf{u}) = \mathbf{b}. $$
<figure><img src="image_9.png" /><figcaption>(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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Observed point cloud (magenta) and deformed mesh computed by our algorithm for four cases.</figcaption></figure>

<table><caption>(c) Target registration error statistics (in millimeters) for all 112 datasets, as returned by the website after submission.</caption><thead><tr><td>Surface Coverage</td><td>Average</td><td>Standard deviation</td><td>Median</td></tr></thead><tbody><tr><td>20-28 %</td><td>3.54</td><td>1.11</td><td>3.47</td></tr><tr><td>28-36 %</td><td>3.27</td><td>0.85</td><td>3.19</td></tr><tr><td>36-44 %</td><td>3.13</td><td>0.82</td><td>3.13</td></tr><tr><td>All data sets</td><td>3.31</td><td>0.94</td><td>3.19</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>(d) Comparison of our registration results with other participants (our column is the second one).</figcaption></figure>

<caption>Figure 5: Results for the *Sparse Data Challenge* dataset [29].</caption>
<figure><img src="image_9.png" /><figcaption>(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.</figcaption></figure>

<table><caption>(b) Observed point cloud (magenta) and deformed mesh computed by our algorithm for four cases.</caption><thead><tr><td>Surface Coverage</td><td>Average</td><td>Standard deviation</td><td>Median</td></tr></thead><tbody><tr><td>20-28 %</td><td>3.54</td><td>1.11</td><td>3.47</td></tr><tr><td>28-36 %</td><td>3.27</td><td>0.85</td><td>3.19</td></tr><tr><td>36-44 %</td><td>3.13</td><td>0.82</td><td>3.13</td></tr><tr><td>All data sets</td><td>3.31</td><td>0.94</td><td>3.19</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>(d) Comparison of our registration results with other participants (our column is the second one).</figcaption></figure>

<caption>Figure 5: Results for the *Sparse Data Challenge* dataset [29].</caption>
It is usually solved using a Newton [65, 66] or quasi-Newton [67, 68] 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)p = \nabla 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. It 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.

# 5 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 [69] 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 $\Gamma$ points. Indeed, we have in practice a cloud of points allowing to reconstruct $\Gamma$. 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.
It is usually solved using a Newton [65, 66] or quasi-Newton [67, 68] 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)p = \nabla 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. It 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.

# 5 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 [69] 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 $\Gamma$ points. Indeed, we have in practice a cloud of points allowing to reconstruct $\Gamma$. 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.
<figure><img src="image_12.png" /><figcaption>Figure 6: Convergence statistics (top), along with the true and reconstructed deformations (bottom) for two test cases with the nonlinear model.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 6: Convergence statistics (top), along with the true and reconstructed deformations (bottom) for two test cases with the nonlinear model.</figcaption></figure>
## A.1 Proof of Theorem 3.1

Consider a minimizing sequence $(g_j)$ of elements of $\mathcal{G}_M$. As $(g_j)$ is bounded in $L^\infty(\partial\Omega_N)$, it converges $\star$-weakly in $L^\infty(\partial\Omega_N)$, towards a limit that we denote $g$. As $\mathcal{G}_M$ is closed for the weak-$\star$ topology, $g$ is an element of $\mathcal{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

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(u_j) : \varepsilon(v) \, dx = \int_{\partial\Omega_N} g_j \cdot v \, ds $$

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(u) : \varepsilon(v) \, dx = \int_{\partial\Omega_N} g \cdot v \, ds. $$

Taking into account the weak-$\star$ convergence of $(g_j)$ towards $g$ leads to

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(u_j) : \varepsilon(v) \, dx \to \int_{\Omega_0} A\varepsilon(u) : \varepsilon(v) \, dx. \quad (11) $$

As $\partial\Omega_D$ has positive measure in $\partial\Omega$, the linear elasticity bilinear form is an inner product on $H_D^1(\Omega_0)$ [70, Theorems 6.3-3 and 6.3-4], and (11) is equivalent to $u_j \to u$ in $H_D^1(\Omega_0)$. In particular, $u_j \to u$ in $L^2(\Omega_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|_{\partial\Omega_N}$ of a function $v \in W_D^{1,q'}(\Omega_0)$ is in $L^q(\partial\Omega_N)$ [71, §5.5, Theorem 1]. This confirms that an element of $L^\infty(\partial\Omega_N) \subset L^{q'}(\partial\Omega_N)'$ represents an element of $W_D^{-1,q}(\Omega_0)$ (with continuous injection). Because elements of $(g_j)$ are in $\mathcal{G}_M$, the sequence is bounded in $W_D^{-1,q}(\Omega_0)$. Then, due to Proposition 3.1, the sequence $(u_j)$ is uniformly bounded in $W_D^{1,q}(\Omega_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(\bar{\Omega}_0)$, we obtain that $J(u) = \lim J(u_j) = \min_{\mathcal{G}_M} \Phi$, and $g$ is a solution to (3).

## A.2 Proof of Proposition 3.2

The continuous version of $J$ reads

$$ J(u) = \int_{\Gamma} j_{y}(u) \, dy \quad \text{where} \quad j_{y}(u) = \frac{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).

**Lemma A.1.** Let $y$ a fixed point in $\Gamma$. The application $j_y$ is locally Lipschitz continuous on $C(\bar{\Omega}_0)$. The Lipschitz constant can be chosen independent from $y$ but still depends on $\Gamma$. In addition, for $u \in C(\bar{\Omega}_0)$, $j_y$ has directional derivatives. The directional derivative in the direction $v \in C(\bar{\Omega}_0)$ reads

$$ d j_y(u)(v) = \min_{x \in P_y(u)} v(x) \cdot (x + u(x) - y). $$

*Proof.* Recalling that $S_u = (\mathrm{Id} + u)(S_0)$, we write $j_y(u) = \min_{x \in S_0} \frac{1}{2} \|x + u(x) - y\|^2$.
## A.1 Proof of Theorem 3.1

Consider a minimizing sequence $(g_j)$ of elements of $\mathcal{G}_M$. As $(g_j)$ is bounded in $L^\infty(\partial\Omega_N)$, it converges $\star$-weakly in $L^\infty(\partial\Omega_N)$, towards a limit that we denote $g$. As $\mathcal{G}_M$ is closed for the weak-$\star$ topology, $g$ is an element of $\mathcal{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

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(u_j) : \varepsilon(v) \, dx = \int_{\partial\Omega_N} g_j \cdot v \, ds $$

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(u) : \varepsilon(v) \, dx = \int_{\partial\Omega_N} g \cdot v \, ds. $$

Taking into account the weak-$\star$ convergence of $(g_j)$ towards $g$ leads to

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(u_j) : \varepsilon(v) \, dx \to \int_{\Omega_0} A\varepsilon(u) : \varepsilon(v) \, dx. \quad (11) $$

As $\partial\Omega_D$ has positive measure in $\partial\Omega$, the linear elasticity bilinear form is an inner product on $H_D^1(\Omega_0)$ [70, Theorems 6.3-3 and 6.3-4], and (11) is equivalent to $u_j \to u$ in $H_D^1(\Omega_0)$. In particular, $u_j \to u$ in $L^2(\Omega_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|_{\partial\Omega_N}$ of a function $v \in W_D^{1,q'}(\Omega_0)$ is in $L^q(\partial\Omega_N)$ [71, §5.5, Theorem 1]. This confirms that an element of $L^\infty(\partial\Omega_N) \subset L^{q'}(\partial\Omega_N)'$ represents an element of $W_D^{-1,q}(\Omega_0)$ (with continuous injection). Because elements of $(g_j)$ are in $\mathcal{G}_M$, the sequence is bounded in $W_D^{-1,q}(\Omega_0)$. Then, due to Proposition 3.1, the sequence $(u_j)$ is uniformly bounded in $W_D^{1,q}(\Omega_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(\bar{\Omega}_0)$, we obtain that $J(u) = \lim J(u_j) = \min_{\mathcal{G}_M} \Phi$, and $g$ is a solution to (3).

## A.2 Proof of Proposition 3.2

The continuous version of $J$ reads

$$ J(u) = \int_{\Gamma} j_{y}(u) \, dy \quad \text{where} \quad j_{y}(u) = \frac{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).

**Lemma A.1.** Let $y$ a fixed point in $\Gamma$. The application $j_y$ is locally Lipschitz continuous on $C(\bar{\Omega}_0)$. The Lipschitz constant can be chosen independent from $y$ but still depends on $\Gamma$. In addition, for $u \in C(\bar{\Omega}_0)$, $j_y$ has directional derivatives. The directional derivative in the direction $v \in C(\bar{\Omega}_0)$ reads

$$ \mathrm{d}j_y(u)(v) = \min_{x \in P_y(u)} v(x) \cdot (x + u(x) - y). $$

*Proof.* Recalling that $S_u = (\mathrm{Id} + u)(S_0)$, we write $j_y(u) = \min_{x \in S_0} \frac{1}{2} \|x + u(x) - y\|^2$.
First, we check Lipschitz continuity. Let $u_1, u_2 \in C(\bar{\Omega}_0)$. As $S_0$ is compact, there is a $x_1 \in S_0$ with
$j_y(u_1) = \frac{1}{2}\|x_1 - u_1(x_1) - y\|^2$. Then,

$$
\begin{aligned}
j_y(u_2) - j_y(u_1) &\le \frac{1}{2}\|x_1 + u_2(x_1) - y\|^2 - \frac{1}{2}\|x_1 + u_1(x_1) - y\|^2 \\
&= \frac{1}{2}(2x_1 + u_2(x_1) + u_1(x_1) - 2y) \cdot (u_2(x_1) - u_1(x_1)) \\
&\le \left( \max_{x \in S_0} \|x - y\| + \frac{1}{2}\|u_2 + u_1\|_{L^\infty(\Omega_0)} \right) \|u_2 - u_1\|_{L^\infty(\Omega_0)},
\end{aligned}
$$

As $\Gamma$ is compact, we obtain the Lipschitz constant

$$
L = \max_{y \in \Gamma} \left( \max_{x \in S_0} \|x - y\| \right) + \frac{1}{2} \|u_2 + u_1\|_{L^{\infty}(\Omega_0)}.
$$

To compute directional derivatives, we define the application

$$
f : \mathbb{R}_+ \times S_0 \to \mathbb{R} \\
(t, x) \mapsto \frac{1}{2} \|x + u(x) + tv(x) - y\|^2,
$$

and we note that $j_y(u + tv) = \min_{x \in P_y(u+tv)} f(t, x)$. The partial derivative $\partial_t f$ exists and is continuous, 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 \searrow 0} \frac{j_y(u + tv) - j_y(u)}{t} = \min_{x \in P_y(u)} \partial_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

$$
\forall v \in B \quad \forall t \in (0, 1) \quad \left| \frac{j_y(u + tv) - j_y(u)}{t} \right| \le L \|v\|_{L^{\infty}(\Omega_0)} \le Lr.
$$

Using the dominated convergence theorem, we pass to the limit in the expression

$$
\frac{J(u + tv) - J(u)}{t} = \int_{\Gamma} \frac{j_y(u + tv) - j_y(u)}{t} dy \xrightarrow{t \searrow 0} \int_{\Gamma} dj_y(u)(v) dy
$$

to obtain the derivative.

Now, if $P_y(u)$ is a singleton $\{x_y\}$ for almost every $y \in \Gamma$, we obtain

$$
\langle dJ(u), v \rangle = \int_{\Gamma} v(x_y) \cdot (x_y + u(x_y) - y) \, dy \le |\Gamma| \left( \max_{y \in \Gamma} d(y, S_u) \right) \|v\|_{L^{\infty}(\Omega_0)}
$$

i.e. $dJ(u)$ is a continuous linear form.

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 \in \Gamma$, and we denote $z = x_y + u(x_y) - y \in \mathbb{R}^3$. We look for an elementary adjoint state $p_{\ell,y}$ such
that

$$
\forall h \in L^q(\partial \Omega_N) \quad \int_{\partial \Omega_N} p_{\ell,y} \cdot h \, ds = w_h(x_y) \cdot z. \quad (12)
$$
First, we check Lipschitz continuity. Let $u_1, u_2 \in C(\bar{\Omega}_0)$. As $S_0$ is compact, there is a $x_1 \in S_0$ with
$j_y(u_1) = \frac{1}{2}\|x_1 - u_1(x_1) - y\|^2$. Then,

$$
\begin{aligned}
j_y(u_2) - j_y(u_1) &\le \frac{1}{2}\|x_1 + u_2(x_1) - y\|^2 - \frac{1}{2}\|x_1 + u_1(x_1) - y\|^2 \\
&= \frac{1}{2}(2x_1 + u_2(x_1) + u_1(x_1) - 2y) \cdot (u_2(x_1) - u_1(x_1)) \\
&\le \left( \max_{x \in S_0} \|x - y\| + \frac{1}{2}\|u_2 + u_1\|_{L^\infty(\Omega_0)} \right) \|u_2 - u_1\|_{L^\infty(\Omega_0)},
\end{aligned}
$$

As $\Gamma$ is compact, we obtain the Lipschitz constant

$$
L = \max_{y \in \Gamma} \left( \max_{x \in S_0} \|x - y\| \right) + \frac{1}{2} \|u_2 + u_1\|_{L^{\infty}(\Omega_0)}.
$$

To compute directional derivatives, we define the application

$$
f: \mathbb{R}_+ \times S_0 \to \mathbb{R} \\
(t, x) \mapsto \frac{1}{2} \|x + u(x) + tv(x) - y\|^2,
$$

and we note that $j_y(u + tv) = \min_{x \in P_y(u+tv)} f(t, x)$. The partial derivative $\partial_t f$ exists and is continuous, 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 \searrow 0} \frac{j_y(u + tv) - j_y(u)}{t} = \min_{x \in P_y(u)} \partial_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

$$
\forall v \in B \quad \forall t \in (0, 1) \quad \left| \frac{j_y(u + tv) - j_y(u)}{t} \right| \le L \|v\|_{L^{\infty}(\Omega_0)} \le Lr.
$$

Using the dominated convergence theorem, we pass to the limit in the expression

$$
\frac{J(u + tv) - J(u)}{t} = \int_{\Gamma} \frac{j_y(u + tv) - j_y(u)}{t} dy \xrightarrow{t \searrow 0} \int_{\Gamma} dj_y(u)(v) dy
$$

to obtain the derivative.

Now, if $P_y(u)$ is a singleton $\{x_y\}$ for almost every $y \in \Gamma$, we obtain

$$
\langle dJ(u), v \rangle = \int_{\Gamma} v(x_y) \cdot (x_y + u(x_y) - y) \, dy \le |\Gamma| \left( \max_{y \in \Gamma} d(y, S_u) \right) \|v\|_{L^{\infty}(\Omega_0)}
$$

i.e. $dJ(u)$ is a continuous linear form.

**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 \in \Gamma$, and we denote $z = x_y + u(x_y) - y \in \mathbb{R}^3$. We look for an elementary adjoint state $p_{\ell,y}$ such
that

$$
\forall h \in L^q(\partial \Omega_N) \quad \int_{\partial \Omega_N} p_{\ell,y} \cdot h \, ds = w_h(x_y) \cdot z. \quad (12)
$$
The regularity assumption on $\partial\Omega_N$ yields to the existence of a family of bounded open subsets $\omega_\delta$ of $\partial\Omega_N$, and such that $(\omega_\delta)_{\delta>0}$ shrinks to $\{x_y\}$ as $\delta \searrow 0$. Setting $f_\delta = \mathbf{1}_{\omega_\delta}/\mathcal{H}^2(\omega_\delta)$, where $\mathcal{H}^2(\omega_\delta)$ is the two-dimensional Hausdorff measure of $\omega_\delta$, we introduce $p_\delta$ as the solution of the PDE

$$ 
\left\{ 
\begin{array}{lll}
\mathrm{div}(A\varepsilon(p_\delta)) & = & 0 \quad \text{in } \Omega_0 \\
p_\delta & = & 0 \quad \text{on } \partial\Omega_N \\
A\varepsilon(p_\delta) \cdot n & = & f_\delta z \quad \text{on } \partial\Omega_N,
\end{array}
\right. 
$$

According to Proposition 3.1, one has $p_\delta \in W_D^{1,q}(\Omega_0)$, and using the Green formula, one gets

$$ \forall h \in L^q(\partial\Omega_N) \quad \int_{\partial\Omega_N} p_\delta \cdot h \,ds = \frac{1}{\mathcal{H}^2(\omega_\delta)} \int_{\omega_\delta} w_h \cdot z \,ds. $$

Combining this identity to Proposition 3.1 yields

$$ \left| \int_{\partial\Omega_N} p_\delta \cdot h \,ds \right| \le \|z\| \|w_h\|_{C(\bar{\Omega}_0)} \le C' \|z\| \|w_h\|_{W^{1,q}(\Omega_0)} \le C'' \|z\| \|h\|_{L^q(\partial\Omega_N)}, $$

where the constants $C'$ and $C''$ are independent from $\delta$ (and from $y$). This means that $(p_\delta)_{\delta>0}$ is bounded as a family in $L^q'(\partial\Omega_N)$. One can thus extract from $(p_\delta|_{\partial\Omega_N})_{\delta>0}$ a subsequence $(p_k)_{k \in \mathbb{N}}$ converging to some $p^* \in L^q'(\partial\Omega_N)$ for the weak-topology of $L^q'(\partial\Omega_N)$. Since

$$ \lim_{\delta \searrow 0} \frac{1}{\mathcal{H}^2(\omega_\delta)} \int_{\omega_\delta} w_h \cdot z \,ds = w_h(x_y) \cdot z, $$

we conclude that $p_{\ell,y} = p^*$ is a solution to (12). Now, if we define $p_\ell = \int_\Gamma p_{\ell,y} \,dy$, we can use the Tonelli and Fubini theorems [49, Theorems IV.4 and IV.5] to check that

$$ \forall h \in L^q(\partial\Omega_N) \quad \int_{\partial\Omega_N} p_\ell \cdot h \,ds = \int_\Gamma \int_{\partial\Omega_N} p_{\ell,y} \cdot h \,ds \,dy = \int_\Gamma w_h(x_y) \cdot (x_y + u(x_y) - y) \,dy, $$

which makes $p_\ell$ 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

$$ \int_{\partial\Omega_N} (p_1 - p_2) \cdot h \,ds = 0 $$

for all $h \in L^q(\partial\Omega_N)$, hence the result.

## A.4 Proof of Theorem 3.2

We first give a justification for the structure of $d\Phi(g)$. For a direction $h \in L^\infty(\partial\Omega_N)$, we consider the associated displacement field perturbation $w_h$, which satisfies

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(w) : \varepsilon(v) \,dx = \int_{\partial\Omega_N} h \cdot v \,ds. $$

By exploiting the notion of solution in the sense of transposition for Problem (7), we obtain for $\ell \in L(u)$

$$ \langle \ell, w_h \rangle = \int_{\partial\Omega_N} p_\ell \cdot h \,ds, $$

hence

$$ d\Phi(g)(h) = \min_{\ell \in L(u)} \int_{\partial\Omega_N} p_\ell \cdot h \,ds. $$
The regularity assumption on $\partial\Omega_N$ yields to the existence of a family of bounded open subsets $\omega_\delta$ of $\partial\Omega_N$, and such that $(\omega_\delta)_{\delta>0}$ shrinks to $\{x_y\}$ as $\delta \searrow 0$. Setting $f_\delta = \mathbf{1}_{\omega_\delta}/\mathcal{H}^2(\omega_\delta)$, where $\mathcal{H}^2(\omega_\delta)$ is the two-dimensional Hausdorff measure of $\omega_\delta$, we introduce $p_\delta$ as the solution of the PDE

$$ 
\left\{ 
\begin{array}{lll}
\mathrm{div}(A\varepsilon(p_\delta)) & = & 0 \quad \text{in } \Omega_0 \\
p_\delta & = & 0 \quad \text{on } \partial\Omega_N \\
A\varepsilon(p_\delta) \cdot n & = & f_\delta z \quad \text{on } \partial\Omega_N,
\end{array}
\right. 
$$

According to Proposition 3.1, one has $p_\delta \in W_D^{1,q}(\Omega_0)$, and using the Green formula, one gets

$$ \forall h \in L^q(\partial\Omega_N) \quad \int_{\partial\Omega_N} p_\delta \cdot h \, ds = \frac{1}{\mathcal{H}^2(\omega_\delta)} \int_{\omega_\delta} w_h \cdot z \, ds. $$

Combining this identity to Proposition 3.1 yields

$$ \left| \int_{\partial\Omega_N} p_\delta \cdot h \, ds \right| \le \|z\| \|w_h\|_{C(\bar{\Omega}_0)} \le C' \|z\| \|w_h\|_{W^{1,q}(\Omega_0)} \le C'' \|z\| \|h\|_{L^q(\partial\Omega_N)}, $$

where the constants $C'$ and $C''$ are independent from $\delta$ (and from $y$). This means that $(p_\delta)_{\delta>0}$ is bounded as a family in $L^q'(\partial\Omega_N)$. One can thus extract from $(p_\delta|_{\partial\Omega_N})_{\delta>0}$ a subsequence $(p_k)_{k \in \mathbb{N}}$ converging to some $p^* \in L^q'(\partial\Omega_N)$ for the weak-topology of $L^q'(\partial\Omega_N)$. Since

$$ \lim_{\delta \searrow 0} \frac{1}{\mathcal{H}^2(\omega_\delta)} \int_{\omega_\delta} w_h \cdot z \, ds = w_h(x_y) \cdot z, $$

we conclude that $p_{\ell,y} = p^*$ is a solution to (12). Now, if we define $p_\ell = \int_\Gamma p_{\ell,y} \, dy$, we can use the Tonelli and Fubini theorems [49, Theorems IV.4 and IV.5] to check that

$$ \forall h \in L^q(\partial\Omega_N) \quad \int_{\partial\Omega_N} p_\ell \cdot h \, ds = \int_\Gamma \int_{\partial\Omega_N} p_{\ell,y} \cdot h \, ds \, dy = \int_\Gamma w_h(x_y) \cdot (x_y + u(x_y) - y) \, dy, $$

which makes $p_\ell$ 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

$$ \int_{\partial\Omega_N} (p_1 - p_2) \cdot h \, ds = 0 $$

for all $h \in L^q(\partial\Omega_N)$, hence the result.

## A.4 Proof of Theorem 3.2

We first give a justification for the structure of $d\Phi(g)$. For a direction $h \in L^\infty(\partial\Omega_N)$, we consider the associated displacement field perturbation $w_h$, which satisfies

$$ \forall v \in H_D^1(\Omega) \quad \int_{\Omega_0} A\varepsilon(w) : \varepsilon(v) \, dx = \int_{\partial\Omega_N} h \cdot v \, ds. $$

By exploiting the notion of solution in the sense of transposition for Problem (7), we obtain for $\ell \in L(u)$

$$ \langle \ell, w_h \rangle = \int_{\partial\Omega_N} p_\ell \cdot h \, ds, $$

hence

$$ d\Phi(g)(h) = \min_{\ell \in L(u)} \int_{\partial\Omega_N} p_\ell \cdot h \, ds. $$
We now characterize the cone of admissible directions at $g$. Let us introduce the active set $A_g = \{x \in \partial\Omega_N | \|g(x)\| = M\}$. As the pointwise norm $\|g(x)\|$ should not exceed $M$ on $\partial\Omega_N$, admissible directions $h \in L^\infty(\partial\Omega_N)$ are characterized (see for instance [73]) by

$$h(x) \cdot g(x) \le 0 \quad \text{almost everywhere in } A_g.$$

An alternative formulation for the cone of admissible directions is the intersection $L^\infty(\partial\Omega_N) \cap K$, where $K$ is the closed cone

$$K = \left\{ h \in L^2(\partial\Omega_N) \mid \forall \lambda \in L^2(\partial\Omega_N, \mathbb{R}_+) \quad \int_{A_g} \lambda(h \cdot g) \, ds \le 0 \right\}.$$

Let us briefly check this expression. First, it is clear that a direction $h$ which satisfies $h \cdot g \le 0$ almost everywhere in $A_g$ is in $K$. On the other hand, if there is a subset $S \subset A_g$ with positive Lebesgue measure where $h \cdot g > 0$, then, by choosing $\lambda = \mathbf{1}_S$, we obtain $h \notin K$.

We notice that this new expression for $K$ is the expression of a polar cone, namely

$$K = Q^\circ = \left\{ h \in L^2(\partial\Omega_N) \mid \forall q \in Q \quad \langle h, q \rangle \le 0 \right\},$$

where the closed cone $Q$ is defined by

$$Q = \{\mathbf{1}_{A_g} \lambda g, \lambda \in L^2(\partial\Omega_N, \mathbb{R}_+)\}.$$

We denote by $-p_\ell = p_K + p_Q$ the so-called Moreau decomposition of $-p_\ell$ [74] with $p_K \in K$, $p_Q \in Q$ and $\langle p_K, p_Q \rangle = 0$.

Now coming back to the optimality condition, as $g$ is a local minimizer of $\Phi$, the Euler inequation stipulates that $d\Phi(g)(h) \ge 0$ for every element $h$ of the cone of admissible directions. In particular, we consider the admissible direction

$$h = \frac{p_K}{\max(1, \|p_K\|)} \in L^\infty(\partial\Omega_N) \cap K.$$

It results from the Euler inequation that

$$0 \le d\Phi(g)(h) \le \int_{\partial\Omega_N} p_\ell \cdot h \, ds = - \int_{\partial\Omega_N} \min(\|p_K\|, \|p_K\|^2) \, ds.$$

Finally, we obtain $p_K(\cdot) = 0$, which yields the first-order optimality condition $p_\ell + p_Q = 0$. Observing that this condition rewrites as in the statement of Theorem 3.2, we are done.

## References

[1] Haouchine N, Dequidt J, Peterlík I, Kerrien E, Berger M, Cotin S. 2013 Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery. In *2013 IEEE International Symposium on Mixed and Augmented Reality (ISMAR)* pp. 199–208.

[2] Arguillère S, Trélat E, Trouvé A, Younes L. 2015 Shape deformation analysis from the optimal control viewpoint. *Journal de Mathématiques Pures et Appliquées* **104**, 139–178.

[3] Charon N, Charlier B, Trouvé A. 2018 Metamorphoses of functional shapes in Sobolev spaces. *Foundations of Computational Mathematics* **18**, 1535–1596.

[4] Bauer M, Charon N, Younes L. 2019 Chapter 16 - Metric registration of curves and surfaces using optimal control. In Kimmel R, Tai XC, editors, *Processing, Analyzing and Learning of Images, Shapes, and Forms: Part 2* vol. 20*Handbook of Numerical Analysis* pp. 613–646. Elsevier.

[5] Antonsanti PL, Glaunès J, Benseghir T, Jugnon V, Kaltenmark I. 2021 Partial Matching in the Space of Varifolds. In Feragen A, Sommer S, Schnabel J, Nielsen M, editors, *Information Processing in Medical Imaging* pp. 123–135 Cham. Springer International Publishing.
We now characterize the cone of admissible directions at $g$. Let us introduce the active set $A_g = \{x \in \partial\Omega_N | \|g(x)\| = M\}$. As the pointwise norm $\|g(x)\|$ should not exceed $M$ on $\partial\Omega_N$, admissible directions $h \in L^\infty(\partial\Omega_N)$ are characterized (see for instance [73]) by

$$h(x) \cdot g(x) \le 0 \quad \text{almost everywhere in } A_g.$$

An alternative formulation for the cone of admissible directions is the intersection $L^\infty(\partial\Omega_N) \cap K$, where $K$ is the closed cone

$$K = \left\{ h \in L^2(\partial\Omega_N) \mid \forall \lambda \in L^2(\partial\Omega_N, \mathbb{R}_+) \quad \int_{A_g} \lambda(h \cdot g) \, ds \le 0 \right\}.$$

Let us briefly check this expression. First, it is clear that a direction $h$ which satisfies $h \cdot g \le 0$ almost everywhere in $A_g$ is in $K$. On the other hand, if there is a subset $S \subset A_g$ with positive Lebesgue measure where $h \cdot g > 0$, then, by choosing $\lambda = \mathbf{1}_S$, we obtain $h \notin K$.

We notice that this new expression for $K$ is the expression of a polar cone, namely

$$K = Q^\circ = \left\{ h \in L^2(\partial\Omega_N) \mid \forall q \in Q \quad \langle h, q \rangle \le 0 \right\},$$

where the closed cone $Q$ is defined by

$$Q = \{\mathbf{1}_{A_g} \lambda g, \lambda \in L^2(\partial\Omega_N, \mathbb{R}_+)\}.$$

We denote by $-p_\ell = p_K + p_Q$ the so-called Moreau decomposition of $-p_\ell$ [74] with $p_K \in K$, $p_Q \in Q$ and $\langle p_K, p_Q \rangle = 0$.

Now coming back to the optimality condition, as $g$ is a local minimizer of $\Phi$, the Euler inequation stipulates that $d\Phi(g)(h) \ge 0$ for every element $h$ of the cone of admissible directions. In particular, we consider the admissible direction

$$h = \frac{p_K}{\max(1, \|p_K\|)} \in L^\infty(\partial\Omega_N) \cap K.$$

It results from the Euler inequation that

$$0 \le d\Phi(g)(h) \le \int_{\partial\Omega_N} p_\ell \cdot h \, ds = - \int_{\partial\Omega_N} \min(\|p_K\|, \|p_K\|^2) \, ds.$$

Finally, we obtain $p_K(\cdot) = 0$, which yields the first-order optimality condition $p_\ell + p_Q = 0$. Observing that this condition rewrites as in the statement of Theorem 3.2, we are done.

## References

[1] Haouchine N, Dequidt J, Peterlík I, Kerrien E, Berger M, Cotin S. 2013 Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery. In *2013 IEEE International Symposium on Mixed and Augmented Reality (ISMAR)* pp. 199–208.

[2] Arguillère S, Trélat E, Trouvé A, Younes L. 2015 Shape deformation analysis from the optimal control viewpoint. *Journal de Mathématiques Pures et Appliquées* **104**, 139–178.

[3] Charon N, Charlier B, Trouvé A. 2018 Metamorphoses of functional shapes in Sobolev spaces. *Foundations of Computational Mathematics* **18**, 1535–1596.

[4] Bauer M, Charon N, Younes L. 2019 Chapter 16 - Metric registration of curves and surfaces using optimal control. In Kimmel R, Tai XC, editors, *Processing, Analyzing and Learning of Images, Shapes, and Forms: Part 2* vol. 20*Handbook of Numerical Analysis* pp. 613–646. Elsevier.

[5] Antonsanti PL, Glaunès J, Benseghir T, Jugnon V, Kaltenmark I. 2021 Partial Matching in the Space of Varifolds. In Feragen A, Sommer S, Schnabel J, Nielsen M, editors, *Information Processing in Medical Imaging* pp. 123–135 Cham. Springer International Publishing.
[6] Charon N, Younes L. 2023 pp. 1929–1958. In *Shape Spaces: From Geometry to Biological Plausibility*, pp. 1929–1958. Cham: Springer International Publishing.

[7] Marchesseau S, Chatelin S, Delingette H. 2017 Nonlinear Biomechanical Model of the Liver. In Payan Y, Ohayon J, editors, *Biomechanics of Living Organs* vol. 1 *Translational Epigenetics* pp. 243–265. Oxford: Academic Press.

[8] Peterlík I, Duriez C, Cotin S. 2012 Modeling and Real-Time Simulation of a Vascularized Liver Tissue. In Ayache N, Delingette H, Golland P, Mori K, editors, *Medical Image Computing and Computer-Assisted Intervention – MICCAI 2012* pp. 50–57 Berlin, Heidelberg. Springer Berlin Heidelberg.

[9] Haouchine N, Cotin S, Peterlík I, Dequidt J, Lopez MS, Kerrien E, Berger M. 2015 Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery. *IEEE Transactions on Visualization and Computer Graphics* **21**, 584–597.

[10] Ozkan E, Goksel O. 2018 Compliance boundary conditions for patient-specific deformation simulation using the finite element method. *Biomedical Physics & Engineering Express* **4**, 025003.

[11] Nikolaev S, Cotin S. 2020 Estimation of boundary conditions for patient-specific liver simulation during augmented surgery. *International Journal of Computer Assisted Radiology and Surgery* **15**, 1107–1115.

[12] Sotiras A, Davatzikos C, Paragios N. 2013 Deformable medical image registration: A survey. *IEEE transactions on medical imaging* **32**, 1153.

[13] Burger M, Modersitzki J, Ruthotto L. 2013 A hyperelastic regularization energy for image registration. *SIAM Journal on Scientific Computing* **35**, B132–B148.

[14] de Buhan M, Dapogny C, Frey P, Nardoni C. 2016 An optimization method for elastic shape matching. *C. R. Math. Acad. Sci. Paris* **354**, 783–787.

[15] Bauer M, Charon N, Harms P, Hsieh HW. 2021 A numerical framework for elastic surface matching, comparison, and interpolation. *International Journal of Computer Vision* **129**, 2425–2444.

[16] Kurtek S, Jermyn IH, Xie Q, Klassen E, Laga H. 2016 pp. 257–277. In *Elastic Shape Analysis of Surfaces and Images*, pp. 257–277. Cham: Springer International Publishing.

[17] Hartman E, Sukurdeep Y, Klassen E, Charon N, Bauer M. 2023 Elastic Shape Analysis of Surfaces with Second-Order Sobolev Metrics: A Comprehensive Numerical Framework. *International Journal of Computer Vision* **131**, 1183–1209.

[18] Besl PJ, McKay ND. 1992 Method for registration of 3-D shapes. In Schenker PS, editor, *Sensor Fusion IV: Control Paradigms and Data Structures* vol. 1611 pp. 586 – 606. International Society for Optics and Photonics SPIE.

[19] Suwelack S, Röhl S, Bodenstedt S, Reichard D, Dillmann R, dos Santos T, Maier-Hein L, Wagner M, Wünscher J, Kenngott H, Müller BP, Speidel S. 2014 Physics-based shape matching for intraoperative image guidance. *Medical Physics* **41**, 111901.

[20] Plantefève R, Peterlík I, Haouchine N, Cotin S. 2016 Patient-Specific Biomechanical Modeling for Guidance During Minimally-Invasive Hepatic Surgery. *Annals of Biomedical Engineering* **44**, 139–153.

[21] Peterlík I, Courtecuisse H, Rohling R, Abolmaesumi P, Nguan C, Cotin S, Salcudean S. 2018 Fast elastic registration of soft tissues under large deformations. *Medical Image Analysis* **45**, 24–40.

[22] Courtecuisse H, Jiang Z, Mayeur O, Witz JF, Lecomte-Grosbras P, Cosson M, Brieu M, Cotin S. 2020 Three-dimensional physics-based registration of pelvic system using 2D dynamic magnetic resonance imaging slices. *Strain* **56**, e12339.

[23] Morin F, Courtecuisse H, Reinertsen I, Le Lann F, Palombi O, Payan Y, Chabanas M. 2017 Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation. *Medical image analysis* **40**, 133–153.

[24] Rucker DC, Wu Y, Clements LW, Ondrake JE, Pfeiffer TS, Simpson AL, Jarnagin WR, Miga MI. 2014 A Mechanics-Based Nonrigid Registration Method for Liver Surgery Using Sparse Intraoperative Data. *IEEE Transactions on Medical Imaging* **33**, 147–158.

[25] Heiselman JS, Clements LW, Collins JA, Weis JA, Simpson AL, Geevarghese SK, Kingham TP, Jarnagin WR, Miga MI. 2017 Characterization and correction of intraoperative soft tissue deformation in image-guided laparoscopic liver surgery. *Journal of Medical Imaging* **5**, 1–12.

[26] Heiselman JS, Jarnagin WR, Miga MI. 2020 Intraoperative Correction of Liver Deformation Using Sparse Surface and Vascular Features via Linearized Iterative Boundary Reconstruction. *IEEE Transactions on Medical Imaging* **39**, 2223–2234.

[27] Özgür E, Koo B, Le Roy B, Buc E, Bartoli A. 2018 Preoperative liver registration for augmented monocular laparoscopy using backward-forward biomechanical simulation. *International Journal of Computer Assisted Radiology and Surgery* **13**, 1629–1640.

[28] Xie H, Song J, Zhong Y, Gu C, Choi KS. 2022 Constrained finite element method for runtime modeling of soft tissue deformation. *Applied Mathematical Modelling* **109**, 599–612.
[6] Charon N, Younes L. 2023 pp. 1929–1958. In *Shape Spaces: From Geometry to Biological Plausibility*, pp. 1929–1958. Cham: Springer International Publishing.

[7] Marchesseau S, Chatelin S, Delingette H. 2017 Nonlinear Biomechanical Model of the Liver. In Payan Y, Ohayon J, editors, *Biomechanics of Living Organs* vol. 1 *Translational Epigenetics* pp. 243–265. Oxford: Academic Press.

[8] Peterlík I, Duriez C, Cotin S. 2012 Modeling and Real-Time Simulation of a Vascularized Liver Tissue. In Ayache N, Delingette H, Golland P, Mori K, editors, *Medical Image Computing and Computer-Assisted Intervention – MICCAI 2012* pp. 50–57 Berlin, Heidelberg. Springer Berlin Heidelberg.

[9] Haouchine N, Cotin S, Peterlík I, Dequidt J, Lopez MS, Kerrien E, Berger M. 2015 Impact of Soft Tissue Heterogeneity on Augmented Reality for Liver Surgery. *IEEE Transactions on Visualization and Computer Graphics* **21**, 584–597.

[10] Ozkan E, Goksel O. 2018 Compliance boundary conditions for patient-specific deformation simulation using the finite element method. *Biomedical Physics & Engineering Express* **4**, 025003.

[11] Nikolaev S, Cotin S. 2020 Estimation of boundary conditions for patient-specific liver simulation during augmented surgery. *International Journal of Computer Assisted Radiology and Surgery* **15**, 1107–1115.

[12] Sotiras A, Davatzikos C, Paragios N. 2013 Deformable medical image registration: A survey. *IEEE transactions on medical imaging* **32**, 1153.

[13] Burger M, Modersitzki J, Ruthotto L. 2013 A hyperelastic regularization energy for image registration. *SIAM Journal on Scientific Computing* **35**, B132–B148.

[14] de Buhan M, Dapogny C, Frey P, Nardoni C. 2016 An optimization method for elastic shape matching. *C. R. Math. Acad. Sci. Paris* **354**, 783–787.

[15] Bauer M, Charon N, Harms P, Hsieh HW. 2021 A numerical framework for elastic surface matching, comparison, and interpolation. *International Journal of Computer Vision* **129**, 2425–2444.

[16] Kurtek S, Jermyn IH, Xie Q, Klassen E, Laga H. 2016 pp. 257–277. In *Elastic Shape Analysis of Surfaces and Images*, pp. 257–277. Cham: Springer International Publishing.

[17] Hartman E, Sukurdeep Y, Klassen E, Charon N, Bauer M. 2023 Elastic Shape Analysis of Surfaces with Second-Order Sobolev Metrics: A Comprehensive Numerical Framework. *International Journal of Computer Vision* **131**, 1183–1209.

[18] Besl PJ, McKay ND. 1992 Method for registration of 3-D shapes. In Schenker PS, editor, *Sensor Fusion IV: Control Paradigms and Data Structures* vol. 1611 pp. 586 – 606. International Society for Optics and Photonics SPIE.

[19] Suwelack S, Röhl S, Bodenstedt S, Reichard D, Dillmann R, dos Santos T, Maier-Hein L, Wagner M, Wünscher J, Kenngott H, Müller BP, Speidel S. 2014 Physics-based shape matching for intraoperative image guidance. *Medical Physics* **41**, 111901.

[20] Plantefève R, Peterlík I, Haouchine N, Cotin S. 2016 Patient-Specific Biomechanical Modeling for Guidance During Minimally-Invasive Hepatic Surgery. *Annals of Biomedical Engineering* **44**, 139–153.

[21] Peterlík I, Courtecuisse H, Rohling R, Abolmaesumi P, Nguan C, Cotin S, Salcudean S. 2018 Fast elastic registration of soft tissues under large deformations. *Medical Image Analysis* **45**, 24–40.

[22] Courtecuisse H, Jiang Z, Mayeur O, Witz JF, Lecomte-Grosbras P, Cosson M, Brieu M, Cotin S. 2020 Three-dimensional physics-based registration of pelvic system using 2D dynamic magnetic resonance imaging slices. *Strain* **56**, e12339.

[23] Morin F, Courtecuisse H, Reinertsen I, Le Lann F, Palombi O, Payan Y, Chabanas M. 2017 Brain-shift compensation using intraoperative ultrasound and constraint-based biomechanical simulation. *Medical image analysis* **40**, 133–153.

[24] Rucker DC, Wu Y, Clements LW, Ondrake JE, Pfeiffer TS, Simpson AL, Jarnagin WR, Miga MI. 2014 A Mechanics-Based Nonrigid Registration Method for Liver Surgery Using Sparse Intraoperative Data. *IEEE Transactions on Medical Imaging* **33**, 147–158.

[25] Heiselman JS, Clements LW, Collins JA, Weis JA, Simpson AL, Geevarghese SK, Kingham TP, Jarnagin WR, Miga MI. 2017 Characterization and correction of intraoperative soft tissue deformation in image-guided laparoscopic liver surgery. *Journal of Medical Imaging* **5**, 1–12.

[26] Heiselman JS, Jarnagin WR, Miga MI. 2020 Intraoperative Correction of Liver Deformation Using Sparse Surface and Vascular Features via Linearized Iterative Boundary Reconstruction. *IEEE Transactions on Medical Imaging* **39**, 2223–2234.

[27] Özgür E, Koo B, Le Roy B, Buc E, Bartoli A. 2018 Preoperative liver registration for augmented monocular laparoscopy using backward-forward biomechanical simulation. *International Journal of Computer Assisted Radiology and Surgery* **13**, 1629–1640.

[28] Xie H, Song J, Zhong Y, Gu C, Choi KS. 2022 Constrained finite element method for runtime modeling of soft tissue deformation. *Applied Mathematical Modelling* **109**, 599–612.
[29] Mestdagh G, Cotin S. 2022 An Optimal Control Problem for Elastic Registration and Force Estimation in Augmented Surgery. In Wang L, Dou Q, Fletcher PT, Speidel S, Li S, editors, *Medical Image Computing and Computer Assisted Intervention – MICCAI 2022* pp. 74–83 Cham. Springer Nature Switzerland.

[30] Kaltenmark I, Trouvé A. 2019 Estimation of a growth development with partial diffeomorphic mappings. *Quarterly of Applied Mathematics* **77**, 227–267.

[31] Hsieh DN, Arguillère S, Charon N, Younes L. 2022 Mechanistic modeling of longitudinal shape changes: equations of motion and inverse problems. *SIAM J. Appl. Dyn. Syst.* **21**, 80–101.

[32] Hogea C, Davatzikos C, Biros G. 2008 An image-driven parameter estimation problem for a reaction–diffusion glioma growth model with mass effects. *Journal of Mathematical Biology* **56**, 793–825.

[33] Gooya A, Pohl KM, Bilello M, Cirillo L, Biros G, Melhem ER, Davatzikos C. 2012 GLISTR: Glioma Image Segmentation and Registration. *IEEE Transactions on Medical Imaging* **31**, 1941–1954.

[34] Scheufele K, Mang A, Gholami A, Davatzikos C, Biros G, Mehl M. 2019 Coupling brain-tumor biophysical models and diffeomorphic image registration. *Computer Methods in Applied Mechanics and Engineering* **347**, 533–567.

[35] Sundar H, Davatzikos C, Biros G. 2009 Biomechanically-Constrained 4D Estimation of Myocardial Motion. In Yang GZ, Hawkes D, Rueckert D, Noble A, Taylor C, editors, *Medical Image Computing and Computer-Assisted Intervention – MICCAI 2009* pp. 257–265 Berlin, Heidelberg. Springer Berlin Heidelberg.

[36] Ortigosa R, Martínez-Frutos J, Mora-Corral C, Pedregal P, Periago F. 2020 Optimal Control of Soft Materials Using a Hausdorff Distance Functional. *SIAM Journal on Control and Optimization*.

[37] Miller K, Joldes G, Lance D, Wittek A. 2007 Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation. *Communications in Numerical Methods in Engineering* **23**, 121–134.

[38] Nesme M, Payan Y, Faure F. 2005 Efficient, Physically Plausible Finite Elements. In *Eurographics* Short papers Dublin, Ireland.

[39] van Kaick O, Zhang H, Hamarneh G, Cohen-Or D. 2011 A Survey on Shape Correspondence. *Computer Graphics Forum* **30**, 1681–1707.

[40] Sahillioğlu Y. 2020 Recent advances in shape correspondence. *The Visual Computer* **36**, 1705–1721.

[41] Kaltenmark I, Charlier B, Charon N. 2017 A General Framework for Curve and Surface Comparison and Registration With Oriented Varifolds. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*.

[42] Rodolá E, Cosmo L, Bronstein MM, Torsello A, Cremers D. 2017 Partial Functional Correspondence. *Computer Graphics Forum* **36**, 222–236.

[43] Gröger K. 1989 A $W^{1,p}$-estimate for solutions to mixed boundary value problems for second order elliptic differential equations. *Math. Ann.* **283**, 679–687.

[44] Haller-Dintelmann R, Meyer C, Rehberg J, Schiela A. 2009 Hölder continuity and optimal control for nonsmooth elliptic problems. *Appl. Math. Optim.* **60**, 397–428.

[45] Gröger K, Rehberg J. 1989 Resolvent estimates in $W^{-1,p}$ for second order elliptic differential operators in case of mixed boundary conditions. *Math. Ann.* **285**, 105–113.

[46] Droniou J. 2000 Solving convection-diffusion equations with mixed, Neumann and Fourier boundary conditions and measures as data, by a duality method. *Adv. Differential Equations* **5**, 1341–1396.

[47] Shi P, Wright S. 1994 Higher integrability of the gradient in linear elasticity. *Math. Ann.* **299**, 435–448.

[48] Herzog R, Meyer C, Wachsmuth G. 2011 Integrability of displacement and stresses in linear and nonlinear elasticity with mixed boundary conditions. *J. Math. Anal. Appl.* **382**, 802–813.

[49] Brezis H. 1983 *Analyse fonctionnelle*. Collection Mathématiques Appliquées pour la Maîtrise. [Collection of Applied Mathematics for the Master’s Degree]. Masson, Paris. Théorie et applications. [Theory and applications].

[50] Krumbiegel K, Rehberg J. 2013 Second Order Sufficient Optimality Conditions for Parabolic Optimal Control Problems with Pointwise State Constraints. *SIAM Journal on Control and Optimization* **51**, 304–331.

[51] Disser K, Kaiser HC, Rehberg J. 2015 Optimal Sobolev Regularity for Linear Second-Order Divergence Elliptic Operators Occurring in Real-World Problems. *SIAM Journal on Mathematical Analysis* **47**, 1719–1746.

[52] Mestdagh G. 2022 *An optimal control formulation for organ registration in augmented surgery*. Theses Université de Strasbourg.

[53] Danskin JM. 1967 *The theory of max-min and its application to weapons allocation problems*. Econometrics and Operations Research, Vol. V. Springer-Verlag New York, Inc., New York.

[54] Ponce AC. 2016 *Elliptic PDEs, measures and capacities* vol. 23EMS Tracts in Mathematics. European Mathematical Society (EMS), Zürich. From the Poisson equations to nonlinear Thomas-Fermi problems.
[29] Mestdagh G, Cotin S. 2022 An Optimal Control Problem for Elastic Registration and Force Estimation in Augmented Surgery. In Wang L, Dou Q, Fletcher PT, Speidel S, Li S, editors, *Medical Image Computing and Computer Assisted Intervention – MICCAI 2022* pp. 74–83 Cham. Springer Nature Switzerland.

[30] Kaltenmark I, Trouvé A. 2019 Estimation of a growth development with partial diffeomorphic mappings. *Quarterly of Applied Mathematics* **77**, 227–267.

[31] Hsieh DN, Arguillère S, Charon N, Younes L. 2022 Mechanistic modeling of longitudinal shape changes: equations of motion and inverse problems. *SIAM J. Appl. Dyn. Syst.* **21**, 80–101.

[32] Hogea C, Davatzikos C, Biros G. 2008 An image-driven parameter estimation problem for a reaction–diffusion glioma growth model with mass effects. *Journal of Mathematical Biology* **56**, 793–825.

[33] Gooya A, Pohl KM, Bilello M, Cirillo L, Biros G, Melhem ER, Davatzikos C. 2012 GLISTR: Glioma Image Segmentation and Registration. *IEEE Transactions on Medical Imaging* **31**, 1941–1954.

[34] Scheufele K, Mang A, Gholami A, Davatzikos C, Biros G, Mehl M. 2019 Coupling brain-tumor biophysical models and diffeomorphic image registration. *Computer Methods in Applied Mechanics and Engineering* **347**, 533–567.

[35] Sundar H, Davatzikos C, Biros G. 2009 Biomechanically-Constrained 4D Estimation of Myocardial Motion. In Yang GZ, Hawkes D, Rueckert D, Noble A, Taylor C, editors, *Medical Image Computing and Computer-Assisted Intervention – MICCAI 2009* pp. 257–265 Berlin, Heidelberg. Springer Berlin Heidelberg.

[36] Ortigosa R, Martínez-Frutos J, Mora-Corral C, Pedregal P, Periago F. 2020 Optimal Control of Soft Materials Using a Hausdorff Distance Functional. *SIAM Journal on Control and Optimization*.

[37] Miller K, Joldes G, Lance D, Wittek A. 2007 Total Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation. *Communications in Numerical Methods in Engineering* **23**, 121–134.

[38] Nesme M, Payan Y, Faure F. 2005 Efficient, Physically Plausible Finite Elements. In *Eurographics* Short papers Dublin, Ireland.

[39] van Kaick O, Zhang H, Hamarneh G, Cohen-Or D. 2011 A Survey on Shape Correspondence. *Computer Graphics Forum* **30**, 1681–1707.

[40] Sahillioğlu Y. 2020 Recent advances in shape correspondence. *The Visual Computer* **36**, 1705–1721.

[41] Kaltenmark I, Charlier B, Charon N. 2017 A General Framework for Curve and Surface Comparison and Registration With Oriented Varifolds. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*.

[42] Rodolá E, Cosmo L, Bronstein MM, Torsello A, Cremers D. 2017 Partial Functional Correspondence. *Computer Graphics Forum* **36**, 222–236.

[43] Gröger K. 1989 A $W^{1,p}$-estimate for solutions to mixed boundary value problems for second order elliptic differential equations. *Math. Ann.* **283**, 679–687.

[44] Haller-Dintelmann R, Meyer C, Rehberg J, Schiela A. 2009 Hölder continuity and optimal control for nonsmooth elliptic problems. *Appl. Math. Optim.* **60**, 397–428.

[45] Gröger K, Rehberg J. 1989 Resolvent estimates in $W^{-1,p}$ for second order elliptic differential operators in case of mixed boundary conditions. *Math. Ann.* **285**, 105–113.

[46] Droniou J. 2000 Solving convection-diffusion equations with mixed, Neumann and Fourier boundary conditions and measures as data, by a duality method. *Adv. Differential Equations* **5**, 1341–1396.

[47] Shi P, Wright S. 1994 Higher integrability of the gradient in linear elasticity. *Math. Ann.* **299**, 435–448.

[48] Herzog R, Meyer C, Wachsmuth G. 2011 Integrability of displacement and stresses in linear and nonlinear elasticity with mixed boundary conditions. *J. Math. Anal. Appl.* **382**, 802–813.

[49] Brezis H. 1983 *Analyse fonctionnelle*. Collection Mathématiques Appliquées pour la Maîtrise. [Collection of Applied Mathematics for the Master’s Degree]. Masson, Paris. Théorie et applications. [Theory and applications].

[50] Krumbiegel K, Rehberg J. 2013 Second Order Sufficient Optimality Conditions for Parabolic Optimal Control Problems with Pointwise State Constraints. *SIAM Journal on Control and Optimization* **51**, 304–331.

[51] Disser K, Kaiser HC, Rehberg J. 2015 Optimal Sobolev Regularity for Linear Second-Order Divergence Elliptic Operators Occurring in Real-World Problems. *SIAM Journal on Mathematical Analysis* **47**, 1719–1746.

[52] Mestdagh G. 2022 *An optimal control formulation for organ registration in augmented surgery*. Theses Université de Strasbourg.

[53] Danskin JM. 1967 *The theory of max-min and its application to weapons allocation problems*. Econometrics and Operations Research, Vol. V. Springer-Verlag New York, Inc., New York.

[54] Ponce AC. 2016 *Elliptic PDEs, measures and capacities* vol. 23EMS Tracts in Mathematics. European Mathematical Society (EMS), Zürich. From the Poisson equations to nonlinear Thomas-Fermi problems.
[55] Littman W, Stampacchia G, Weinberger HF. 1963 Regular points for elliptic equations with discontinuous coefficients. *Ann. Scuola Norm. Sup. Pisa Cl. Sci.* (3) **17**, 43–77.

[56] Guttman A. 1984 R-Trees: A Dynamic Index Structure for Spatial Searching. *SIGMOD Rec.* **14**, 47–57.

[57] Allard J, Cotin S, Faure F, Bensoussan PJ, Poyer F, Duriez C, Delingette H, Grisoni L. 2007 SOFA - an Open Source Framework for Medical Simulation. In *MMVR 15 - Medicine Meets Virtual Reality* vol. 125* *Studies in Health Technology and Informatics* pp. 13–18 Palm Beach, United States. IOP Press.

[58] Byrd RH, Lu P, Nocedal J, Zhu C. 1995 A Limited Memory Algorithm for Bound Constrained Optimization. *SIAM Journal on Scientific Computing* **16**, 1190–1208.

[59] Collins JA, Weis JA, Heiselman JS, Clements LW, Simpson AL, Jarnagin WR, Miga MI. 2017 Improving Registration Robustness for Image-Guided Liver Surgery in a Novel Human-to-Phantom Data Framework. *IEEE Transactions on Medical Imaging* **36**, 1502–1510.

[60] Brewer EL, Clements LW, Collins JA, Doss DJ, Heiselman JS, Miga MI, Pavas CD, III EHW. 2019 The image-to-physical liver registration sparse data challenge. In Fei B, Linte CA, editors, *Medical Imaging 2019: Image-Guided Procedures, Robotic Interventions, and Modeling* vol. 10951 pp. 364–370. International Society for Optics and Photonics SPIE.

[61] Jia M, Kyan M. 2021 Improving Intraoperative Liver Registration in Image-Guided Surgery with Learning-Based Reconstruction. In *ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* pp. 1230–1234.

[62] Pfeiffer M, Riediger C, Leger S, Kühn JP, Seppelt D, Hoffmann RT, Weitz J, Speidel S. 2020 Non-Rigid Volume to Surface Registration Using a Data-Driven Biomechanical Model. In Martel AL, Abolmaesumi P, Stoyanov D, Mateus D, Zuluaga MA, Zhou SK, Racoceanu D, Joskowicz L, editors, *Medical Image Computing and Computer Assisted Intervention – MICCAI 2020* pp. 724–734 Cham. Springer International Publishing.

[63] Heiselman JS, Collins JA, Ringel MJ, Jarnagin WR, Miga MI. 2023 Comparison study of sparse data-driven soft tissue registration: preliminary results from the image-to-physical liver registration sparse data challenge. In Linte CA, Siewerdsen JH, editors, *Medical Imaging 2023: Image-Guided Procedures, Robotic Interventions, and Modeling* vol. 12466 p. 124660M. International Society for Optics and Photonics SPIE.

[64] Bonet J, Wood RD. 2008 *Nonlinear continuum mechanics for finite element analysis*. Cambridge University Press, Cambridge second edition.

[65] Gross C, Krause R. 2009 On the Convergence of Recursive Trust-Region Methods for Multiscale Nonlinear Optimization and Applications to Nonlinear Mechanics. *SIAM Journal on Numerical Analysis* **47**, 3044–3069.

[66] Youett J, Sander O, Kornhuber R. 2019 A Globally Convergent Filter-Trust-Region Method for Large Deformation Contact Problems. *SIAM Journal on Scientific Computing* **41**, B114–B138.

[67] Gelin JC, Picart P. 1988 Use of quasi-Newton methods for large strain elastic-plastic finite element computations. *Communications in Applied Numerical Methods* **4**, 457–469.

[68] Yusa Y, Miyauchi S, Okada H. 2021 Performance investigation of quasi-Newton-based parallel nonlinear FEM for large-deformation elastic-plastic analysis over 100 thousand degrees of freedom. *Mechanical Engineering Journal* **8**, 21–00053–21–00053.

[69] Odot A, Mestdagh G, Privat Y, Cotin S. 2023 Real-time elastic partial shape matching using a neural network-based adjoint method. Preprint: https://hal.inria.fr/hal-04019777v1.

[70] Ciarlet PG. 1988 *Mathematical elasticity. Vol. I* vol. 20*Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam. Three-dimensional elasticity.

[71] Evans LC. 2010 *Partial differential equations* vol. 19*Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI second edition.

[72] Adams RA, Fournier JJF. 2003 *Sobolev spaces* vol. 140*Pure and Applied Mathematics (Amsterdam). Elsevier/Academic Press, Amsterdam second edition.

[73] Bednarczuk E, Pierre M, Rouy E, Sokolowski J. 2000 Tangent Sets in Some Functional Spaces. *Nonlinear Anal.* **42**, 871–886.

[74] Moreau JJ. 1962 Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires. *Comptes rendus hebdomadaires des séances de l'Académie des sciences* **255**, 238–240.
[55] Littman W, Stampacchia G, Weinberger HF. 1963 Regular points for elliptic equations with discontinuous coefficients. *Ann. Scuola Norm. Sup. Pisa Cl. Sci.* (3) **17**, 43–77.

[56] Guttman A. 1984 R-Trees: A Dynamic Index Structure for Spatial Searching. *SIGMOD Rec.* **14**, 47–57.

[57] Allard J, Cotin S, Faure F, Bensoussan PJ, Poyer F, Duriez C, Delingette H, Grisoni L. 2007 SOFA - an Open Source Framework for Medical Simulation. In *MMVR 15 - Medicine Meets Virtual Reality* vol. 125* *Studies in Health Technology and Informatics* pp. 13–18 Palm Beach, United States. IOP Press.

[58] Byrd RH, Lu P, Nocedal J, Zhu C. 1995 A Limited Memory Algorithm for Bound Constrained Optimization. *SIAM Journal on Scientific Computing* **16**, 1190–1208.

[59] Collins JA, Weis JA, Heiselman JS, Clements LW, Simpson AL, Jarnagin WR, Miga MI. 2017 Improving Registration Robustness for Image-Guided Liver Surgery in a Novel Human-to-Phantom Data Framework. *IEEE Transactions on Medical Imaging* **36**, 1502–1510.

[60] Brewer EL, Clements LW, Collins JA, Doss DJ, Heiselman JS, Miga MI, Pavas CD, III EHW. 2019 The image-to-physical liver registration sparse data challenge. In Fei B, Linte CA, editors, *Medical Imaging 2019: Image-Guided Procedures, Robotic Interventions, and Modeling* vol. 10951 pp. 364–370. International Society for Optics and Photonics SPIE.

[61] Jia M, Kyan M. 2021 Improving Intraoperative Liver Registration in Image-Guided Surgery with Learning-Based Reconstruction. In *ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)* pp. 1230–1234.

[62] Pfeiffer M, Riediger C, Leger S, Kühn JP, Seppelt D, Hoffmann RT, Weitz J, Speidel S. 2020 Non-Rigid Volume to Surface Registration Using a Data-Driven Biomechanical Model. In Martel AL, Abolmaesumi P, Stoyanov D, Mateus D, Zuluaga MA, Zhou SK, Racoceanu D, Joskowicz L, editors, *Medical Image Computing and Computer Assisted Intervention – MICCAI 2020* pp. 724–734 Cham. Springer International Publishing.

[63] Heiselman JS, Collins JA, Ringel MJ, Jarnagin WR, Miga MI. 2023 Comparison study of sparse data-driven soft tissue registration: preliminary results from the image-to-physical liver registration sparse data challenge. In Linte CA, Siewerdsen JH, editors, *Medical Imaging 2023: Image-Guided Procedures, Robotic Interventions, and Modeling* vol. 12466 p. 124660M. International Society for Optics and Photonics SPIE.

[64] Bonet J, Wood RD. 2008 *Nonlinear continuum mechanics for finite element analysis*. Cambridge University Press, Cambridge second edition.

[65] Gross C, Krause R. 2009 On the Convergence of Recursive Trust-Region Methods for Multiscale Nonlinear Optimization and Applications to Nonlinear Mechanics. *SIAM Journal on Numerical Analysis* **47**, 3044–3069.

[66] Youett J, Sander O, Kornhuber R. 2019 A Globally Convergent Filter-Trust-Region Method for Large Deformation Contact Problems. *SIAM Journal on Scientific Computing* **41**, B114–B138.

[67] Gelin JC, Picart P. 1988 Use of quasi-Newton methods for large strain elastic-plastic finite element computations. *Communications in Applied Numerical Methods* **4**, 457–469.

[68] Yusa Y, Miyauchi S, Okada H. 2021 Performance investigation of quasi-Newton-based parallel nonlinear FEM for large-deformation elastic-plastic analysis over 100 thousand degrees of freedom. *Mechanical Engineering Journal* **8**, 21–00053–21–00053.

[69] Odot A, Mestdagh G, Privat Y, Cotin S. 2023 Real-time elastic partial shape matching using a neural network-based adjoint method. Preprint: https://hal.inria.fr/hal-04019777v1.

[70] Ciarlet PG. 1988 *Mathematical elasticity. Vol. I* vol. 20*Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam. Three-dimensional elasticity.

[71] Evans LC. 2010 *Partial differential equations* vol. 19*Graduate Studies in Mathematics*. American Mathematical Society, Providence, RI second edition.

[72] Adams RA, Fournier JJF. 2003 *Sobolev spaces* vol. 140*Pure and Applied Mathematics (Amsterdam). Elsevier/Academic Press, Amsterdam second edition.

[73] Bednarczuk E, Pierre M, Rouy E, Sokolowski J. 2000 Tangent Sets in Some Functional Spaces. *Nonlinear Anal.* **42**, 871–886.

[74] Moreau JJ. 1962 Décomposition orthogonale d'un espace hilbertien selon deux cônes mutuellement polaires. *Comptes rendus hebdomadaires des séances de l'Académie des sciences* **255**, 238–240.