# Finite Strain Formulation of the Discrete Equilibrium Gap Principle: Application to Mechanically Consistent Regularization for Large Motion Tracking

Martin Genet

► To cite this version:

Martin Genet. Finite Strain Formulation of the Discrete Equilibrium Gap Principle: Application to Mechanically Consistent Regularization for Large Motion Tracking. Comptes Rendus. Mécanique, 2023, 351, pp.429-458. 10.5802/crmeca.228. hal-04132311v4

HAL Id: hal-04132311

https://hal.science/hal-04132311v4

Submitted on 19 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
# Comptes Rendus
## Mécanique

Martin Genet

Finite strain formulation of the discrete equilibrium gap principle: application to
mechanically consistent regularization for large motion tracking

Volume 351 (2023), p. 429-458

Published online: 8 December 2023

https://doi.org/10.5802/crmeca.228

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
# Comptes Rendus
## Mécanique

Martin Genet

Finite strain formulation of the discrete equilibrium gap principle: application to
mechanically consistent regularization for large motion tracking

Volume 351 (2023), p. 429-458

Published online: 8 December 2023

https://doi.org/10.5802/crmeca.228

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# Finite strain formulation of the discrete
equilibrium gap principle: application to
mechanically consistent regularization for
large motion tracking

*Formulation en transformation finie du principe d'écart
d'équilibre discret : application à la régularisation
mécanique pour le suivi de mouvement*

Martin Genet<sup>a, b</sup>

**Abstract.** The equilibrium gap principle offers a good trade-off between robustness and accuracy for regularizing motion tracking, as it simply enforces that the tracked motion corresponds to a body deforming under arbitrary loadings. This paper introduces an extension of the equilibrium gap principle in the large deformation setting, a novel regularization term to control surface tractions, both in the context of finite element motion tracking, and an inverse problem consistent reformulation of the tracking problem. Tracking performance of the proposed method, with displacement resolution down to the pixel size, is demonstrated on synthetic images representing various motions with various signal-to-noise ratios.

**Résumé.** Le principe de l'écart d'équilibre offre un bon compromis entre la robustesse et la précision pour la régularisation du suivi du mouvement, car il impose simplement que le mouvement suivi corresponde à celui d'un corps se déformant sous l'effet de charges arbitraires. Cet article présente une extension du principe de l'écart d'équilibre dans le cadre des grandes déformations, un nouveau terme de régularisation pour contrôler les tractions de surface, les deux dans le contexte du suivi de mouvement par éléments finis, et une reformulation « problème inverse » cohérente du problème de suivi de mouvement avec régularisation mécanique. Les performances de suivi de la méthode proposée, avec une résolution de déplacement allant jusqu'à la taille du pixel de l'image, sont démontrées sur des images synthétiques représentant divers mouvements avec différents rapports signal-sur-bruit.

**Keywords.** Motion tracking, Mechanical regularization, Equilibrium gap principle, Finite element method, Inverse problems.

**Mots-clés.** Suivi de mouvement, Régularisation mécanique, Principe de l'écart d'équilibre, Méthode des éléments finis, Problèmes inverses.

Manuscript received 7 June 2023, revised 22 August 2023, accepted 26 September 2023.

<footnote><sup>a</sup> Laboratoire de Mécanique des Solides, École Polytechnique/IPP/CNRS, Palaiseau, France</footnote>

<footnote><sup>b</sup> Équipe MΞDISIM, INRIA, Palaiseau, France</footnote>

<footnote>E-mail: martin.genet@polytechnique.edu</footnote>
# 1. Introduction

Motion tracking is an important field of image processing, with many application domains from experimental mechanics to biomedical engineering. In experimental mechanics, especially in the content of material parameter identification, it induced a true change of paradigm, as it is no longer necessary to perform delicate experiments with simple kinematics like pure tension or compression; instead rather complex experiments can now be performed, involving potentially many deformation mechanisms, as long as the potentially complex kinematics can be tracked based on surface or volume images [1-4]. In biomedical engineering, it allows for the quantitative analysis of biomedical images, hence to derive objective and quantitative biomarkers for improved diagnosis, either directly based on kinematics [5, 6], or by merging physical models and imaging data into so called digital twins [7, 8].

Many approaches have been developed over the past decades, based on many variants of the many aspects of the method, such as harmonic/Fourier versus intensity/features tracking, local versus global approaches, etc. [9-12]. In this paper, we use an intensity-based global tracking approach, which is the most natural to integrate our novel regularization approach, though other tracking approaches could have been considered as well. Existing intensity-based global approaches differ notably in their motion models (splines, finite elements, etc.), image similarity metrics (mean squared error, structural similarity index, mutual information, etc.), optimization methods (gradient descent, Gauss-Newton, Newton, etc.), etc. [11, 12].

One key question in motion tracking is the regularization, i.e., the *a priori* knowledge introduced in the process to improve the quality of the tracking. It is required by the intrinsic ill-posedness of the problem (we are looking for a vector field—the displacement—from an input scalar field—the image—), as well as image finite resolution, noise and bias. As in any inverse or optimization problem, efficient regularization requires a fine trade-off, here between providing enough constraint to help the tracking quality and robustness, while providing enough freedom so as to not interfere with the actual motion. In some sense regularization allows for a control, through penalization, of the function space into which the solution is sought. Many regularization terms have been proposed in the literature, such as Laplacian smoothing [13], fluid-like mechanical regularization [14], incompressibility [15], hyperelastic energy [16], etc., some of which will be discussed in details in this paper. An optimal trade-off is arguably reached by the so-called equilibrium gap regularization [2, 17], which puts no direct constraint on the kinematics, while enforcing that the motion is close to a solution of a mechanics problem, in a sense that will be specified later in the paper.

The equilibrium gap principle was originally formulated, at the discrete level and in the linear setting, in the context of material parameter identification based on full-field measurement, in [18], and was later used in the context of motion tracking notably in [19]. An extension to the non linear setting was proposed, at the continuous level, in [17,20,21], with multiple applications to biomedical images [22-24]. In this paper, we propose another extension, still in the nonlinear finite strain setting but at the discrete level—hence allowing to better distinguish the equilibrium gaps induced by the motion itself and the finite element discretization, as will be detailed in the paper—, and show that it performs better than all previous formulations.

As already mentioned, and discussed in details later in the paper, the equilibrium gap regularization consists in enforcing that the obtained displacement is close to a solution of a mechanics problem with generic material behavior and arbitrary imposed surface tractions [17, 19]. To actually obtain some regularization of the problem, this arbitrariness must be handled, and the surface tractions must be somehow controlled independently of the discretization of the displacement field [19]. A surface Laplacian of the displacement was used as an additional regularization term in [19], which unfortunately does not generalize to the large motion setting. One
# 1. Introduction

Motion tracking is an important field of image processing, with many application domains from experimental mechanics to biomedical engineering. In experimental mechanics, especially in the content of material parameter identification, it induced a true change of paradigm, as it is no longer necessary to perform delicate experiments with simple kinematics like pure tension or compression; instead rather complex experiments can now be performed, involving potentially many deformation mechanisms, as long as the potentially complex kinematics can be tracked based on surface or volume images [1–4]. In biomedical engineering, it allows for the quantitative analysis of biomedical images, hence to derive objective and quantitative biomarkers for improved diagnosis, either directly based on kinematics [5, 6], or by merging physical models and imaging data into so called digital twins [7, 8].

Many approaches have been developed over the past decades, based on many variants of the many aspects of the method, such as harmonic/Fourier versus intensity/features tracking, local versus global approaches, etc. [9–12]. In this paper, we use an intensity-based global tracking approach, which is the most natural to integrate our novel regularization approach, though other tracking approaches could have been considered as well. Existing intensity-based global approaches differ notably in their motion models (splines, finite elements, etc.), image similarity metrics (mean squared error, structural similarity index, mutual information, etc.), optimization methods (gradient descent, Gauss–Newton, Newton, etc.), etc. [11, 12].

One key question in motion tracking is the regularization, i.e., the *a priori* knowledge introduced in the process to improve the quality of the tracking. It is required by the intrinsic ill-posedness of the problem (we are looking for a vector field—the displacement—from an input scalar field—the image—), as well as image finite resolution, noise and bias. As in any inverse or optimization problem, efficient regularization requires a fine trade-off, here between providing enough constraint to help the tracking quality and robustness, while providing enough freedom so as to not interfere with the actual motion. In some sense regularization allows for a control, through penalization, of the function space into which the solution is sought. Many regularization terms have been proposed in the literature, such as Laplacian smoothing [13], fluid-like mechanical regularization [14], incompressibility [15], hyperelastic energy [16], etc., some of which will be discussed in details in this paper. An optimal trade-off is arguably reached by the so-called equilibrium gap regularization [2, 17], which puts no direct constraint on the kinematics, while enforcing that the motion is close to a solution of a mechanics problem, in a sense that will be specified later in the paper.

The equilibrium gap principle was originally formulated, at the discrete level and in the linear setting, in the context of material parameter identification based on full-field measurement, in [18], and was later used in the context of motion tracking notably in [19]. An extension to the non linear setting was proposed, at the continuous level, in [17, 20, 21], with multiple applications to biomedical images [22–24]. In this paper, we propose another extension, still in the nonlinear finite strain setting but at the discrete level—hence allowing to better distinguish the equilibrium gaps induced by the motion itself and the finite element discretization, as will be detailed in the paper—, and show that it performs better than all previous formulations.

As already mentioned, and discussed in details later in the paper, the equilibrium gap regularization consists in enforcing that the obtained displacement is close to a solution of a mechanics problem with generic material behavior and arbitrary imposed surface tractions [17, 19]. To actually obtain some regularization of the problem, this arbitrariness must be handled, and the surface tractions must be somehow controlled independently of the discretization of the displacement field [19]. A surface Laplacian of the displacement was used as an additional regularization term in [19], which unfortunately does not generalize to the large motion setting. One
option would be to use a separate discretization for the displacement and the tractions, which however would represent a significant technical difficulty. Instead, in this paper, we propose an additional regularization term based on the surface gradient of the normal and tangential components of the surface tractions, and show that it performs as expected.

The rest of the paper is organized as follows. We first formulate the general motion tracking problem (Section 2.1), then we provide a short literature review on mechanical regularization (Section 2.2), and we describe our proposed regularization term, including body (Section 2.3.1) and boundary (Section 2.3.2) terms. Then, we give an inverse problem formulation to the regularized motion tracking problem (Section 2.4), and we describe our numerical strategy for the resolution (Section 2.5). We finish the Methods section with a description of the synthetic images that will be used for the validation of our method (Section 2.6), and then analyse tracking Results in the case of rigid motion (Section 3.1), non rigid but homogeneous motion (Section 3.2) and non homogeneous motion (Section 3.3), as well as in the case of refined meshes (Section 3.4).

## 2. Methods

### 2.1. The motion tracking problem

Let us start by precising the problem setting and notations. We consider $I_0$ and $I$, two images (i.e., image intensity fields) representing the same body $\mathcal{B}$ at two instants $t_0$ and $t$:

$$ I_0: \begin{cases} \square_0 \to \mathbb{R} \\ \underline{X} \mapsto I_0(\underline{X}), \end{cases} \quad I: \begin{cases} \square \to \mathbb{R} \\ \underline{x} \mapsto I(\underline{x}), \end{cases} \quad (1) $$

where $\square_0$ and $\square$ are the image domains at $t_0$ and $t$, which are usually identical. The domains occupied by the body $\mathcal{B}$ at $t_0$ and $t$ are denoted by $\Omega_0$ and $\omega$, respectively. The problem is to find the smooth mapping $\Phi$, or equivalently the smooth displacement field $\underline{U}$, between material points of the reference and deformed domains:

$$ \underline{\Phi}: \begin{cases} \Omega_0 \to \omega \\ \underline{X} \mapsto \underline{x} = \underline{\Phi}(\underline{X}), \end{cases} \quad \underline{U}: \begin{cases} \Omega_0 \to \mathbb{R}^3 \\ \underline{X} \mapsto \underline{U}(\underline{X}) := \underline{\Phi}(\underline{X}) - \underline{X}, \end{cases} \quad (2) $$

where $\underline{X}$ and $\underline{x}$ denote the position of a given material point in the reference and deformed configurations, respectively. Due to its intrinsic ill-posedness, the problem is formulated as a regularized minimization problem:

$$ \text{Find } \underline{U}^{\text{sol}} := \operatorname*{argmin}_{\underline{U}} \left\{ j(\underline{U}) := (1 - \beta) \frac{J^{\text{ima}}(\underline{U})}{J_0^{\text{ima}}} + \beta \frac{J^{\text{reg}}(\underline{U})}{J_0^{\text{reg}}} \right\}, \quad (3) $$

where $J^{\text{ima}}$ is the image similarity metric, or “correlation energy”, $J^{\text{reg}}$ is the regularization energy, $J_0^{\text{ima}}$ and $J_0^{\text{reg}}$ are normalization terms, and $\beta$ defines the regularization strength. The normalization terms $J_0^{\text{ima}}$ and $J_0^{\text{reg}}$ allow for the consistent addition of “energies” with very different physical units, and are typically taken equal to the value of $J^{\text{ima}}$ and $J^{\text{reg}}$ for a chosen displacement field [19], for instance a plane wave displacement with a period of 10 finite element characteristic lengths and a unit magnitude. The correlation energy is assumed to be convex, at least in the neighborhood of the solution, though it is in general not quadratic.

In image intensity-based global approaches, the following correlation energy is generally used:

$$ J^{\text{ima}}(\underline{U}) := \frac{1}{2} \int_{\Omega_0} (I(\underline{X} + \underline{U}(\underline{X})) - I_0(\underline{X}))^2 \, d\Omega_0. \quad (4) $$

Other metrics have been proposed; however, we retain this one notably because it can be differentiated straightforwardly.
option would be to use a separate discretization for the displacement and the tractions, which however would represent a significant technical difficulty. Instead, in this paper, we propose an additional regularization term based on the surface gradient of the normal and tangential components of the surface tractions, and show that it performs as expected.

The rest of the paper is organized as follows. We first formulate the general motion tracking problem (Section 2.1), then we provide a short literature review on mechanical regularization (Section 2.2), and we describe our proposed regularization term, including body (Section 2.3.1) and boundary (Section 2.3.2) terms. Then, we give an inverse problem formulation to the regularized motion tracking problem (Section 2.4), and we describe our numerical strategy for the resolution (Section 2.5). We finish the Methods section with a description of the synthetic images that will be used for the validation of our method (Section 2.6), and then analyse tracking Results in the case of rigid motion (Section 3.1), non rigid but homogeneous motion (Section 3.2) and non homogeneous motion (Section 3.3), as well as in the case of refined meshes (Section 3.4).

## 2. Methods

### 2.1. The motion tracking problem

Let us start by precising the problem setting and notations. We consider $I_0$ and $I$, two images (i.e., image intensity fields) representing the same body $\mathcal{B}$ at two instants $t_0$ and $t$:

$$ I_0: \begin{cases} \square_0 \to \mathbb{R} \\ \underline{X} \mapsto I_0(\underline{X}), \end{cases} \quad I: \begin{cases} \square \to \mathbb{R} \\ \underline{x} \mapsto I(\underline{x}), \end{cases} \quad (1) $$

where $\square_0$ and $\square$ are the image domains at $t_0$ and $t$, which are usually identical. The domains occupied by the body $\mathcal{B}$ at $t_0$ and $t$ are denoted by $\Omega_0$ and $\omega$, respectively. The problem is to find the smooth mapping $\Phi$, or equivalently the smooth displacement field $\underline{U}$, between material points of the reference and deformed domains:

$$ \underline{\Phi}: \begin{cases} \Omega_0 \to \omega \\ \underline{X} \mapsto \underline{x} = \underline{\Phi}(\underline{X}), \end{cases} \quad \underline{U}: \begin{cases} \Omega_0 \to \mathbb{R}^3 \\ \underline{X} \mapsto \underline{U}(\underline{X}) := \underline{\Phi}(\underline{X}) - \underline{X}, \end{cases} \quad (2) $$

where $\underline{X}$ and $\underline{x}$ denote the position of a given material point in the reference and deformed configurations, respectively. Due to its intrinsic ill-posedness, the problem is formulated as a regularized minimization problem:

$$ \text{Find } \underline{U}^{\text{sol}} := \operatorname*{argmin}_{\underline{U}} \left\{ j(\underline{U}) := (1 - \beta) \frac{J^{\text{ima}}(\underline{U})}{J_0^{\text{ima}}} + \beta \frac{J^{\text{reg}}(\underline{U})}{J_0^{\text{reg}}} \right\}, \quad (3) $$

where $J^{\text{ima}}$ is the image similarity metric, or “correlation energy”, $J^{\text{reg}}$ is the regularization energy, $J_0^{\text{ima}}$ and $J_0^{\text{reg}}$ are normalization terms, and $\beta$ defines the regularization strength. The normalization terms $J_0^{\text{ima}}$ and $J_0^{\text{reg}}$ allow for the consistent addition of “energies” with very different physical units, and are typically taken equal to the value of $J^{\text{ima}}$ and $J^{\text{reg}}$ for a chosen displacement field [19], for instance a plane wave displacement with a period of 10 finite element characteristic lengths and a unit magnitude. The correlation energy is assumed to be convex, at least in the neighborhood of the solution, though it is in general not quadratic.

In image intensity-based global approaches, the following correlation energy is generally used:

$$ J^{\text{ima}}(\underline{U}) := \frac{1}{2} \int_{\Omega_0} (I(\underline{X} + \underline{U}(\underline{X})) - I_0(\underline{X}))^2 \, d\Omega_0. \quad (4) $$

Other metrics have been proposed; however, we retain this one notably because it can be differentiated straightforwardly.
We will employ the finite element method to discretize this problem, such that the displacement field is approximated as $\underline{U}(\underline{X}) \approx {}^t\underline{\mathbb{N}}(\underline{X}) \cdot \underline{U}$ with $\underline{\mathbb{N}}$ the array of shape functions. Thus, the problem becomes a finite dimensional problem:

$$ \text{Find } \underline{U}^{\text{sol}} := \text{argmin}_{\underline{U}} \left\{ J(\underline{U}) := (1 - \beta) \frac{J^{\text{ima}}(\underline{U})}{J_0^{\text{ima}}} + \beta \frac{J^{\text{reg}}(\underline{U})}{J_0^{\text{reg}}} \right\}, \quad (5) $$

where

$$ J^{\text{ima}}(\underline{U}) := \frac{1}{2} \int_{\Omega_0} \left( I\left(\underline{X} + {}^t\underline{\mathbb{N}}(\underline{X}) \cdot \underline{U}\right) - I_0(\underline{X}) \right)^2 d\Omega_0. \quad (6) $$

Note that even without mechanical regularization, finite element discretization introduces some kind of kinematical regularization, as the richness of the approximation space is controlled by the mesh size and shape functions degree.

## 2.2. Short literature review on mechanical regularization

Many regularization terms have been proposed in the literature, see for instance [11]. Here we will briefly recall the major classes of proposals with mechanical content. All these approaches require to define a constitutive law, though it does not need to model the actual behavior of the tracked body: generic material laws can be used, and the material stiffness simply controls the strength of the regularization.

Elastic [25] and hyperelastic [16] regularizations have been proposed, which consist in penalizing the strain energy of the body:

$$ J^{\text{reg,el}}(\underline{U}) := \int_{\Omega_0} \rho_0 \Psi(\underline{U}) \, d\Omega_0, \quad (7) $$

where $\rho_0$ is the reference mass density and $\Psi$ the reference specific free energy—in practice, generic laws such as Hooke or neohookean are usually used. This was probably inspired by the fact that the elastostatic problem can often be formulated as a minimization problem with the system potential energy that is the sum of the elastic and loading potential energies. However, by only considering the elastic energy beside the image similarity metric, and no specific loading energy, it is implicitly assumed that the only load applied to the body is a body force associated to image dissimilarity, and no boundary tractions. Moreover, this regularization penalizes strain itself, as only rigid body motions have zero elastic energy.

Another, arguably more mechanically consistent, regularization approach is based on the equilibrium gap principle [18]. It was formulated in the linear setting, and directly at the discrete level, in [19, 26]. We first define $\underline{K} := \int_{\Omega_0} \underline{\mathbb{B}} : \underline{K} : {}^t\underline{\mathbb{B}} \, d\Omega_0$ the system stiffness matrix, with $\underline{K}$ the material stiffness tensor (like for elastic regularization, a generic isotropic Hooke law is usually used, such that $\underline{K} = \lambda\underline{1} \otimes \underline{1} + 2\mu\underline{1}$ with $\lambda$ and $\mu$ the Lamé constants) and $\underline{\mathbb{B}}$ the array of shape functions symmetric gradients (i.e., such that $\underline{\varepsilon}(\underline{X}) := (\underline{\text{Grad}}(\underline{U})(\underline{X}))_{\text{sym}} = (1/2)({}^t\underline{\text{Grad}}(\underline{U})(\underline{X}) + \underline{\text{Grad}}(\underline{U})(\underline{X})) \approx {}^t\underline{\mathbb{B}}(\underline{X}) \cdot \underline{U}$). Then the regularization is expressed as

$$ J^{\text{reg,eq,lin}}(\underline{U}) := \frac{1}{2} {}^t\underline{U} \cdot {}^t\underline{K}^* \cdot \underline{K}^* \cdot \underline{U}, \quad (8) $$

where $\underline{K}^*$ is a modified system stiffness matrix in which all lines associated to boundary degrees of freedom at which the external force in unknown have been set to 0. The implicit assumption here is that the body is in equilibrium with some arbitrary/unknown boundary tractions (though the smoothness of the surface displacement is usually controlled by an additional surface Laplacian term [19]) and no body force (though known body forces such as gravity
We will employ the finite element method to discretize this problem, such that the displacement field is approximated as $\underline{U}(\underline{X}) \approx {}^t\underline{\mathbb{N}}(\underline{X}) \cdot \underline{U}$ with $\underline{\mathbb{N}}$ the array of shape functions. Thus, the problem becomes a finite dimensional problem:

$$ \text{Find } \underline{U}^{\text{sol}} := \text{argmin}_{\underline{U}} \left\{ J(\underline{U}) := (1 - \beta) \frac{J^{\text{ima}}(\underline{U})}{J_0^{\text{ima}}} + \beta \frac{J^{\text{reg}}(\underline{U})}{J_0^{\text{reg}}} \right\}, \quad (5) $$

where

$$ J^{\text{ima}}(\underline{U}) := \frac{1}{2} \int_{\Omega_0} \left( I\left(\underline{X} + {}^t\underline{\mathbb{N}}(\underline{X}) \cdot \underline{U}\right) - I_0(\underline{X}) \right)^2 d\Omega_0. \quad (6) $$

Note that even without mechanical regularization, finite element discretization introduces some kind of kinematical regularization, as the richness of the approximation space is controlled by the mesh size and shape functions degree.

## 2.2. Short literature review on mechanical regularization

Many regularization terms have been proposed in the literature, see for instance [11]. Here we will briefly recall the major classes of proposals with mechanical content. All these approaches require to define a constitutive law, though it does not need to model the actual behavior of the tracked body: generic material laws can be used, and the material stiffness simply controls the strength of the regularization.

Elastic [25] and hyperelastic [16] regularizations have been proposed, which consist in penalizing the strain energy of the body:

$$ J^{\text{reg,el}}(\underline{U}) := \int_{\Omega_0} \rho_0 \Psi(\underline{U}) d\Omega_0, \quad (7) $$

where $\rho_0$ is the reference mass density and $\Psi$ the reference specific free energy—in practice, generic laws such as Hooke or neohookean are usually used. This was probably inspired by the fact that the elastostatic problem can often be formulated as a minimization problem with the system potential energy that is the sum of the elastic and loading potential energies. However, by only considering the elastic energy beside the image similarity metric, and no specific loading energy, it is implicitly assumed that the only load applied to the body is a body force associated to image dissimilarity, and no boundary tractions. Moreover, this regularization penalizes strain itself, as only rigid body motions have zero elastic energy.

Another, arguably more mechanically consistent, regularization approach is based on the equilibrium gap principle [18]. It was formulated in the linear setting, and directly at the discrete level, in [19, 26]. We first define $\underline{K} := \int_{\Omega_0} \underline{\mathbb{B}} : \underline{K} : {}^t\underline{\mathbb{B}} d\Omega_0$ the system stiffness matrix, with $\underline{K}$ the material stiffness tensor (like for elastic regularization, a generic isotropic Hooke law is usually used, such that $\underline{K} = \lambda\underline{1} \otimes \underline{1} + 2\mu\underline{1}$ with $\lambda$ and $\mu$ the Lamé constants) and $\underline{\mathbb{B}}$ the array of shape functions symmetric gradients (i.e., such that $\underline{\varepsilon}(\underline{X}) := (\underline{\text{Grad}}(\underline{U})(\underline{X}))_{\text{sym}} = (1/2)({}^t\underline{\text{Grad}}(\underline{U})(\underline{X}) + \underline{\text{Grad}}(\underline{U})(\underline{X})) \approx {}^t\underline{\mathbb{B}}(\underline{X}) \cdot \underline{U}$). Then the regularization is expressed as

$$ J^{\text{reg,eq,lin}}(\underline{U}) := \frac{1}{2} {}^t\underline{U} \cdot {}^t\underline{K}^* \cdot \underline{K}^* \cdot \underline{U}, \quad (8) $$

where $\underline{K}^*$ is a modified system stiffness matrix in which all lines associated to boundary degrees of freedom at which the external force in unknown have been set to 0. The implicit assumption here is that the body is in equilibrium with some arbitrary/unknown boundary tractions (though the smoothness of the surface displacement is usually controlled by an additional surface Laplacian term [19]) and no body force (though known body forces such as gravity
could easily be taken into account). Thus, it does not penalize strain, only deviation from equilibrium, as any equilibrium solution, even with very large strain, cancels the equilibrium gap "energy".

A first attempt toward extending this principle to the large deformation nonlinear setting has been made in [17]. The main idea is to directly penalize the non verification of the internal linear momentum equilibrium equation ($\underline{\underline{\text{Div}(P)}}$ = 0, where $P$ is the first Piola-Kirchhoff stress tensor associated to the displacement field $\underline{U}$ through the chosen constitutive law) through the following equivalent norm:

$$J^{\text{reg,eq,cont}}(\underline{U}) := \frac{1}{2} \sum_{E} \int_{E} \|\underline{\underline{\text{Div}(P)}}\|^2 dE + \frac{1}{2h} \sum_{F} \int_{F} \|[\underline{P} \cdot \underline{\underline{N}}]\|^2 dF \quad (9)$$

where $E$ is the set of finite elements, $F$ the set of interior faces, $\underline{N}$ the faces normal, and $h$ a characteristic length of the finite element discretization. Note that the internal angular momentum equilibrium ($\underline{tP} \cdot \underline{tF}^{-1} = \underline{F}^{-1} \cdot \underline{P}$, where $\underline{F}$ is the deformation gradient associated to the displacement field $\underline{U}$) is usually exactly verified through the constitutive relation, hence it does not need to be penalized. Thanks to this equivalent norm (inspired by discontinuous Galerkin approaches; it is necessary because under standard finite element discretization the finite element stress field is discontinuous across finite element faces, so that its divergence is not defined there, and thus cannot be directly integrated over the mesh), all terms of the regularization (9) are proper element or face integrals and can be computed using standard finite element integration techniques (for instance, in practice, if one uses first order polynomial shape functions, the stress field is constant by element, so that the first term vanishes, and one simply has to compute the normal stress jump across all internal faces of the mesh). Thus, the underlying assumption is the same as with regularization (8), i.e., that the body is in equilibrium with some arbitrary boundary tractions and no body force, although the geometrically nonlinear formulation allows to correctly handle large deformations, including large rotations. Nevertheless, one limitation of this approach is that the penalization (9) encompasses not just the equilibrium gap arising from the potential non-equilibrium nature of the considered displacement field (i.e., as it might not be an equilibrium solution of a mechanical problem), but also accounts for the equilibrium gap resulting from the finite element discretization itself (i.e., even an actual finite element solution of a mechanical problem leads to a nonzero regularization energy (9), as a finite element solution verifies the equilibrium equation only weakly, not strongly). This calls for an improvement of the approach, which is described in the following section.

## 2.3. The nonlinear discrete equilibrium gap regularization

### 2.3.1. Body stresses regularization

In order to avoid the limitation of the regularization term (9), we propose here a novel formulation of the equilibrium gap principle, still in the general context of nonlinear mechanics, but which allows to completely exclude the discretization-induced equilibrium gap from the regularization term. The general idea is still to penalize the non-verification of the internal linear momentum balance ($\underline{\text{Div}(P)}$ = 0 in $\Omega_0$, since the body force is neglected for the sake of simplicity, though it could be considered as well). However, after standard finite element discretization, the divergence of the stress tensor is not defined on the element faces (in 3D) or edges (in 2D) thus it is not square integrable, so we define its projection $\underline{\underline{\Pi}}_b$ onto a space of square integrable functions, typically a standard continuous finite element space denoted $V^h$:

$$\underline{\underline{\Pi}}_b \in V^h \mid \int_{\Omega_0} \underline{\underline{\Pi}}_b \cdot \underline{\underline{\Pi}}_b^* = - \int_{\Omega_0} \underline{\underline{\text{Div}(P)}} \cdot \underline{\underline{\Pi}}_b^* \quad \forall \underline{\underline{\Pi}}_b^* \in V_0^h, \quad (10)$$
could easily be taken into account). Thus, it does not penalize strain, only deviation from equilibrium, as any equilibrium solution, even with very large strain, cancels the equilibrium gap “energy”.

A first attempt toward extending this principle to the large deformation nonlinear setting has been made in [17]. The main idea is to directly penalize the non verification of the internal linear momentum equilibrium equation ($\underline{\underline{\text{Div}(P)}}$ = 0, where $P$ is the first Piola-Kirchhoff stress tensor associated to the displacement field $\underline{U}$ through the chosen constitutive law) through the following equivalent norm:

$$J^{\text{reg,eq,cont}}(\underline{U}) := \frac{1}{2} \sum_{E} \int_{E} \|\underline{\underline{\text{Div}(P)}}\|^2 dE + \frac{1}{2h} \sum_{F} \int_{F} \|[\underline{P} \cdot \underline{\underline{N}}]\|^2 dF \quad (9)$$

where $E$ is the set of finite elements, $F$ the set of interior faces, $\underline{N}$ the faces normal, and $h$ a characteristic length of the finite element discretization. Note that the internal angular momentum equilibrium ($\underline{\underline{tP}} \cdot \underline{\underline{tF}}^{-1} = \underline{\underline{F}}^{-1} \cdot \underline{P}$, where $\underline{F}$ is the deformation gradient associated to the displacement field $\underline{U}$) is usually exactly verified through the constitutive relation, hence it does not need to be penalized. Thanks to this equivalent norm (inspired by discontinuous Galerkin approaches; it is necessary because under standard finite element discretization the finite element stress field is discontinuous across finite element faces, so that its divergence is not defined there, and thus cannot be directly integrated over the mesh), all terms of the regularization (9) are proper element or face integrals and can be computed using standard finite element integration techniques (for instance, in practice, if one uses first order polynomial shape functions, the stress field is constant by element, so that the first term vanishes, and one simply has to compute the normal stress jump across all internal faces of the mesh). Thus, the underlying assumption is the same as with regularization (8), i.e., that the body is in equilibrium with some arbitrary boundary tractions and no body force, although the geometrically nonlinear formulation allows to correctly handle large deformations, including large rotations. Nevertheless, one limitation of this approach is that the penalization (9) encompasses not just the equilibrium gap arising from the potential non-equilibrium nature of the considered displacement field (i.e., as it might not be an equilibrium solution of a mechanical problem), but also accounts for the equilibrium gap resulting from the finite element discretization itself (i.e., even an actual finite element solution of a mechanical problem leads to a nonzero regularization energy (9), as a finite element solution verifies the equilibrium equation only weakly, not strongly). This calls for an improvement of the approach, which is described in the following section.

## 2.3. The nonlinear discrete equilibrium gap regularization

### 2.3.1. Body stresses regularization

In order to avoid the limitation of the regularization term (9), we propose here a novel formulation of the equilibrium gap principle, still in the general context of nonlinear mechanics, but which allows to completely exclude the discretization-induced equilibrium gap from the regularization term. The general idea is still to penalize the non-verification of the internal linear momentum balance ($\underline{\underline{\text{Div}(P)}}$ = 0 in $\Omega_0$, since the body force is neglected for the sake of simplicity, though it could be considered as well). However, after standard finite element discretization, the divergence of the stress tensor is not defined on the element faces (in 3D) or edges (in 2D) thus it is not square integrable, so we define its projection $\underline{\underline{\Pi}}_b$ onto a space of square integrable functions, typically a standard continuous finite element space denoted $V^h$:

$$\underline{\underline{\Pi}}_b \in V^h \mid \int_{\Omega_0} \underline{\underline{\Pi}}_b \cdot \underline{\underline{\Pi}}_b^* = - \int_{\Omega_0} \underline{\underline{\text{Div}(P)}} \cdot \underline{\underline{\Pi}}_b^* \quad \forall \underline{\underline{\Pi}}_b^* \in V_0^h, \quad (10)$$
where $V_0^h$ denotes the space of functions of $V^h$ that vanish at the boundary, i.e., $\Pi_b^*(\partial\Omega_0) = 0$. Indeed, for now we are only interested in body equilibrium gap. Thus, after integration by parts, we obtain:

$$ \underline{\Pi}_b \in V^h \mid \int_{\Omega_0} \underline{\Pi}_b \cdot \underline{\Pi}_b^* = \int_{\Omega_0} \underline{P} : \underline{\text{Grad}}(\underline{\Pi}_b^*) \quad \forall \underline{\Pi}_b^* \in V_0^h. \quad (11) $$

Similar projections, though often performed element-by-element, are used in a *posteriori* error estimation methods, to quantify the distance between the finite element solution and an actual equilibrium solution, i.e., the discretization error [27,28]. However, the objective here is opposite, i.e., we want to discard the discretization error and quantify the intrinsic equilibrium gap of the considered finite element field.

Nevertheless, thanks to the proposed projection, we can actually define the regularization term in a consistent manner:

$$ J_b^{\text{reg,eq}}(\underline{U}) := \frac{1}{2} \int_{\Omega_0} \underline{\Pi}_b \cdot \underline{\Pi}_b \, d\Omega_0. \quad (12) $$

A key question is the calculation of this term in the finite element context. Projection (11) simply leads to the following linear system:

$$ \underline{M} \cdot \underline{\Pi}_b = \underline{R}_b \quad (13) $$

with $\underline{M} := \int_{\Omega_0} \underline{N} \cdot \underline{t} \underline{N} d\Omega_0$ the mass matrix, $\underline{\Pi}_b$ such that $\underline{\Pi}_b = \underline{t} \underline{N} \cdot \underline{\Pi}_b$, and

$$ (\underline{R}_b)_i := \begin{cases} \int_{\Omega_0} \underline{P} : \underline{\text{Grad}}(\underline{N}_i) \, d\Omega_0 & \text{if } i \text{ body d.o.f.} \\ 0 & \text{if } i \text{ boundary d.o.f.,} \end{cases} \quad (14) $$

where $\underline{N}_i$ is the (vector) shape function associated to the degree of freedom (d.o.f.) $i$, i.e., the $i$th line of the $\underline{N}$ array. Thus, the norm (12) can be expressed as

$$ J_b^{\text{reg,eq}}(\underline{U}) = \frac{1}{2} \underline{t} \underline{\Pi}_b \cdot \underline{M} \cdot \underline{\Pi}_b = \frac{1}{2} \underline{t} \underline{R}_b \cdot \underline{M}^{-1} \cdot \underline{R}_b. \quad (15) $$

Note that if we linearize this expression, we obtain $\underline{R}_b \approx \underline{K}^* \cdot \underline{U}$ and thus $J_b^{\text{reg,eq}} \approx (1/2)^t \underline{U} \cdot \underline{t} \underline{K}^* \cdot \underline{M}^{-1} \cdot \underline{K}^* \cdot \underline{U}$, i.e., an expression similar to (8), the mass matrix allowing to make the term consistent when refining the mesh.

### 2.3.2. Boundary tractions regularization

The regularization term proposed in Section 2.3.1 basically enforces that the motion solution corresponds to the motion of a body in equilibrium with some arbitrary tractions applied on its boundary. The arbitrary nature of these tractions can be problematic, especially as the computational mesh is refined and the variations of these tractions is not controlled. In [19], it was proposed to add a penalization term corresponding to the Laplacian of the boundary displacement, which however does not generalize to the large motion context. Thus, here we propose a new term, consistent with the body term that involves internal stresses directly, which consists in penalizing the surface tractions gradients, an idea already considered in [29].

It is important, however, to note that surface tractions can vary intrinsically but also because of the surface curvature. For instance, a simple homogeneous pressure applied onto a curved surface corresponds to a vector that varies in space. To avoid penalizing the surface curvature itself, we propose to penalize the surface gradient of the normal and tangential components of the boundary tractions separately. In 2D, the tangential part is scalar; in 3D we propose to take the norm of the tangential force vector.

Let us start with the regularization term associated to the normal traction, i.e., $F_n(\underline{U}) := \underline{t} \underline{N} \cdot \underline{P}(\underline{U}) \cdot \underline{N}$ where $\underline{N}$ is the body outward normal. The problem is the same as for the equilibrium equation, i.e., the surface gradient of the boundary tractions associated to a standard finite element displacement field is not defined on the boundary elements edges (in 3D) or points
where $V_0^h$ denotes the space of functions of $V^h$ that vanish at the boundary, i.e., $\Pi_b^*(\partial\Omega_0) = 0$. Indeed, for now we are only interested in body equilibrium gap. Thus, after integration by parts, we obtain:

$$ \underline{\Pi}_b \in V^h \mid \int_{\Omega_0} \underline{\Pi}_b \cdot \underline{\Pi}_b^* = \int_{\Omega_0} \underline{P} : \underline{\text{Grad}}(\underline{\Pi}_b^*) \quad \forall \underline{\Pi}_b^* \in V_0^h. \quad (11) $$

Similar projections, though often performed element-by-element, are used in a *posteriori* error estimation methods, to quantify the distance between the finite element solution and an actual equilibrium solution, i.e., the discretization error [27,28]. However, the objective here is opposite, i.e., we want to discard the discretization error and quantify the intrinsic equilibrium gap of the considered finite element field.

Nevertheless, thanks to the proposed projection, we can actually define the regularization term in a consistent manner:

$$ J_b^{\text{reg,eq}}(\underline{U}) := \frac{1}{2} \int_{\Omega_0} \underline{\Pi}_b \cdot \underline{\Pi}_b \, d\Omega_0. \quad (12) $$

A key question is the calculation of this term in the finite element context. Projection (11) simply leads to the following linear system:

$$ \underline{M} \cdot \underline{\Pi}_b = \underline{R}_b \quad (13) $$

with $\underline{M} := \int_{\Omega_0} \underline{N} \cdot \underline{t} \underline{N} d\Omega_0$ the mass matrix, $\underline{\Pi}_b$ such that $\underline{\Pi}_b = \underline{t} \underline{N} \cdot \underline{\Pi}_b$, and

$$ (\underline{R}_b)_i := \begin{cases} \int_{\Omega_0} \underline{P} : \underline{\text{Grad}}(\underline{N}_i) \, d\Omega_0 & \text{if } i \text{ body d.o.f.} \\ 0 & \text{if } i \text{ boundary d.o.f.,} \end{cases} \quad (14) $$

where $\underline{N}_i$ is the (vector) shape function associated to the degree of freedom (d.o.f.) $i$, i.e., the $i$th line of the $\underline{N}$ array. Thus, the norm (12) can be expressed as

$$ J_b^{\text{reg,eq}}(\underline{U}) = \frac{1}{2} \underline{t} \underline{\Pi}_b \cdot \underline{M} \cdot \underline{\Pi}_b = \frac{1}{2} \underline{t} \underline{R}_b \cdot \underline{M}^{-1} \cdot \underline{R}_b. \quad (15) $$

Note that if we linearize this expression, we obtain $\underline{R}_b \approx \underline{K}^* \cdot \underline{U}$ and thus $J_b^{\text{reg,eq}} \approx (1/2)^t \underline{U} \cdot \underline{t} \underline{K}^* \cdot \underline{M}^{-1} \cdot \underline{K}^* \cdot \underline{U}$, i.e., an expression similar to (8), the mass matrix allowing to make the term consistent when refining the mesh.

### 2.3.2. Boundary tractions regularization

The regularization term proposed in Section 2.3.1 basically enforces that the motion solution corresponds to the motion of a body in equilibrium with some arbitrary tractions applied on its boundary. The arbitrary nature of these tractions can be problematic, especially as the computational mesh is refined and the variations of these tractions is not controlled. In [19], it was proposed to add a penalization term corresponding to the Laplacian of the boundary displacement, which however does not generalize to the large motion context. Thus, here we propose a new term, consistent with the body term that involves internal stresses directly, which consists in penalizing the surface tractions gradients, an idea already considered in [29].

It is important, however, to note that surface tractions can vary intrinsically but also because of the surface curvature. For instance, a simple homogeneous pressure applied onto a curved surface corresponds to a vector that varies in space. To avoid penalizing the surface curvature itself, we propose to penalize the surface gradient of the normal and tangential components of the boundary tractions separately. In 2D, the tangential part is scalar; in 3D we propose to take the norm of the tangential force vector.

Let us start with the regularization term associated to the normal traction, i.e., $F_n(\underline{U}) := \underline{t} \underline{N} \cdot \underline{P}(\underline{U}) \cdot \underline{N}$ where $\underline{N}$ is the body outward normal. The problem is the same as for the equilibrium equation, i.e., the surface gradient of the boundary tractions associated to a standard finite element displacement field is not defined on the boundary elements edges (in 3D) or points
(in 2D) thus it is not integrable. Hence we propose to use the same technique, i.e., we first define its projection onto a space of square integrable functions, again a typical finite element space, continuous on the domain boundary and denoted $\partial V^h$:

$$ \underline{\Pi}_n \in \partial V^h \left| \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{\Pi}_n^* = - \int_{\partial \Omega_0} \underline{\mathrm{Grad}}_s(F_n) \cdot \underline{\Pi}_n^* \quad \forall \underline{\Pi}_n^* \in \partial V^h, \right. \quad (16) $$

where $\underline{\mathrm{Grad}}_s(F_n) = \underline{\Pi} \cdot \underline{\mathrm{Grad}}(F_n)$ with $\underline{\Pi} := \underline{1} - \underline{N} \otimes \underline{N}$ the projection operator onto the domain boundary [30]. After integration by parts of the right hand side we obtain:

$$ \underline{\Pi}_n \in \partial V^h \left| \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{\Pi}_n^* = \int_{\partial \Omega_0} F_n \cdot \mathrm{Div}_s(\underline{\Pi}_n^*) \quad \forall \underline{\Pi}_n^* \in \partial V^h, \right. \quad (17) $$

where $\mathrm{Div}_s(\underline{\Pi}_n^*) = \mathrm{tr}(\underline{\Pi} \cdot \underline{\mathrm{Grad}}(\underline{\Pi}_n^*) \underline{\Pi})$ [30]. Then, the regularization term is actually defined as:

$$ J_n^{\mathrm{reg,eq}}(\underline{U}) := \frac{1}{2} \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{\Pi}_n d\partial \Omega_0. \quad (18) $$

The discretization procedure is similar to the bulk term (15):

$$ J_n^{\mathrm{reg,eq}}(\underline{U}) = \frac{1}{2} t \underline{\mathbb{R}}_n \cdot \underline{\mathbb{M}}_{\partial \Omega_0}^{-1} \cdot \underline{\mathbb{R}}_n \quad (19) $$

where $\underline{\mathbb{M}}_{\partial \Omega_0} := \int_{\partial \Omega_0} \underline{\mathbb{N}} \cdot t \underline{\mathbb{N}} d\partial \Omega_0$ is the “mass” matrix of the domain boundary, and

$$ (\underline{\mathbb{R}}_n)_i := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \int_{\partial \Omega_0} F_n \cdot \mathrm{Div}_s(\underline{\mathbb{N}}_i) d\partial \Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (20) $$

Regarding the regularization term associated to tangential tractions, in 2D the tangential force is a scalar, defined as $F_{t,2D}(\underline{U}) := t \underline{\underline{T}} \cdot \underline{\underline{P}}(\underline{U}) \cdot \underline{\underline{N}}$ where $\underline{\underline{T}}$ is the body tangential vector, so it is the same formulation as the normal traction (19):

$$ J_{t,2D}^{\mathrm{reg,eq}}(\underline{U}) = \frac{1}{2} t \underline{\underline{\mathbb{R}}}_{t,2D} \cdot \underline{\underline{\mathbb{M}}}_{\partial \Omega_0}^{-1} \cdot \underline{\underline{\mathbb{R}}}_{t,2D} \quad (21) $$

where

$$ (\underline{\mathbb{R}}_{t,2D})_i := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \int_{\partial \Omega_0} F_{t,2D} \cdot \mathrm{Div}_s(\underline{\mathbb{N}}_i) d\partial \Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (22) $$

In 3D the tangential force is a vector, defined as $F_{t,3D}(\underline{U}) := \underline{\underline{\Pi}} \cdot (\underline{\underline{P}}(\underline{U}) \cdot \underline{\underline{N}})$, and we propose to simply penalize the gradient of its norm, which leads to:

$$ J_{t,3D}^{\mathrm{reg,eq}}(\underline{U}) = \frac{1}{2} t \underline{\underline{\mathbb{R}}}_{t,3D} \cdot \underline{\underline{\mathbb{M}}}_{\partial \Omega_0}^{-1} \cdot \underline{\underline{\mathbb{R}}}_{t,3D} \quad (23) $$

where

$$ (\underline{\mathbb{R}}_{t,3D})_i := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \int_{\partial \Omega_0} \|F_{t,3D}\| \cdot \mathrm{Div}_s(\underline{\mathbb{N}}_i) d\partial \Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (24) $$

## 2.4. Inverse problem reformulation

Because the proposed regularization has a strong mechanical sense, the regularized tracking problem can be reformulated as an inverse problem, where the unknown is the traction field applied on the domain boundary, and which can be formulated as a constrained optimization problem:

$$ \text{Find } \underline{T}^{\mathrm{sol}} := \underset{(\underline{T})}{\operatorname{argmin}} \{ J(\underline{T}) := (1-\beta)\|\underline{I} \circ \underline{\Phi}(\underline{T}) - \underline{I}_0\|_{\Omega_0}^2 + \beta\|\underline{T}\|_{\partial \Omega_0}^2 \}, \quad (25) $$

where $\underline{\Phi}(\underline{T})$ is the mapping solution associated to the traction field $\underline{T}$, i.e., which verifies constitutive and equilibrium equations, and $\beta$ is the regularization strength. Norms will be specified later on. The meaning of this formulation is that we search for a traction field which generates
(in 2D) thus it is not integrable. Hence we propose to use the same technique, i.e., we first define its projection onto a space of square integrable functions, again a typical finite element space, continuous on the domain boundary and denoted $\partial V^h$:

$$ \underline{\Pi}_n \in \partial V^h \left| \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{\Pi}_n^* = - \int_{\partial \Omega_0} \underline{\mathrm{Grad}}_s(F_n) \cdot \underline{\Pi}_n^* \quad \forall \underline{\Pi}_n^* \in \partial V^h, \quad (16) $$

where $\underline{\mathrm{Grad}}_s(F_n) = \underline{\Pi} \cdot \underline{\mathrm{Grad}}(F_n)$ with $\underline{\Pi} := \underline{1} - \underline{N} \otimes \underline{N}$ the projection operator onto the domain boundary [30]. After integration by parts of the right hand side we obtain:

$$ \underline{\Pi}_n \in \partial V^h \left| \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{\Pi}_n^* = \int_{\partial \Omega_0} F_n \cdot \mathrm{Div}_s(\underline{\Pi}_n^*) \quad \forall \underline{\Pi}_n^* \in \partial V^h, \quad (17) $$

where $\mathrm{Div}_s(\underline{\Pi}_n^*) = \mathrm{tr}(\underline{\Pi} \cdot \underline{\mathrm{Grad}}(\underline{\Pi}_n^*) \cdot \underline{\Pi})$ [30]. Then, the regularization term is actually defined as:

$$ J_n^{\mathrm{reg,eq}}(\underline{U}) := \frac{1}{2} \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{\Pi}_n d\partial \Omega_0. \quad (18) $$

The discretization procedure is similar to the bulk term (15):

$$ J_n^{\mathrm{reg,eq}}(\underline{U}) = \frac{1}{2} \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{M}^{-1}_{\partial \Omega_0} \cdot \underline{R}_n \quad (19) $$

where $\underline{M}_{\partial \Omega_0} := \int_{\partial \Omega_0} \underline{N} \cdot \underline{t} \underline{N} d\partial \Omega_0$ is the “mass” matrix of the domain boundary, and

$$ (\underline{\mathbb{R}}_n)_i := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \int_{\partial \Omega_0} F_n \cdot \mathrm{Div}_s(\underline{\mathbb{N}}_i) d\partial \Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (20) $$

Regarding the regularization term associated to tangential tractions, in 2D the tangential force is a scalar, defined as $F_{t,2D}(\underline{U}) := \underline{t} \underline{T} \cdot \underline{P}(\underline{U}) \cdot \underline{N}$ where $\underline{T}$ is the body tangential vector, so it is the same formulation as the normal traction (19):

$$ J_{t,2D}^{\mathrm{reg,eq}}(\underline{U}) = \frac{1}{2} \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{M}^{-1}_{\partial \Omega_0} \cdot \underline{R}_{t,2D} \quad (21) $$

where

$$ (\underline{\mathbb{R}}_{t,2D})_i := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \int_{\partial \Omega_0} F_{t,2D} \cdot \mathrm{Div}_s(\underline{\mathbb{N}}_i) d\partial \Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (22) $$

In 3D the tangential force is a vector, defined as $F_{t,3D}(\underline{U}) := \underline{\Pi} \cdot (\underline{P}(\underline{U}) \cdot \underline{N})$, and we propose to simply penalize the gradient of its norm, which leads to:

$$ J_{t,3D}^{\mathrm{reg,eq}}(\underline{U}) = \frac{1}{2} \int_{\partial \Omega_0} \underline{\Pi}_n \cdot \underline{M}^{-1}_{\partial \Omega_0} \cdot \underline{R}_{t,3D} \quad (23) $$

where

$$ (\underline{\mathbb{R}}_{t,3D})_i := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \int_{\partial \Omega_0} \|F_{t,3D}\| \cdot \mathrm{Div}_s(\underline{\mathbb{N}}_i) d\partial \Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (24) $$

## 2.4. Inverse problem reformulation

Because the proposed regularization has a strong mechanical sense, the regularized tracking problem can be reformulated as an inverse problem, where the unknown is the traction field applied on the domain boundary, and which can be formulated as a constrained optimization problem:

$$ \text{Find } \underline{T}^{\mathrm{sol}} := \underset{\underline{T}}{\mathrm{argmin}} \{ J(\underline{T}) := (1-\beta)\|\underline{I} \circ \underline{\Phi}(\underline{T}) - \underline{I}_0\|_{\Omega_0}^2 + \beta\|\underline{T}\|_{\partial \Omega_0}^2 \}, \quad (25) $$

where $\underline{\Phi}(\underline{T})$ is the mapping solution associated to the traction field $\underline{T}$, i.e., which verifies constitutive and equilibrium equations, and $\beta$ is the regularization strength. Norms will be specified later on. The meaning of this formulation is that we search for a traction field which generates
a displacement field (through a chosen generic mechanical behavior and standard mechanical equilibrium) that allows to match the two image intensities, the regularization term allowing to control the smoothness of the traction field. In optimization, such constraints can be enforced strongly through Lagrange multipliers, or approximately through penalization [31]. Here we propose to simply use penalization, such that the problem becomes:

$$ \text{Find } \underline{T}^{\text{sol}}, \underline{U}^{\text{sol}} := \operatorname*{argmin}_{\underline{T}, \underline{U}} \left\{ J(\underline{T}, \underline{U}) := (1-\beta)\|\underline{I} \circ \underline{\Phi} - \underline{I}_0\|_{\underline{\Omega}_0}^2 \right. \\ \left. + \beta \|\underline{T}\|_{\underline{\partial}\underline{\Omega}_0}^2 + \gamma \left( \|\underline{\mathrm{Div}}(\underline{P})\|_{\underline{\Omega}_0}^2 + \|{}^t\underline{P} \cdot {}^t\underline{F}^{-1} - \underline{F}^{-1} \cdot \underline{P}\|_{\underline{\Omega}_0}^2 + \|\underline{P} \cdot \underline{N} - \underline{T}\|_{\underline{\partial}\underline{\Omega}_0}^2 \right) \right\}, \quad (26) $$

where $\gamma$ is the penalization coefficient, while the three additional terms represent the balance of linear momentum (without imposed force, though it could be introduced straightforwardly), the balance of angular momentum (which is usually verified exactly thanks to the constitutive framework), and the balance with applied tractions, respectively. In this formulation, $\underline{\Phi}$, $\underline{F}$ and $\underline{P}$ are the mapping, deformation gradient and first Piola-Kirchhoff stress tensor associated to the displacement field $\underline{U}$. Again, the norms will be specified later on. To simplify the formulation, we can remove the traction field variable by enforcing strongly the last balance term, i.e., $\underline{T} = \underline{P}(\underline{U}) \cdot \underline{N}$, leading to the following formulation:

$$ \text{Find } \underline{U}^{\text{sol}} := \operatorname*{argmin}_{\underline{U}} \left\{ J(\underline{U}) := (1-\beta)\|\underline{I} \circ \underline{\Phi} - \underline{I}_0\|_{\underline{\Omega}_0}^2 + \beta\|\underline{P} \cdot \underline{N}\|_{\underline{\partial}\underline{\Omega}_0}^2 + \gamma\|\underline{\mathrm{Div}}(\underline{P})\|_{\underline{\Omega}_0}^2 \right\}. \quad (27) $$

This formulation corresponds formally (i.e., using the right norms) to the original formulation (3), when using $J_b^{\text{reg,eq}}$ for the bulk regularization and $J_{n/t}^{\text{reg,eq}}$ for the surface terms. This gives another point of view on each term, i.e., that the bulk term can be seen as a penalization term, while the actual regularization comes from the surface terms. Also, it is important to notice that this reformulation only makes sense for “equilibrium gap” regularization terms: if one would consider the motion tracking problem (3) with elastic regularization term (7) as an inverse problem, the associated direct problem would not have a proper physical meaning, as it would correspond to a body deforming under a body force induced by the image mismatch and no boundary force.

## 2.5. Numerical resolution

After describing existing and new regularization terms, we recall that the general motion tracking problem formulation was given by Equation (5), with the image correlation term $J^{\text{ima}}$ given by Equation (6), and where the regularization term $J^{\text{reg}}$ may take the form $J^{\text{reg,el}}$ (Equation (7)), $J^{\text{reg,eq,cont}}$ (Equation (9)), $J_b^{\text{reg,eq}}$ (Equation (12)), $J_n^{\text{reg,eq}}$ (Equation (18)), $J_t^{\text{reg,eq}}$ (Equation (22) in 2D, Equation (24) in 3D), or combinations of such terms. Indeed, as already discussed, it is often necessary to combine bulk and boundary terms. In such cases one might want to use different weights for each regularization term [19]; however, to simplify the analysis of the performance of the various terms considered here, we propose to simply sum them, such that they have the same weight, and there is only one parameter controlling the regularization strength, namely $\beta$.

Minimization problem (5) is actually formulated as a root finding problem:

$$ \underline{U} | \underline{\nabla} \underline{J}(\underline{U}) = (1 - \beta) \underline{\nabla} \underline{J}^{\text{ima}}(\underline{U}) + \beta \underline{\nabla} \underline{J}^{\text{reg}}(\underline{U}) = 0 \quad (28) $$

where the gradient of the image term is simply:

$$ \underline{\nabla} \underline{J}^{\text{ima}}(\underline{U}) := \int_{\underline{\Omega}_0} \left( I \left( \underline{X} + {}^t\underline{N} \cdot \underline{U} \right) - \underline{I}_0 \right) \frac{\partial I}{\partial \underline{x}} \left( \underline{X} + {}^t\underline{N} \cdot \underline{U} \right) \cdot \underline{N} \, d\underline{\Omega}_0, \quad (29) $$

and the gradients of the proposed regularization terms are:

$$ \underline{\nabla}_{\underline{b}}^{\text{reg,eq}}(\underline{U}) := {}^t\underline{d}\underline{R}_{\underline{b}}(\underline{U}) \cdot \underline{M}^{-1} \cdot \underline{R}_{\underline{b}}(\underline{U}) \quad (30) $$
a displacement field (through a chosen generic mechanical behavior and standard mechanical equilibrium) that allows to match the two image intensities, the regularization term allowing to control the smoothness of the traction field. In optimization, such constraints can be enforced strongly through Lagrange multipliers, or approximately through penalization [31]. Here we propose to simply use penalization, such that the problem becomes:

$$ \text{Find } \underline{T}^{\text{sol}}, \underline{U}^{\text{sol}} := \operatorname*{argmin}_{\underline{T}, \underline{U}} \left\{ \begin{aligned} & J(\underline{T}, \underline{U}) := (1-\beta)\|\underline{I} \circ \underline{\Phi} - \underline{I}_0\|_{\underline{\Omega}_0}^2 \\ & + \beta \|\underline{T}\|_{\underline{\partial}\underline{\Omega}_0}^2 + \gamma \left( \|\underline{\mathrm{Div}}(\underline{P})\|_{\underline{\Omega}_0}^2 + \|{}^t\underline{P} \cdot {}^t\underline{F}^{-1} - \underline{F}^{-1} \cdot \underline{P}\|_{\underline{\Omega}_0}^2 + \|\underline{P} \cdot \underline{N} - \underline{T}\|_{\underline{\partial}\underline{\Omega}_0}^2 \right) \end{aligned} \right\}, \quad (26) $$

where $\gamma$ is the penalization coefficient, while the three additional terms represent the balance of linear momentum (without imposed force, though it could be introduced straightforwardly), the balance of angular momentum (which is usually verified exactly thanks to the constitutive framework), and the balance with applied tractions, respectively. In this formulation, $\underline{\Phi}$, $\underline{F}$ and $\underline{P}$ are the mapping, deformation gradient and first Piola-Kirchhoff stress tensor associated to the displacement field $\underline{U}$. Again, the norms will be specified later on. To simplify the formulation, we can remove the traction field variable by enforcing strongly the last balance term, i.e., $\underline{T} = \underline{P}(\underline{U}) \cdot \underline{N}$, leading to the following formulation:

$$ \text{Find } \underline{U}^{\text{sol}} := \operatorname*{argmin}_{\underline{U}} \left\{ J(\underline{U}) := (1-\beta)\|\underline{I} \circ \underline{\Phi} - \underline{I}_0\|_{\underline{\Omega}_0}^2 + \beta\|\underline{P} \cdot \underline{N}\|_{\underline{\partial}\underline{\Omega}_0}^2 + \gamma\|\underline{\mathrm{Div}}(\underline{P})\|_{\underline{\Omega}_0}^2 \right\}. \quad (27) $$

This formulation corresponds formally (i.e., using the right norms) to the original formulation (3), when using $J_b^{\text{reg,eq}}$ for the bulk regularization and $J_{n/t}^{\text{reg,eq}}$ for the surface terms. This gives another point of view on each term, i.e., that the bulk term can be seen as a penalization term, while the actual regularization comes from the surface terms. Also, it is important to notice that this reformulation only makes sense for “equilibrium gap” regularization terms: if one would consider the motion tracking problem (3) with elastic regularization term (7) as an inverse problem, the associated direct problem would not have a proper physical meaning, as it would correspond to a body deforming under a body force induced by the image mismatch and no boundary force.

## 2.5. Numerical resolution

After describing existing and new regularization terms, we recall that the general motion tracking problem formulation was given by Equation (5), with the image correlation term $J^{\text{ima}}$ given by Equation (6), and where the regularization term $J^{\text{reg}}$ may take the form $J^{\text{reg,el}}$ (Equation (7)), $J^{\text{reg,eq,cont}}$ (Equation (9)), $J_b^{\text{reg,eq}}$ (Equation (12)), $J_n^{\text{reg,eq}}$ (Equation (18)), $J_t^{\text{reg,eq}}$ (Equation (22) in 2D, Equation (24) in 3D), or combinations of such terms. Indeed, as already discussed, it is often necessary to combine bulk and boundary terms. In such cases one might want to use different weights for each regularization term [19]; however, to simplify the analysis of the performance of the various terms considered here, we propose to simply sum them, such that they have the same weight, and there is only one parameter controlling the regularization strength, namely $\beta$.

Minimization problem (5) is actually formulated as a root finding problem:

$$ \underline{U} | \underline{\nabla} \underline{J}(\underline{U}) = (1 - \beta) \underline{\nabla} \underline{J}^{\text{ima}}(\underline{U}) + \beta \underline{\nabla} \underline{J}^{\text{reg}}(\underline{U}) = 0 \quad (28) $$

where the gradient of the image term is simply:

$$ \underline{\nabla} \underline{J}^{\text{ima}}(\underline{U}) := \int_{\underline{\Omega}_0} \left( I \left( \underline{X} + {}^t\underline{N} \cdot \underline{U} \right) - \underline{I}_0 \right) \frac{\partial I}{\partial \underline{x}} \left( \underline{X} + {}^t\underline{N} \cdot \underline{U} \right) \cdot \underline{N} \, d\underline{\Omega}_0, \quad (29) $$

and the gradients of the proposed regularization terms are:

$$ \underline{\nabla}_{\underline{b}}^{\text{reg,eq}}(\underline{U}) := {}^t d\underline{\mathbb{R}}_{\underline{b}}(\underline{U}) \cdot \underline{M}^{-1} \cdot \underline{\mathbb{R}}_{\underline{b}}(\underline{U}) \quad (30) $$
with

$$ (d\mathbb{R}_b)_{ij} := \begin{cases} \iint_{\Omega_0} \underline{\text{Grad}}(\underline{\mathcal{N}_i}) : \underline{\frac{\partial P}{\partial F}} : \underline{\text{Grad}}(\underline{\mathcal{N}_j}) d\Omega_0 & \text{if } i \text{ body d.o.f.} \\ 0 & \text{if } i \text{ boundary d.o.f.,} \end{cases} \quad (31) $$

and

$$ \underline{\nabla \mathcal{J}^{\text{reg,eq}}}_{n/t}(\underline{\mathcal{U}}) := t \underline{d \mathbb{R}_{n/t}}(\underline{\mathcal{U}}) \cdot \underline{M_{\partial \Omega_0}^{-1}} \cdot \underline{\mathbb{R}_{n/t}}(\underline{\mathcal{U}}) \quad (32) $$

with

$$ (d\mathbb{R}_{n/t})_{ij} := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \iint_{\partial\Omega_0} \text{Div}_s(\underline{\mathcal{N}_i}) \frac{\partial F_{n/t}}{\partial \underline{\underline{F}}} : \underline{\text{Grad}}(\underline{\mathcal{N}_j}) d\partial\Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (33) $$

We propose to solve this problem using a Gauss-Newton method, which was shown to have good convergence properties for the image correlation problem in the thorough compared analysis of the performance of multiple variants of the method proposed in [32]. Details were given in [17] in the context of the regularization term $J^{\text{reg,eq,cont}}$ (Equation (9)). The key point is that the Jacobian associated to the image term only contains the image gradient product term, not the image hessian term which usually degrades the convergence due to the double derivative of the noise; in principle the second derivation of the regularization terms could be computed exactly, however it is tedious and not necessary, so we employ the same approach as for the image term. Thus, at each Newton iteration we solve the following linear system:

$$ \underline{\Delta \mathcal{U}} | \left( (1 - \beta) \underline{\nabla \nabla \mathcal{J}^{\text{ima}}} + \beta \underline{\nabla \nabla \mathcal{J}^{\text{reg}}} \right) \underline{\Delta \mathcal{U}} = - \left( (1 - \beta) \underline{\nabla \mathcal{J}^{\text{ima}}} + \beta \underline{\nabla \mathcal{J}^{\text{reg}}} \right), \quad (34) $$

where

$$ \underline{\nabla \nabla \mathcal{J}^{\text{ima}}} := \int_{\Omega_0} \underline{\mathcal{N}} \cdot \left( \frac{\partial I}{\partial \underline{x}} \cdot \frac{t \partial I}{\partial \underline{x}} \right) \cdot t \underline{\mathcal{N}} d\Omega_0, \quad (35) $$

which corresponds to a mass matrix weighted by the deformed image gradients, and

$$ \underline{\nabla \nabla \mathcal{J}^{\text{reg,eq}}}_{b/n/t} := t \underline{d \mathbb{R}_{b/n/t}} \cdot \underline{M_{\Omega_0/\partial\Omega_0}^{-1}} \cdot \underline{d \mathbb{R}_{b/n/t}} \quad (36) $$

As a stopping criterion for the Newton iterations, we propose to simply use the relative displacement with a prescribed tolerance, which is well adapted for the motion tracking problem:

$$ \text{err} := \frac{\|\underline{\Delta \mathcal{U}}\|}{\|\underline{\mathcal{U}}\|} < \text{tol.} \quad (37) $$

Nevertheless, in practice, we found that such Gauss-Newton iterations do not always converge toward a solution for such a highly stiff problem. Indeed, they often lead to inverted elements, for which the mechanical model cannot be evaluated. This is a situation the regularization is supposed to prevent through the use of a proper energy potential with an infinite energy barrier (for instance, the Ogden-Ciarlet-Geymonat has a $-\ln(J)$ term, where $J$ denotes the volume change [33, 34]), but the Newton iterations can sometimes pass this energy barrier. Indeed, it is possible, and it is actually the case for the Ogden-Ciarlet-Geymonat potential, that even though the potential is only defined for non-inverted elements, the expression of its derivative is well defined for inverted elements. This is a common issue in large deformation computation, which is usually solved using adaptation time stepping [35, 36]. However, it is not possible here because the “time” increment is controlled by the image temporal discretization. Hence, we augment the Newton iterations with a backtracking line search [37] that prevents the passing of energy barriers.

Finally, the image series integrator is described in Algorithm 1, the nonlinear solver in Algorithm 2, and the line search in Algorithm 3.

The richness of the displacement solution space is controlled by the finite element mesh and interpolation degree. When using a coarse mesh, with many pixels per element, the mesh
with

$$ (d\mathbb{R}_b)_{ij} := \begin{cases} \iint_{\Omega_0} \underline{\text{Grad}}(\underline{\mathcal{N}_i}) : \underline{\frac{\partial P}{\partial F}} : \underline{\text{Grad}}(\underline{\mathcal{N}_j}) d\Omega_0 & \text{if } i \text{ body d.o.f.} \\ 0 & \text{if } i \text{ boundary d.o.f.,} \end{cases} \quad (31) $$

and

$$ \underline{\nabla \mathcal{J}^{\text{reg,eq}}}_{n/t}(\underline{\mathcal{U}}) := t \underline{d\mathbb{R}_{n/t}}(\underline{\mathcal{U}}) \cdot \underline{M_{\partial \Omega_0}^{-1}} \cdot \underline{\mathbb{R}_{n/t}}(\underline{\mathcal{U}}) \quad (32) $$

with

$$ (d\mathbb{R}_{n/t})_{ij} := \begin{cases} 0 & \text{if } i \text{ body d.o.f.} \\ \iint_{\partial\Omega_0} \text{Div}_s(\underline{\mathcal{N}_i}) \frac{\partial F_{n/t}}{\partial \underline{\underline{F}}} : \underline{\text{Grad}}(\underline{\mathcal{N}_j}) d\partial\Omega_0 & \text{if } i \text{ boundary d.o.f.} \end{cases} \quad (33) $$

We propose to solve this problem using a Gauss–Newton method, which was shown to have good convergence properties for the image correlation problem in the thorough compared analysis of the performance of multiple variants of the method proposed in [32]. Details were given in [17] in the context of the regularization term $J^{\text{reg,eq,cont}}$ (Equation (9)). The key point is that the Jacobian associated to the image term only contains the image gradient product term, not the image hessian term which usually degrades the convergence due to the double derivative of the noise; in principle the second derivation of the regularization terms could be computed exactly, however it is tedious and not necessary, so we employ the same approach as for the image term. Thus, at each Newton iteration we solve the following linear system:

$$ \underline{\Delta \mathcal{U}} | \left( (1 - \beta) \underline{\nabla \nabla \mathcal{J}^{\text{ima}}} + \beta \underline{\nabla \nabla \mathcal{J}^{\text{reg}}} \right) \underline{\Delta \mathcal{U}} = - \left( (1 - \beta) \underline{\nabla \mathcal{J}^{\text{ima}}} + \beta \underline{\nabla \mathcal{J}^{\text{reg}}} \right), \quad (34) $$

where

$$ \underline{\nabla \nabla \mathcal{J}^{\text{ima}}} := \int_{\Omega_0} \underline{\mathcal{N}} \cdot \left( \frac{\partial I}{\partial \underline{x}} \cdot \frac{t}{\underline{x}} \right) \cdot \underline{t} \underline{\mathcal{N}} d\Omega_0, \quad (35) $$

which corresponds to a mass matrix weighted by the deformed image gradients, and

$$ \underline{\nabla \nabla \mathcal{J}^{\text{reg,eq}}}_{b/n/t} := t \underline{d\mathbb{R}_{b/n/t}} \cdot \underline{M_{\Omega_0/\partial\Omega_0}^{-1}} \cdot \underline{d\mathbb{R}_{b/n/t}} \quad (36) $$

As a stopping criterion for the Newton iterations, we propose to simply use the relative displacement with a prescribed tolerance, which is well adapted for the motion tracking problem:

$$ \text{err} := \frac{\|\underline{\Delta \mathcal{U}}\|}{\|\underline{\mathcal{U}}\|} < \text{tol.} \quad (37) $$

Nevertheless, in practice, we found that such Gauss–Newton iterations do not always converge toward a solution for such a highly stiff problem. Indeed, they often lead to inverted elements, for which the mechanical model cannot be evaluated. This is a situation the regularization is supposed to prevent through the use of a proper energy potential with an infinite energy barrier (for instance, the Ogden–Ciarlet–Geymonat has a $-\ln(J)$ term, where $J$ denotes the volume change [33, 34]), but the Newton iterations can sometimes pass this energy barrier. Indeed, it is possible, and it is actually the case for the Ogden–Ciarlet–Geymonat potential, that even though the potential is only defined for non-inverted elements, the expression of its derivative is well defined for inverted elements. This is a common issue in large deformation computation, which is usually solved using adaptation time stepping [35, 36]. However, it is not possible here because the “time” increment is controlled by the image temporal discretization. Hence, we augment the Newton iterations with a backtracking line search [37] that prevents the passing of energy barriers.

Finally, the image series integrator is described in Algorithm 1, the nonlinear solver in Algorithm 2, and the line search in Algorithm 3.

The richness of the displacement solution space is controlled by the finite element mesh and interpolation degree. When using a coarse mesh, with many pixels per element, the mesh
**Initialisation**

Read initial image $I_0$

Initiate finite element solution: $\underline{U} \leftarrow 0$

**foreach** frame index $t = 1, 2, \dots$ **do**

Read current image $I_t$

Compute current displacement: $\underline{U} \leftarrow \text{nonlinear\_solver}(I_0, I_t)$

**end**

**Algorithm 1:** Image series integrator. The nonlinear solver is detailed Algorithm 2. Note that in this basic version, the Newton iterations at frame $t$ are naturally initialized with the converged solution at frame $t-1$. A multi-level version is presented in Algorithm 4.

**while** err > tol **do**

Assemble residual, Jacobian: $\underline{\nabla J}, \overline{\nabla \nabla J}$

Compute solution increment: $\underline{\Delta U} \leftarrow \text{linear\_solver}(\underline{\nabla J}, \overline{\nabla \nabla J})$

Update error: $\text{err} \leftarrow \frac{\|\underline{\Delta U}\|}{\|\underline{U}\|}$

Compute relaxation: $\alpha \leftarrow \text{line\_search}(\underline{U}, \underline{\Delta U})$

Update solution: $\underline{U} \leftarrow \underline{U} + \alpha \underline{\Delta U}$

**end**

**Algorithm 2:** Nonlinear solver. In practice we use a tolerance of tol = 1%.

**Initialisation**

Compute initial energy: $J_0 \leftarrow J(\underline{U})$

Initialise relaxation counter: $k \leftarrow 0$

**while** $J > J_0$ **do**

Define current relaxation: $\alpha \leftarrow f^k$

Compute current energy: $J \leftarrow J(\underline{U} + \alpha \underline{\Delta U})$

Update relaxation counter: $k \leftarrow k + 1$

**end**

**Algorithm 3:** Backtracking line search. In practice we use a backtracking factor of $f = 0.5$. Here $J$ refers to the full cost function, encompassing both the image and regularization terms, as defined Problem (5).

itself acts as some kind of kinematical regularization, limiting the size of the displacement solution space. Conversely, when solving for very fine meshes, for instance when elements reach the size of the image pixels, convergence becomes problematic as the displacement solution space becomes huge and there is very little information per element; in this case mechanical regularization compensates for the lack of information, but it might not be enough to obtain a robust convergence of the nonlinear iterations. Thus, multi-resolution is often necessary to perform motion tracking on very fine meshes [9, 19], which consists in performing tracking on successively refined meshes, initiating the nonlinear iterations of a given frame and given refinement level by the converged solution obtained at the same frame but at the previous refinement level, instead of the converged solution obtained at the same refinement level but at the previous frame, as it is naturally the case in the single-refinement-level Algorithm 1. Since two triangulations of the same geometric domain might not overlap, for instance if the domain has a curved boundary, and especially if the triangulations have different characteristic sizes, one cannot simply interpolate the displacement field from the coarse grid to the fine one. Instead we propose to use a projection. The multi-resolution frame integrator is detailed in Algorithm 4.
**Initialisation**

Read initial image $I_0$

Instantiate finite element solution: $\underline{U} \leftarrow 0$

**foreach** frame index $t = 1, 2, \dots$ **do**

Read current image $I_t$

Compute current displacement: $\underline{U} \leftarrow \text{nonlinear\_solver}(I_0, I_t)$

**end**

**Algorithm 1:** Image series integrator. The nonlinear solver is detailed Algorithm 2. Note that in this basic version, the Newton iterations at frame $t$ are naturally initialized with the converged solution at frame $t-1$. A multi-level version is presented in Algorithm 4.

**while** err > tol **do**

Assemble residual, Jacobian: $\underline{\nabla J}, \overline{\nabla \nabla J}$

Compute solution increment: $\underline{\Delta U} \leftarrow \text{linear\_solver}(\underline{\nabla J}, \overline{\nabla \nabla J})$

Update error: $\text{err} \leftarrow \frac{\|\underline{\Delta U}\|}{\|\underline{U}\|}$

Compute relaxation: $\alpha \leftarrow \text{line\_search}(\underline{U}, \underline{\Delta U})$

Update solution: $\underline{U} \leftarrow \underline{U} + \alpha \underline{\Delta U}$

**end**

**Algorithm 2:** Nonlinear solver. In practice we use a tolerance of tol = 1%.

**Initialisation**

Compute initial energy: $J_0 \leftarrow J(\underline{U})$

Initialise relaxation counter: $k \leftarrow 0$

**while** $J > J_0$ **do**

Define current relaxation: $\alpha \leftarrow f^k$

Compute current energy: $J \leftarrow J(\underline{U} + \alpha \underline{\Delta U})$

Update relaxation counter: $k \leftarrow k + 1$

**end**

**Algorithm 3:** Backtracking line search. In practice we use a backtracking factor of $f = 0.5$. Here $J$ refers to the full cost function, encompassing both the image and regularization terms, as defined Problem (5).

itself acts as some kind of kinematical regularization, limiting the size of the displacement solution space. Conversely, when solving for very fine meshes, for instance when elements reach the size of the image pixels, convergence becomes problematic as the displacement solution space becomes huge and there is very little information per element; in this case mechanical regularization compensates for the lack of information, but it might not be enough to obtain a robust convergence of the nonlinear iterations. Thus, multi-resolution is often necessary to perform motion tracking on very fine meshes [9, 19], which consists in performing tracking on successively refined meshes, initiating the nonlinear iterations of a given frame and given refinement level by the converged solution obtained at the same frame but at the previous refinement level, instead of the converged solution obtained at the same refinement level but at the previous frame, as it is naturally the case in the single-refinement-level Algorithm 1. Since two triangulations of the same geometric domain might not overlap, for instance if the domain has a curved boundary, and especially if the triangulations have different characteristic sizes, one cannot simply interpolate the displacement field from the coarse grid to the fine one. Instead we propose to use a projection. The multi-resolution frame integrator is detailed in Algorithm 4.
foreach refinement level $k = 1, 2, \dots$ do

Initialisation

Read initial image $I_0$

Instantiate finite element solution $\mathcal{U}^k$

foreach frame index $t = 1, 2, \dots$ do

Read current image $I_t$

if $k > 1$ then

Initialize displacement: $\mathcal{U}^k \leftarrow \mathcal{U}_t^{k-1}$

Compute current displacement: $\mathcal{U}^k \leftarrow \text{nonlinear solver}(I_0, I_t)$

Save displacement: $\mathcal{U}_t^k \leftarrow \mathcal{U}^k$

end

end

**Algorithm 4:** Image series integrator with multi-resolution. Compared to the single-refinement-level Algorithm 1, here the Newton iterations at frame $t$ of refinement level $k$ are initialized with the converged solution at frame $t$ of refinement level $k-1$.

In order to speed up the multi-resolution algorithm, pyramidal filtering can be applied to the images, such that the tracking of successively refined meshes is performed on successively more resolved images, until the final mesh is used to track the original high-resolution images [38]; however, for the sake of simplicity, here we use the full images at all mesh refinement levels.

These algorithms have been implemented in an open-source library [39] written in python and based on the FEniCS [40, 41] and VTK [42] libraries. It is currently freely available online at https://gitlab.inria.fr/mgenet/dolfin_warp. We also provide the code to reproduce the results of this paper under the form of jupyter notebooks [43]: static versions are given in the Supplementary material of the paper while interactive versions are currently available online at https://mgenet.gitlabpages.inria.fr/N-DEG-paper-demos/index.html.

## 2.6. Synthetic data

In order to establish the tracking performance of the proposed method, we generated synthetic images corresponding to various objects (simple square, cardiac-like ring) and motions (rigid translation and rotation, homogeneous compression and shear, cardiac-like contraction and twist) with various noise levels. To focus on the regularization term itself, we limited ourselves to highly resolved (though the impact of image resolution could be investigated as well [44]), and textured (with used a tagged-MRI-like pattern, see [45]) images. The images occupy the spatial domain $0; 1^2$ (arbitrary unit), the temporal domain $0; 1$ (arbitrary unit), are discretized with $100 \times 100$ pixels spatially and 21 frames temporally.

For the simple (analytical) motion of the square, we define the initial domain as $\Omega_0^t := 0.1; 0.7 \times 0.2; 0.8$ for the translation case and $\Omega_0^{r,c,s} := 0.2; 0.8^2$ for the rotation, compression and shear cases. The motion models are given by

$$ \left\{ \begin{array}{ll} \underline{x}^t(\underline{X}, t) := \underline{X} + t \underline{D} & \text{with } \underline{D} := \begin{pmatrix} 0.2 \\ 0 \end{pmatrix} \\ \underline{x}^r(\underline{X}, t) := \underline{X}_0 + \underline{R} \cdot (\underline{X} - \underline{X}_0) & \text{with } \underline{X}_0 := \begin{pmatrix} 0.5 \\ 0.5 \end{pmatrix}, \underline{R} := \begin{pmatrix} +\cos(\theta) & -\sin(\theta) \\ +\sin(\theta) & +\cos(\theta) \end{pmatrix}, \theta = \frac{t\pi}{4} \\ \underline{x}^c(\underline{X}, t) := \underline{X}_0 + \underline{F}^c \cdot (\underline{X} - \underline{X}_0) & \text{with } \underline{X}_0 := \begin{pmatrix} 0.5 \\ 0.5 \end{pmatrix}, \underline{F}^c := \begin{pmatrix} \sqrt{1-2t} & 0 \\ 0 & 1 \end{pmatrix} \\ \underline{x}^s(\underline{X}, t) := \underline{X}_0 + \underline{F}^s \cdot (\underline{X} - \underline{X}_0) & \text{with } \underline{X}_0 := \begin{pmatrix} 0.5 \\ 0.5 \end{pmatrix}, \underline{F}^s := \begin{pmatrix} 1 & 0.2 \\ 0 & 1 \end{pmatrix}, \end{array} \right. \quad (38) $$

which represent a 0.2 (arbitrary unit) translation, a $\pi/4$ rad rotation, a 20% compression and a 20% shear, respectively.
Then, for each case, the image generation consisted in the following steps: for each frame (time t), for each pixel (position x), the initial (i.e., first time frame) position of the pixel (X(x, t)) was determined through inversion of the mappings (38) and the intensity was given by

$$I(\underline{x}, t) := I_0(\underline{X}(\underline{x}, t)) \quad \text{with} \quad I_0(\underline{X}) := \sqrt{\left|\sin\left(\frac{\pi\underline{X}}{s}\right)\right|\left|\sin\left(\frac{\pi\underline{Y}}{s}\right)\right|} \quad (39)$$

Note that a more complex imaging model [44] could have been considered as well.

For the cardiac-like case, a ring (center $X_0 := (0.5, 0.5)$, internal radius 0.2, external radius 0.4 (arbitrary units)) was defined, and an hyperelastic (neo-Hookean and Ogden-Ciarlet-Geymonat potentials with unit Young modulus and 0.3 Poisson coefficient) finite element model with very fine mesh (element size equal to the image pixel size) was run with prescribed displacement applied to the internal (inward displacement of 0.1 and rotation of $-\pi/4$ rad) and external (inward displacement of 0.05 and rotation of $-\pi/8$ rad) edges, mimicking the in-plane motion of a cardiac slice. And to generate the images, we employed the following approach: for each time frame (time t) the reference mesh was warped by applying the computed displacement field, and the displacement field was projected onto the image, such that for each pixel (position x) the reference position of the pixel was computed as $X(x, t) = x - U(x, t)$. Finally, the texture model (39) was applied.

Noise was eventually added to the images, characterized by the signal-to-noise ratio (SNR). The magnitude of the signal here is 1, and we added random Gaussian noise with zero mean and standard deviation of 0.1, 0.2 and 0.3, corresponding to SNR of 10, 5 and 3.3, respectively. Note that, like for the imaging model, more complex noise models, for instance including spatial and/or temporal correlations [44], could have been considered.

## 3. Results and discussion

To establish the tracking performance of the proposed method (described in Section 2.3), we now present tracking results on various synthetic images (described in Section 2.6), for various regularization terms, namely elastic (Equation (7), continuous version of the equilibrium gap (Equation (9)), discrete version of the equilibrium gap (Equation (12)) and discrete version of the nonlinear equilibrium gap including surface traction regularization terms (Equations (12), (18), (22) and/or (24))). For each regularization term, we considered both a small strain approximation with a Hooke strain energy potential, and a large strain formulation with the neo-Hookean and Ogden-Ciarlet-Geymonat potentials. For all models, we considered a unit Young modulus and null Poisson ratio.

The tracking performance is evaluated in terms of a normalized tracking error defined as

$$\text{err} := \frac{\sqrt{\frac{1}{T} \int_{0}^{T} \frac{1}{|\Omega_{0}|} \int_{\Omega_{0}} \|\underline{U} - \underline{U}^{\text{ex}}\|^{2}}}{\sqrt{\frac{1}{T} \int_{0}^{T} \frac{1}{|\Omega_{0}|} \int_{\Omega_{0}} \|\underline{U}^{\text{ex}}\|^{2}}}, \quad (40)$$

where $\underline{U}$ is the tracked displacement and $\underline{U}^{\text{ex}}$ is the ground truth used to generate the images.

### 3.1. Rigid body motion

We start with simple rigid motions, namely pure translation and pure rotation. Figures 1 and 2 show normalized tracking error as a function of regularization strength $\beta \in 0;1$, for various levels of image SNR and for various regularization terms. In the plots, dots represent noise realizations, and lines represent their average.
Then, for each case, the image generation consisted in the following steps: for each frame (time t), for each pixel (position x), the initial (i.e., first time frame) position of the pixel (X(x, t)) was determined through inversion of the mappings (38) and the intensity was given by

$$I(\underline{x}, t) := I_0(\underline{X}(\underline{x}, t)) \quad \text{with} \quad I_0(\underline{X}) := \sqrt{\left|\sin\left(\frac{\pi X}{s}\right)\right|\left|\sin\left(\frac{\pi Y}{s}\right)\right|} \qquad (39)$$

Note that a more complex imaging model [44] could have been considered as well.

For the cardiac-like case, a ring (center $X_0 := (0.5, 0.5)$, internal radius 0.2, external radius 0.4 (arbitrary units)) was defined, and an hyperelastic (neo-Hookean and Ogden-Ciarlet-Geymonat potentials with unit Young modulus and 0.3 Poisson coefficient) finite element model with very fine mesh (element size equal to the image pixel size) was run with prescribed displacement applied to the internal (inward displacement of 0.1 and rotation of $-\pi/4$ rad) and external (inward displacement of 0.05 and rotation of $-\pi/8$ rad) edges, mimicking the in-plane motion of a cardiac slice. And to generate the images, we employed the following approach: for each time frame (time t) the reference mesh was warped by applying the computed displacement field, and the displacement field was projected onto the image, such that for each pixel (position x) the reference position of the pixel was computed as $X(x, t) = x - U(x, t)$. Finally, the texture model (39) was applied.

Noise was eventually added to the images, characterized by the signal-to-noise ratio (SNR). The magnitude of the signal here is 1, and we added random Gaussian noise with zero mean and standard deviation of 0.1, 0.2 and 0.3, corresponding to SNR of 10, 5 and 3.3, respectively. Note that, like for the imaging model, more complex noise models, for instance including spatial and/or temporal correlations [44], could have been considered.

## 3. Results and discussion

To establish the tracking performance of the proposed method (described in Section 2.3), we now present tracking results on various synthetic images (described in Section 2.6), for various regularization terms, namely elastic (Equation (7), continuous version of the equilibrium gap (Equation (9)), discrete version of the equilibrium gap (Equation (12)) and discrete version of the nonlinear equilibrium gap including surface traction regularization terms (Equations (12), (18), (22) and/or (24))). For each regularization term, we considered both a small strain approximation with a Hooke strain energy potential, and a large strain formulation with the neo-Hookean and Ogden-Ciarlet-Geymonat potentials. For all models, we considered a unit Young modulus and null Poisson ratio.

The tracking performance is evaluated in terms of a normalized tracking error defined as

$$\text{err} := \frac{\sqrt{\frac{1}{T} \int_{0}^{T} \frac{1}{|\Omega_0|} \int_{\Omega_0} \|\underline{U} - \underline{U}^{\text{ex}}\|^2}}{\sqrt{\frac{1}{T} \int_{0}^{T} \frac{1}{|\Omega_0|} \int_{\Omega_0} \|\underline{U}^{\text{ex}}\|^2}}, \qquad (40)$$

where $\underline{U}$ is the tracked displacement and $\underline{U}^{\text{ex}}$ is the ground truth used to generate the images.

### 3.1. Rigid body motion

We start with simple rigid motions, namely pure translation and pure rotation. Figures 1 and 2 show normalized tracking error as a function of regularization strength $\beta \in 0;1$, for various levels of image SNR and for various regularization terms. In the plots, dots represent noise realizations, and lines represent their average.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Translation case. (a) Normalized displacement error (40) as a function of regularization strength ($\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12), $J_n^{\text{reg,eq}}$ (18) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). For such a simple motion, basically all regularization terms allow to reduce the tracking error. (Continued on next page.)</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Translation case. (a) Normalized displacement error (40) as a function of regularization strength ($\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12), $J_n^{\text{reg,eq}}$ (18) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). For such a simple motion, basically all regularization terms allow to reduce the tracking error. (Continued on next page.)</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 1 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (18) and (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels.</figcaption></figure>

Let us first discuss the pure translation case, i.e., Figure 1. For noiseless images, the tracking error is below 0.1% for all regularization terms and all regularization strengths. Increasing noise level leads to an increase of both the tracking error mean and dispersion. For all considered regularization terms, increasing the regularization strength decreases the mean error, though the dispersion is little impacted. Basically, for such a simple motion, all regularization terms perform rather well. This is explained by the fact that a rigid translation cancels all considered regularization terms, allowing to filter noise-induced spurious motions without interfering with the tracking itself.

One can see that elastic (first row in Figure 1a) and continuous equilibrium gap (second row in Figure 1a) terms behave better than the discrete equilibrium gap terms (third and fourth rows in Figure 1a). This is due to the fact that these terms basically prevent the mesh from deforming, or from deforming in a non affine way, which is compatible with the exact solution here, hence the good tracking performance. However, as we will see later, this constraint will prove problematic for more complex motions.

Focusing on the regularization term introduced in this paper (third and fourth rows in Figure 1a), one can see that the error dispersion becomes significant for images with low SNR. This is due to the fact that these terms represent a much lighter constraint on the displacement field,
<figure><img src="image_4.png" /><figcaption>**Figure 1 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (18) and (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels.</figcaption></figure>

Let us first discuss the pure translation case, i.e., Figure 1. For noiseless images, the tracking error is below 0.1% for all regularization terms and all regularization strengths. Increasing noise level leads to an increase of both the tracking error mean and dispersion. For all considered regularization terms, increasing the regularization strength decreases the mean error, though the dispersion is little impacted. Basically, for such a simple motion, all regularization terms perform rather well. This is explained by the fact that a rigid translation cancels all considered regularization terms, allowing to filter noise-induced spurious motions without interfering with the tracking itself.

One can see that elastic (first row in Figure 1a) and continuous equilibrium gap (second row in Figure 1a) terms behave better than the discrete equilibrium gap terms (third and fourth rows in Figure 1a). This is due to the fact that these terms basically prevent the mesh from deforming, or from deforming in a non affine way, which is compatible with the exact solution here, hence the good tracking performance. However, as we will see later, this constraint will prove problematic for more complex motions.

Focusing on the regularization term introduced in this paper (third and fourth rows in Figure 1a), one can see that the error dispersion becomes significant for images with low SNR. This is due to the fact that these terms represent a much lighter constraint on the displacement field,
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Rotation case. (a) Normalized displacement error (40) as a function of regularization strength ($\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12), $J_n^{\text{reg,eq}}$ (18) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The linear elastic regularization term interferes with the tracking because finite rotation lead to nonzero infinitesimal strain, thus generating spurious elastic energy. (Continued on next page.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Rotation case. (a) Normalized displacement error (40) as a function of regularization strength ($\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12), $J_n^{\text{reg,eq}}$ (18) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The linear elastic regularization term interferes with the tracking because finite rotation lead to nonzero infinitesimal strain, thus generating spurious elastic energy. (Continued on next page.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (18) and (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

solely enforcing that it is close to an equilibrium solution; this is, however, the very reason why it performs well on basically any motion as we will see later. Nevertheless, one can also see that adding the boundary traction terms (fourth row in Figure 1a), which penalizes the non smoothness of the normal and tangential tractions at the edges, helps decreasing the tracking error mean and dispersion compared to bulk terms only (third row in Figure 1a).

Let us now discuss the pure rotation case, i.e., Figure 2. For the nonlinear elastic (upper right in Figure 2a) and all equilibrated (second, third and fourth rows in Figure 2a) regularization terms, the conclusions are the same as for the pure translation case. For the linear elastic regularization term (upper left in Figure 2a), however, the result is quite different. Indeed, if increasing slightly the regularization strength allows to decrease the tracking error, increasing it further completely degrades the tracking. This is due to the fact that finite rotations generate nonzero infinitesimal strains, thus nonzero elastic energy; hence the algorithm, which tries to minimize this energy, tends to underestimate rotations (see Figure 2c). This is not the case for hyperelastic regularization (upper right in Figure 2a), as in the finite strain setting, rigid rotations do not generate any strain. Interestingly, this is also not the case for other regularization terms with linearized kinematics and behaviors (left column, second to fourth rows in Figure 2a),
<figure><img src="image_6.png" /><figcaption>**Figure 2 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (18) and (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

solely enforcing that it is close to an equilibrium solution; this is, however, the very reason why it performs well on basically any motion as we will see later. Nevertheless, one can also see that adding the boundary traction terms (fourth row in Figure 1a), which penalizes the non smoothness of the normal and tangential tractions at the edges, helps decreasing the tracking error mean and dispersion compared to bulk terms only (third row in Figure 1a).

Let us now discuss the pure rotation case, i.e., Figure 2. For the nonlinear elastic (upper right in Figure 2a) and all equilibrated (second, third and fourth rows in Figure 2a) regularization terms, the conclusions are the same as for the pure translation case. For the linear elastic regularization term (upper left in Figure 2a), however, the result is quite different. Indeed, if increasing slightly the regularization strength allows to decrease the tracking error, increasing it further completely degrades the tracking. This is due to the fact that finite rotations generate nonzero infinitesimal strains, thus nonzero elastic energy; hence the algorithm, which tries to minimize this energy, tends to underestimate rotations (see Figure 2c). This is not the case for hyperelastic regularization (upper right in Figure 2a), as in the finite strain setting, rigid rotations do not generate any strain. Interestingly, this is also not the case for other regularization terms with linearized kinematics and behaviors (left column, second to fourth rows in Figure 2a),
<figure><img src="image_6.png" /><figcaption>**Figure 2 (cont.).** (c) Tracking solutions for the linear elastic regularization term (7), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The linear elastic regularization term tends to reduce rotations because they generate spurious elastic energy, which degrades the tracking.</figcaption></figure>

because here the spurious strain induced by the finite rotation is homogeneous, hence the associated stress is also homogeneous and thus equilibrated, i.e., it does not generate any spurious equilibrium gap.

## 3.2. Nonrigid homogeneous deformation

We now consider nonrigid—though still homogeneous—deformations, namely compression (Figure 3) and shear (Figure 4). Conclusions are the same as for rigid transformations for all regularization terms based on the equilibrium gap principle (i.e., second, third and fourth rows in Figures 3a and 4a), which perform very well. Regularization terms based on the elastic energy (i.e., first row in Figures 3a and 4a) are, on the contrary, very problematic as they tend to prevent the mesh from deforming. This illustrates very well the fact that such regularization terms should not be used when tracking nonrigid deformations, as they interfere pathologically with the tracking.

One can also notice that the continuous formulation of the equilibrium gap regularization (second row in Figures 3a and 4a) seems to perform better than the discrete formulation (third and fourth rows in Figures 3a and 4a); as we will see later, this is actually due to the fact that
<figure><img src="image_7.png" /><figcaption>**Figure 2 (cont.).** (c) Tracking solutions for the linear elastic regularization term (7), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The linear elastic regularization term tends to reduce rotations because they generate spurious elastic energy, which degrades the tracking.</figcaption></figure>

because here the spurious strain induced by the finite rotation is homogeneous, hence the associated stress is also homogeneous and thus equilibrated, i.e., it does not generate any spurious equilibrium gap.

## 3.2. Nonrigid homogeneous deformation

We now consider nonrigid—though still homogeneous—deformations, namely compression (Figure 3) and shear (Figure 4). Conclusions are the same as for rigid transformations for all regularization terms based on the equilibrium gap principle (i.e., second, third and fourth rows in Figures 3a and 4a), which perform very well. Regularization terms based on the elastic energy (i.e., first row in Figures 3a and 4a) are, on the contrary, very problematic as they tend to prevent the mesh from deforming. This illustrates very well the fact that such regularization terms should not be used when tracking nonrigid deformations, as they interfere pathologically with the tracking.

One can also notice that the continuous formulation of the equilibrium gap regularization (second row in Figures 3a and 4a) seems to perform better than the discrete formulation (third and fourth rows in Figures 3a and 4a); as we will see later, this is actually due to the fact that
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Compression case. (a) Normalized displacement error (40) as a function of regularization strength (parameter $\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The elastic regularization terms interfere with the tracking because they basically prevent the mesh from deforming. (Continued on next page.)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Compression case. (a) Normalized displacement error (40) as a function of regularization strength (parameter $\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The elastic regularization terms interfere with the tracking because they basically prevent the mesh from deforming. (Continued on next page.)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 3 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

$J^{\text{reg,eq,cont}}$ (9) contains the equilibrium gap induced by the discretization in addition to the equilibrium gap induced by the images, and tends to minimize it, i.e., to maintain the stress as homogeneous as possible, which helps the tracking here because the exact solution has indeed an homogeneous stress field, but which will prove highly pathological for non homogeneous cases.

It is important to notice also that for the compression case, only the tangential term (22) was included in the boundary traction terms (fourth row in Figure 3a). Indeed, the exact solution has highly non-smooth tractions (which are nonzero on left and right edges, zero on bottom and top edges), so that penalizing the surface gradient of the normal surface traction would not make sense. Similarly, in the shear case, only the normal term (18) was included in the boundary tractions terms (fourth row in Figure 4a).

### 3.3. *Nonrigid heterogeneous deformation*

Let us now consider more realistic images, mimicking cardiac tagged magnetic resonance imaging slices. Tracking results are presented in Figure 5. Elastic regularization (first row of Figure 5a)
<figure><img src="image_9.png" /><figcaption>**Figure 3 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

$J^{\text{reg,eq,cont}}$ (9) contains the equilibrium gap induced by the discretization in addition to the equilibrium gap induced by the images, and tends to minimize it, i.e., to maintain the stress as homogeneous as possible, which helps the tracking here because the exact solution has indeed an homogeneous stress field, but which will prove highly pathological for non homogeneous cases.

It is important to notice also that for the compression case, only the tangential term (22) was included in the boundary traction terms (fourth row in Figure 3a). Indeed, the exact solution has highly non-smooth tractions (which are nonzero on left and right edges, zero on bottom and top edges), so that penalizing the surface gradient of the normal surface traction would not make sense. Similarly, in the shear case, only the normal term (18) was included in the boundary tractions terms (fourth row in Figure 4a).

### 3.3. *Nonrigid heterogeneous deformation*

Let us now consider more realistic images, mimicking cardiac tagged magnetic resonance imaging slices. Tracking results are presented in Figure 5. Elastic regularization (first row of Figure 5a)
<figure><img src="image_9.png" /><figcaption>**Figure 3 (cont.).** (c) Tracking solutions for the elastic regularization term (7), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The elastic regularization terms tend to prevent the mesh from deforming, which degrades the tracking.</figcaption></figure>

has the same limitation as for homogeneous deformations, in that it tends to prevent the mesh from deforming and thus interferes with the tracking. Interestingly, the continuous formulation of the equilibrium gap principle (second row of Figure 5a) presents similar limitations, though somewhat weaker for large regularization strengths. This can be explained by the fact that this “energy” contains the equilibrium gap induced by the finite element discretization, thus it tends to minimize the discretization error, i.e., it tends to force the mesh to deform homogeneously, which was compatible with the ground truth of the simple examples of Sections 3.1 and 3.2, but not in more realistic cases. The discrete version of the equilibrium gap principle (third and fourth rows of Figure 5a), on the other hand, which characterizes the equilibrium gap induced by the motion itself but not the one induced by the discretization thanks to the projection step, does not have such pathological behavior, and allows the reduction of the tracking error basically for all regularization levels. One can see that the nonlinear version (based on the finite strain framework and neoHookean and Ogden-Ciarlet-Geymonat potentials) performs slightly better than the linear version (based on the infinitesimal strain framework ant the Hooke law). Also, the surface traction regularization terms (fourth row in Figure 5a) help reducing the tracking error mean and dispersion compared to the tracking with the bulk term only (third row in Figure 5a).
<figure><img src="image_10.png" /><figcaption>**Figure 3 (cont.).** (c) Tracking solutions for the elastic regularization term (7), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The elastic regularization terms tend to prevent the mesh from deforming, which degrades the tracking.</figcaption></figure>

has the same limitation as for homogeneous deformations, in that it tends to prevent the mesh from deforming and thus interferes with the tracking. Interestingly, the continuous formulation of the equilibrium gap principle (second row of Figure 5a) presents similar limitations, though somewhat weaker for large regularization strengths. This can be explained by the fact that this “energy” contains the equilibrium gap induced by the finite element discretization, thus it tends to minimize the discretization error, i.e., it tends to force the mesh to deform homogeneously, which was compatible with the ground truth of the simple examples of Sections 3.1 and 3.2, but not in more realistic cases. The discrete version of the equilibrium gap principle (third and fourth rows of Figure 5a), on the other hand, which characterizes the equilibrium gap induced by the motion itself but not the one induced by the discretization thanks to the projection step, does not have such pathological behavior, and allows the reduce the tracking error basically for all regularization levels. One can see that the nonlinear version (based on the finite strain framework and neoHookean and Ogden-Ciarlet-Geymonat potentials) performs slightly better than the linear version (based on the infinitesimal strain framework ant the Hooke law). Also, the surface traction regularization terms (fourth row in Figure 5a) help reducing the tracking error mean and dispersion compared to the tracking with the bulk term only (third row in Figure 5a).
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Shear case. (a) Normalized displacement error (40) as a function of regularization strength ($\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12) and $J_n^{\text{reg,eq}}$ (18)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The elastic regularization terms interfere with the tracking because they basically prevent the mesh from deforming. (Continued on next page.)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Shear case. (a) Normalized displacement error (40) as a function of regularization strength ($\beta$), for various levels of image noise, for various regularization terms (first row: $J^{\text{reg,el}}$ (7); second row: $J^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12) and $J_n^{\text{reg,eq}}$ (18)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The elastic regularization terms interfere with the tracking because they basically prevent the mesh from deforming. (Continued on next page.)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 4 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (18), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

### 3.4. *Impact of mesh size*

We now investigate in details the impact of mesh size on tracking. Indeed, there is a fundamental trade-off between the richness of the motion model and the robustness of the tracking—the finite element discretization of the tracking already represents some kind of “kinematical” regularization, in the sense that the mesh controls the size of the finite dimensional functional space of the displacement field. Figure 6 (cont.) shows tracking results for the cardiac-like problem, for various mesh sizes (from 0.1, which corresponds to the characteristic size of the image texture, like in the previous examples, to $0.1/2^4 = 0.00625$, which is smaller than the image pixel), for various levels of noise (0 and 0.1) and regularization strengths (0 and 0.1), and for two different mesh refinement strategies: either we simply track the images with the different meshes independently (left in Figure 6 (cont.)a), or we track with successively refined meshes while initializing the tracking at a given mesh refinement level with the converged solution of the previous mesh refinement level (right in Figure 6 (cont.)a). First, one can see that directly tracking images with fine meshes is impossible—mechanical regularization helps a little in some cases, but not enough to obtain satisfying solutions. Thus, to obtain fine solutions, one needs to perform multi-resolution. However, one can also see that multi-resolution is not enough, and mechanical regularization is necessary to obtain fine solutions.
<figure><img src="image_12.png" /><figcaption>**Figure 4 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization (18), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

### 3.4. *Impact of mesh size*

We now investigate in details the impact of mesh size on tracking. Indeed, there is a fundamental trade-off between the richness of the motion model and the robustness of the tracking—the finite element discretization of the tracking already represents some kind of “kinematical” regularization, in the sense that the mesh controls the size of the finite dimensional functional space of the displacement field. Figure 6 (cont.) shows tracking results for the cardiac-like problem, for various mesh sizes (from 0.1, which corresponds to the characteristic size of the image texture, like in the previous examples, to $0.1/2^4 = 0.00625$, which is smaller than the image pixel), for various levels of noise (0 and 0.1) and regularization strengths (0 and 0.1), and for two different mesh refinement strategies: either we simply track the images with the different meshes independently (left in Figure 6 (cont.)a), or we track with successively refined meshes while initializing the tracking at a given mesh refinement level with the converged solution of the previous mesh refinement level (right in Figure 6 (cont.)a). First, one can see that directly tracking images with fine meshes is impossible—mechanical regularization helps a little in some cases, but not enough to obtain satisfying solutions. Thus, to obtain fine solutions, one needs to perform multi-resolution. However, one can also see that multi-resolution is not enough, and mechanical regularization is necessary to obtain fine solutions.
<figure><img src="image_12.png" /><figcaption>**Figure 4 (cont.).** (c) Tracking solutions for the elastic regularization term (7), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The elastic regularization terms tend to prevent the mesh from deforming, which degrades the tracking.</figcaption></figure>

**4. Conclusion**

In this paper, we introduced a consistent formulation, in the nonlinear large deformation setting, of the discrete equilibrium gap principle, and used it as a regularization for the large motion tracking problem. This principle enforces that the tracked motion corresponds to the motion of a body in equilibrium with some tractions applied on its boundary. We also introduced a novel regularization of the boundary tractions involved in the equilibrium gap principle, which naturally applies to both the finite and infinitesimal strain settings.

In summary, we validated our approach by generating images representing various motions and with different signal-to-noise ratios, and then running our motion tracking algorithm with various regularization terms and regularization strengths. Basically we concluded that the finite strain framework is necessary when large motion—especially large rotation—is involved. We showed that elastic regularization interferes pathologically with nonrigid motion tracking, and should essentially not be used in such cases. In the case of heterogeneous motion, which is the most representative of actual applications of motion tracking [46, 47], we also concluded that our previous “continuous” formulation of the equilibrium gap principle interferes with the tracking, while our new “discrete” formulation performs well. Finally, we also showed that, combining
<figure><img src="image_13.png" /><figcaption>**Figure 4 (cont.).** (c) Tracking solutions for the elastic regularization term (7), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The elastic regularization terms tend to prevent the mesh from deforming, which degrades the tracking.</figcaption></figure>

## 4. Conclusion

In this paper, we introduced a consistent formulation, in the nonlinear large deformation setting, of the discrete equilibrium gap principle, and used it as a regularization for the large motion tracking problem. This principle enforces that the tracked motion corresponds to the motion of a body in equilibrium with some tractions applied on its boundary. We also introduced a novel regularization of the boundary tractions involved in the equilibrium gap principle, which naturally applies to both the finite and infinitesimal strain settings.

In summary, we validated our approach by generating images representing various motions and with different signal-to-noise ratios, and then running our motion tracking algorithm with various regularization terms and regularization strengths. Basically we concluded that the finite strain framework is necessary when large motion—especially large rotation—is involved. We showed that elastic regularization interferes pathologically with nonrigid motion tracking, and should essentially not be used in such cases. In the case of heterogeneous motion, which is the most representative of actual applications of motion tracking [46, 47], we also concluded that our previous “continuous” formulation of the equilibrium gap principle interferes with the tracking, while our new “discrete” formulation performs well. Finally, we also showed that, combining
<figure><img src="image_13.png" /><figcaption>**Figure 5.** Cardiac-like case. (a) Normalized displacement error (40) as a function of regularization strength (parameter $\beta$), for various levels of image noise, for various regularization terms (first row: $J_b^{\text{reg,el}}$ (7); second row: $J_b^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12), $J_h^{\text{reg,eq}}$ (18) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The elastic regularization terms interfere with the tracking because they basically prevent the mesh from deforming. (Continued on next page.)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 5.** Cardiac-like case. (a) Normalized displacement error (40) as a function of regularization strength (parameter $\beta$), for various levels of image noise, for various regularization terms (first row: $J_b^{\text{reg,el}}$ (7); second row: $J_b^{\text{reg,eq,cont}}$ (9); third row: $J_b^{\text{reg,eq}}$ (12); fourth row: $J_b^{\text{reg,eq}}$ (12), $J_h^{\text{reg,eq}}$ (18) and $J_{t,2D}^{\text{reg,eq}}$ (22)), and for various constitutive laws (left column: small strain, Hooke law; right column: large strain, neo-Hookean and Ogden-Ciarlet-Geymonat law). The elastic regularization terms interfere with the tracking because they basically prevent the mesh from deforming. (Continued on next page.)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 5 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization terms (18) and (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

multi-level resolution and our mechanical regularization, we were able to track rather large motion with a displacement discretization as fine as the image discretization itself, so that all features of the images motion can be tracked.

In order to promote open and reproducible science, the general motion tracking algorithm as well as all regularization terms have been implemented in an open-source finite element motion tracking library [39] (currently at https://gitlab.inria.fr/mgenet/dolfin_warp), and the specific code to reproduce the results of this paper is also freely available [43] (in the supplementary file and currently at https://mgenet.gitlabpages.inria.fr/N-DEG-paper-demos/index.html).

There are multiple limitations to our current approach, at various levels. The formulation of the equilibrium gap principle introduced here only considers a single homogeneous domain, which might be limiting when tracking motion of highly heterogeneous structures such as diseases organs [8, 48]. In such case, one might need to extend the formulation to multiple zones. Another option would be to estimate the material parameters of the mechanical model used for the regularization term at the same time as the tracking [49]. Actually, the equilibrium gap principle could also be used as a cost function for direct (i.e., without the need for a resolution of the direct mechanics problem) material parameter estimation from displacement measured by any method—the small stain formulation has already been used successfully [18], but in large
<figure><img src="image_15.png" /><figcaption>**Figure 5 (cont.)**. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization terms (18) and (22), with regularization strength $\beta = 0.1$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). Tracking is visually satisfying for all noise levels. (Continued on next page.)</figcaption></figure>

multi-level resolution and our mechanical regularization, we were able to track rather large motion with a displacement discretization as fine as the image discretization itself, so that all features of the images motion can be tracked.

In order to promote open and reproducible science, the general motion tracking algorithm as well as all regularization terms have been implemented in an open-source finite element motion tracking library [39] (currently at https://gitlab.inria.fr/mgenet/dolfin_warp), and the specific code to reproduce the results of this paper is also freely available [43] (in the supplementary file and currently at https://mgenet.gitlabpages.inria.fr/N-DEG-paper-demos/index.html).

There are multiple limitations to our current approach, at various levels. The formulation of the equilibrium gap principle introduced here only considers a single homogeneous domain, which might be limiting when tracking motion of highly heterogeneous structures such as diseases organs [8, 48]. In such case, one might need to extend the formulation to multiple zones. Another option would be to estimate the material parameters of the mechanical model used for the regularization term at the same time as the tracking [49]. Actually, the equilibrium gap principle could also be used as a cost function for direct (i.e., without the need for a resolution of the direct mechanics problem) material parameter estimation from displacement measured by any method—the small stain formulation has already been used successfully [18], but in large
<figure><img src="image_15.png" /><figcaption>**Figure 5 (cont.)**. (c) Tracking solutions for the continuous version of the nonlinear equilibrium gap regularization term (9), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The continuous version of the nonlinear equilibrium gap regularization term tends to prevent the mesh from deforming non homogeneously (which generates a discretization-induced equilibrium gap), which degrades the tracking.</figcaption></figure>

strain only the virtual fields method has been used, which represents an approximation of the equilibrium gap method as it only satisfies the equilibrium in a weak sense for a selected member of test functions [50].

Another limitation, especially in terms of applications, is that the performance of the method comes at a significant computational cost, limiting its use notably in the clinics. Thus, our method could be used to generate ground truth and/or validation data for machine learning algorithms, which could then perform tracking almost instantaneously. This is an active field of research today [51–54].

An important perspective of our motion tracking algorithm is toward low resolution images, which have been proven to drastically impact the tracking quality [55]. One option is to combine images from multiple modalities and combine the strength of both images [44]. Another option to consider is to introduce a model of the imaging modality (e.g., MRI [56]), in order to control the bias induced by the image discretization. Thus, by combining models of the imaging process and the mechanical deformation, we will be able to perform high quality tracking even with low quality images.
<figure><img src="image_16.png" /><figcaption>**Figure 5 (cont.).** (c) Tracking solutions for the continuous version of the nonlinear equilibrium gap regularization term (9), with regularization strength $\beta = 0.8$, for various noise levels: 0 (red), 0.1 (blue), 0.2 (green), 0.3 (purple). The continuous version of the nonlinear equilibrium gap regularization term tends to prevent the mesh from deforming non homogeneously (which generates a discretization-induced equilibrium gap), which degrades the tracking.</figcaption></figure>

strain only the virtual fields method has been used, which represents an approximation of the equilibrium gap method as it only satisfies the equilibrium in a weak sense for a selected member of test functions [50].

Another limitation, especially in terms of applications, is that the performance of the method comes at a significant computational cost, limiting its use notably in the clinics. Thus, our method could be used to generate ground truth and/or validation data for machine learning algorithms, which could then perform tracking almost instantaneously. This is an active field of research today [51–54].

An important perspective of our motion tracking algorithm is toward low resolution images, which have been proven to drastically impact the tracking quality [55]. One option is to combine images from multiple modalities and combine the strength of both images [44]. Another option to consider is to introduce a model of the imaging modality (e.g., MRI [56]), in order to control the bias induced by the image discretization. Thus, by combining models of the imaging process and the mechanical deformation, we will be able to perform high quality tracking even with low quality images.
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 6.** Caption continued on next page.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 6.** Caption continued on next page.</figcaption></figure>
**Figure 6 (cont.)**. Cardiac-like case, impact of mesh size. (a) Normalized displacement error (40) as a function of mesh size, for various levels of image noise, for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization terms (18) and (22), with various regularization strength. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization terms (18) and (22), with regularization strength $\beta = 0.1$.

**Declaration of interests**

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

**Acknowledgements**

MG would like to thank Philippe Moireau, Patrick Le Tallec and Pierre Kerfriden for fruitful discussions.

**Supplementary data**

Supporting information for this article is available on the journal's website under https://doi.org/10.5802/crmeca.228 or from the author.

**References**

[1] T. C. Chu, W. F. Ranson, M. A. Sutton, "Applications of digital-image-correlation techniques to experimental mechanics", *Exp. Mech.* **25** (1985), no. 3, p. 232-244.

[2] F. Hild, S. Roux, "Digital image correlation: from displacement measurement to identification of elastic properties – a review", *Strain* **42** (2006), no. 2, p. 69-80.

[3] N. Lenoir, M. Bornert, J. Desrues, P. Bésuelle, G. Viggiani, "Volumetric digital image correlation applied to X-ray microtomography images from triaxial compression tests on argillaceous rock", *Strain* **43** (2007), no. 3, p. 193-205.

[4] N. Tueni, J. Vizet, M. Genet, A. Pierangelo, J.-M. Allain, "Microstructural deformation observed by Mueller polarimetry during traction assay on myocardium samples", *Sci. Rep.* **10** (2020), no. 1, article no. 20531.

[5] J. Garot, D. A. Bluemke, N. F. Osman, C. E. Rochitte, E. R. McVeigh, E. A. Zerhouni, J. L. Prince, J. A. C. Lima, "Fast determination of regional myocardial strain fields from tagged cardiac images using harmonic phase MRI", *Circulation* **101** (2000), no. 9, p. 981-988.

[6] H. Zou, C. Xi, X. Zhao, A. S. Koh, F. Gao, Y. Su, R.-S. Tan, J. Allen, L. C. Lee, M. Genet, L. Zhong, "Quantification of biventricular strains in heart failure with preserved ejection fraction patient using hyperelastic warping method", *Front. Physiol.* **9** (2018), article no. 1295.

[7] N. P. Smith, A. Vecchi, M. McCormick, D. A. Nordsletten, O. Camara, A. F. Frangi, H. Delingette, M. Sermesant, J. Relan, N. Ayache, M. W. Krueger, W. H. W. Schulze, R. Hose, I. Valverde, P. Beerbaum, C. Staicu, M. Siebes, J. Spaan, P. J. Hunter, J. Weese, H. Lehmann, D. Chapelle, R. Rezavi, "euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling", *Interface Focus* **1** (2011), no. 3, p. 349-364.

[8] C. Patte, P.-Y. Brillet, C. Fetita, T. Gille, J.-F. Bernaudin, H. Nunes, D. Chapelle, M. Genet, "Estimation of regional pulmonary compliance in idiopathic pulmonary fibrosis based on personalized lung poromechanical modeling", *J. Biomech. Eng.* **144** (2022), no. 9, article no. 091008.

[9] M. Bornert, F. Brémand, P. Doumalin, J.-C. Dupré, M. Fazzini, M. Grédiac, F. Hild, S. Mistou, J. Molimard, J.-J. Orteu, L. Robert, Y. Surrel, P. Vacher, B. Wattrisse, "Assessment of digital image correlation measurement errors: methodology and results", *Exp. Mech.* **49** (2009), no. 3, p. 353-370.

[10] F. Hild, S. Roux, "Comparison of local and global approaches to digital image correlation", *Exp. Mech.* **52** (2012), no. 9, p. 1503-1519.

[11] A. Sotiras, C. Davatzikos, N. Paragios, "Deformable medical image registration: a survey", *IEEE Trans. Med. Imaging* **32** (2013), no. 7, p. 1153-1190.
**Figure 6 (cont.)**. Cardiac-like case, impact of mesh size. (a) Normalized displacement error (40) as a function of mesh size, for various levels of image noise, for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization terms (18) and (22), with various regularization strength. (b) Tracking solutions for the discrete version of the nonlinear equilibrium gap (12) including surface traction regularization terms (18) and (22), with regularization strength $\beta = 0.1$.

**Declaration of interests**

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

**Acknowledgements**

MG would like to thank Philippe Moireau, Patrick Le Tallec and Pierre Kerfriden for fruitful discussions.

**Supplementary data**

Supporting information for this article is available on the journal's website under https://doi.org/10.5802/crmeca.228 or from the author.

**References**

[1] T. C. Chu, W. F. Ranson, M. A. Sutton, "Applications of digital-image-correlation techniques to experimental mechanics", *Exp. Mech.* **25** (1985), no. 3, p. 232-244.

[2] F. Hild, S. Roux, "Digital image correlation: from displacement measurement to identification of elastic properties – a review", *Strain* **42** (2006), no. 2, p. 69-80.

[3] N. Lenoir, M. Bornert, J. Desrues, P. Bésuelle, G. Viggiani, "Volumetric digital image correlation applied to X-ray microtomography images from triaxial compression tests on argillaceous rock", *Strain* **43** (2007), no. 3, p. 193-205.

[4] N. Tueni, J. Vizet, M. Genet, A. Pierangelo, J.-M. Allain, "Microstructural deformation observed by Mueller polarimetry during traction assay on myocardium samples", *Sci. Rep.* **10** (2020), no. 1, article no. 20531.

[5] J. Garot, D. A. Bluemke, N. F. Osman, C. E. Rochitte, E. R. McVeigh, E. A. Zerhouni, J. L. Prince, J. A. C. Lima, "Fast determination of regional myocardial strain fields from tagged cardiac images using harmonic phase MRI", *Circulation* **101** (2000), no. 9, p. 981-988.

[6] H. Zou, C. Xi, X. Zhao, A. S. Koh, F. Gao, Y. Su, R.-S. Tan, J. Allen, L. C. Lee, M. Genet, L. Zhong, "Quantification of biventricular strains in heart failure with preserved ejection fraction patient using hyperelastic warping method", *Front. Physiol.* **9** (2018), article no. 1295.

[7] N. P. Smith, A. Vecchi, M. McCormick, D. A. Nordsletten, O. Camara, A. F. Frangi, H. Delingette, M. Sermesant, J. Relan, N. Ayache, M. W. Krueger, W. H. W. Schulze, R. Hose, I. Valverde, P. Beerbaum, C. Staicu, M. Siebes, J. Spaan, P. J. Hunter, J. Weese, H. Lehmann, D. Chapelle, R. Rezavi, "euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling", *Interface Focus* **1** (2011), no. 3, p. 349-364.

[8] C. Patte, P.-Y. Brillet, C. Fetita, T. Gille, J.-F. Bernaudin, H. Nunes, D. Chapelle, M. Genet, "Estimation of regional pulmonary compliance in idiopathic pulmonary fibrosis based on personalized lung poromechanical modeling", *J. Biomech. Eng.* **144** (2022), no. 9, article no. 091008.

[9] M. Bornert, F. Brémand, P. Doumalin, J.-C. Dupré, M. Fazzini, M. Grédiac, F. Hild, S. Mistou, J. Molimard, J.-J. Orteu, L. Robert, Y. Surrel, P. Vacher, B. Wattrisse, "Assessment of digital image correlation measurement errors: methodology and results", *Exp. Mech.* **49** (2009), no. 3, p. 353-370.

[10] F. Hild, S. Roux, "Comparison of local and global approaches to digital image correlation", *Exp. Mech.* **52** (2012), no. 9, p. 1503-1519.

[11] A. Sotiras, C. Davatzikos, N. Paragios, "Deformable medical image registration: a survey", *IEEE Trans. Med. Imaging* **32** (2013), no. 7, p. 1153-1190.
[12] C. Tobon-Gomez, M. De Craene, K. McLeod, L. Tautz, W. Shi, A. Hennemuth, A. Prakosa, H. Wang, G. S. Carr-White, S. Kapetanakis, A. Lutz, V. Rasche, T. Schaeffter, C. Butakoff, O. Friman, T. Mansi, M. Sermesant, X. Zhuang, S. Ourselin, H.-O. Peitgen, X. Pennec, R. Razavi, D. Rueckert, A. F. Frangi, K. S. Rhode, "Benchmarking framework for myocardial tracking and deformation algorithms: an open access database", *Med. Image Anal.* **17** (2013), no. 6, p. 632-648.

[13] J.-C. Passieux, J.-N. Périé, "High resolution digital image correlation using proper generalized decomposition: PGD-DIC", *Int. J. Numer. Methods Eng.* **92** (2012), no. 6, p. 531-550.

[14] G. E. Christensen, R. D. Rabbitt, M. I. Miller, "Deformable templates using large deformation kinematics", *IEEE Trans. Image Process.: A Publ. IEEE Signal Process. Soc.* **5** (1996), no. 10, p. 1435-1447.

[15] T. Mansi, X. Pennec, M. Sermesant, H. Delingette, N. Ayache, "iLogDemons: a demons-based registration algorithm for tracking incompressible elastic biological tissues", *Int. J. Comput. Vis.* **92** (2011), no. 1, p. 92-111.

[16] A. I. Veress, G. T. Gullberg, J. A. Weiss, "Measurement of strain in the left ventricle during diastole with cine-MRI and deformable image registration", *J. Biomech. Eng.* **127** (2005), no. 7, p. 1195-1207.

[17] M. Genet, C. T. Stoeck, C. Deuster, L. C. Lee, S. Kozerke, "Equilibrated warping: finite element image registration with finite strain equilibrium gap regularization", *Med. Image Anal.* **50** (2018), p. 1-22.

[18] D. Claire, F. Hild, S. Roux, "A finite element formulation to identify damage fields: the equilibrium gap method", *Int. J. Numer. Methods Eng.* **61** (2004), no. 2, p. 189-208.

[19] H. Leclerc, J.-N. Périé, S. Roux, F. Hild, "Voxel-scale digital volume correlation", *Exp. Mech.* **51** (2010), no. 4, p. 479-490.

[20] L. C. Lee, M. Genet, "Validation of equilibrated warping—image registration with mechanical regularization—on 3D ultrasound images", in *Functional Imaging and Modeling of the Heart (FIMH)* (Y. Coudière, V. Ozenne, E. Vigmond, N. Zemzemi, eds.), vol. 11504, Springer International Publishing, Cham, 2019, p. 334-341.

[21] E. Berberoğlu, C. T. Stoeck, P. Moireau, S. Kozerke, M. Genet, "Validation of finite element image registration-based cardiac strain estimation from magnetic resonance images", *Proc. Appl. Math. Mech.* **19** (2019), no. 1, article no. e201900418.

[22] C. Xi, C. Latnie, X. Zhao, J. L. Tan, S. T. Wall, M. Genet, L. Zhong, L. C. Lee, "Patient-specific computational analysis of ventricular mechanics in pulmonary arterial hypertension", *J. Biomech. Eng.* **138** (2016), no. 11, article no. 111001.

[23] H. Zou, S. Leng, C. Xi, X. Zhao, A. S. Koh, F. Gao, J. L. Tan, R.-S. Tan, J. C. Allen, L. C. Lee, M. Genet, L. Zhong, "Three-dimensional biventricular strains in pulmonary arterial hypertension patients using hyperelastic warping", *Comput. Methods Programs Biomed.* **189** (2020), article no. 105345.

[24] D. A. Castellanos, K. Škardová, A. Bhattaru, E. Berberoğlu, G. Greil, A. Tandon, J. Dillenbeck, B. Burkhardt, T. Hussain, M. Genet, R. Chabiniok, "Left ventricular torsion obtained using equilibrated warping in patients with repaired tetralogy of fallot", *Pediatr. Cardiol.* **42** (2021), p. 1275-1283.

[25] M. I. Miller, G. E. Christensen, Y. Amit, U. Grenander, "Mathematical textbook of deformable neuroanatomies", *Proc. Natl. Acad. Sci. USA* **90** (1993), no. 24, p. 11944-11948.

[26] J. Réthoré, S. Roux, F. Hild, "An extended and integrated digital image correlation technique applied to the analysis of fractured samples: The equilibrium gap method as a mechanical filter", *Eur. J. Comput. Mech.* **18** (2009), no. 3-4, p. 285-306.

[27] P. Ladevèze, J. P. Pelle, *Mastering Calculations in Linear and Nonlinear Mechanics*, Mechanical Engineering Series, Springer Science, New York, 2005.

[28] O. Zienkiewicz, R. Taylor, J. Zhu, "Ch. 15—Errors, recovery processes, and error estimates", in *The Finite Element Method: its Basis and Fundamentals*, Elsevier, Amsterdam, 2013, p. 493-543.

[29] A. Mendoza, J. Neggers, F. Hild, S. Roux, "Complete mechanical regularization applied to digital image and volume correlation", *Comput. Methods Appl. Mech. Eng.* **355** (2019), p. 27-43.

[30] P. Brandner, T. Jankuhn, S. Praetorius, A. Reuksen, A. Voigt, "Finite element discretization methods for velocity-pressure and stream function formulations of surface Stokes Equations", *SIAM J. Sci. Comput.* **44** (2022), no. 4, p. A1807-A1832.

[31] G. Allaire, *Conception Optimale Des Structures*, Mathématiques & Applications, vol. 58, Springer, Berlin, 2007.

[32] J. Passieux, R. Bouclier, "Classic and inverse compositional Gauss-Newton in global DIC", *Int. J. Numer. Methods Eng.* **119** (2019), no. 6, p. 453-468.

[33] R. W. Ogden, "Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids", *Proc. R. Soc. Lond. A. Math. Phys. Sci.* **328** (1972), no. 1575, p. 567-583.

[34] P. G. Ciarlet, G. Geymonat, "Sur Les Lois de Comportement En Élasticité Non-Linéaire Compressible", *C. R. Acad. Sci. Sér. II* **295** (1982), p. 423-426.

[35] P. Le Tallec, "Numerical methods for nonlinear elasticity", in *Handbook of Numerical Analysis*, vol. 3, Elsevier, 1994, p. 465-622.

[36] M. Genet, "A relaxed growth modeling framework for controlling growth-induced residual stresses", *Clin. Biomech.* **70** (2019), p. 270-277.
[12] C. Tobon-Gomez, M. De Craene, K. McLeod, L. Tautz, W. Shi, A. Hennemuth, A. Prakosa, H. Wang, G. S. Carr-White, S. Kapetanakis, A. Lutz, V. Rasche, T. Schaeffter, C. Butakoff, O. Friman, T. Mansi, M. Sermesant, X. Zhuang, S. Ourselin, H.-O. Peitgen, X. Pennec, R. Razavi, D. Rueckert, A. F. Frangi, K. S. Rhode, "Benchmarking framework for myocardial tracking and deformation algorithms: an open access database", *Med. Image Anal.* **17** (2013), no. 6, p. 632-648.

[13] J.-C. Passieux, J.-N. Périé, "High resolution digital image correlation using proper generalized decomposition: PGD-DIC", *Int. J. Numer. Methods Eng.* **92** (2012), no. 6, p. 531-550.

[14] G. E. Christensen, R. D. Rabbitt, M. I. Miller, "Deformable templates using large deformation kinematics", *IEEE Trans. Image Process.: A Publ. IEEE Signal Process. Soc.* **5** (1996), no. 10, p. 1435-1447.

[15] T. Mansi, X. Pennec, M. Sermesant, H. Delingette, N. Ayache, "iLogDemons: a demons-based registration algorithm for tracking incompressible elastic biological tissues", *Int. J. Comput. Vis.* **92** (2011), no. 1, p. 92-111.

[16] A. I. Veress, G. T. Gullberg, J. A. Weiss, "Measurement of strain in the left ventricle during diastole with cine-MRI and deformable image registration", *J. Biomech. Eng.* **127** (2005), no. 7, p. 1195-1207.

[17] M. Genet, C. T. Stoeck, C. Deuster, L. C. Lee, S. Kozerke, "Equilibrated warping: finite element image registration with finite strain equilibrium gap regularization", *Med. Image Anal.* **50** (2018), p. 1-22.

[18] D. Claire, F. Hild, S. Roux, "A finite element formulation to identify damage fields: the equilibrium gap method", *Int. J. Numer. Methods Eng.* **61** (2004), no. 2, p. 189-208.

[19] H. Leclerc, J.-N. Périé, S. Roux, F. Hild, "Voxel-scale digital volume correlation", *Exp. Mech.* **51** (2010), no. 4, p. 479-490.

[20] L. C. Lee, M. Genet, "Validation of equilibrated warping—image registration with mechanical regularization—on 3D ultrasound images", in *Functional Imaging and Modeling of the Heart (FIMH)* (Y. Coudière, V. Ozenne, E. Vigmond, N. Zemzemi, eds.), vol. 11504, Springer International Publishing, Cham, 2019, p. 334-341.

[21] E. Berberoğlu, C. T. Stoeck, P. Moireau, S. Kozerke, M. Genet, "Validation of finite element image registration-based cardiac strain estimation from magnetic resonance images", *Proc. Appl. Math. Mech.* **19** (2019), no. 1, article no. e201900418.

[22] C. Xi, C. Latnie, X. Zhao, J. L. Tan, S. T. Wall, M. Genet, L. Zhong, L. C. Lee, "Patient-specific computational analysis of ventricular mechanics in pulmonary arterial hypertension", *J. Biomech. Eng.* **138** (2016), no. 11, article no. 111001.

[23] H. Zou, S. Leng, C. Xi, X. Zhao, A. S. Koh, F. Gao, J. L. Tan, R.-S. Tan, J. C. Allen, L. C. Lee, M. Genet, L. Zhong, "Three-dimensional biventricular strains in pulmonary arterial hypertension patients using hyperelastic warping", *Comput. Methods Programs Biomed.* **189** (2020), article no. 105345.

[24] D. A. Castellanos, K. Škardová, A. Bhattaru, E. Berberoğlu, G. Greil, A. Tandon, J. Dillenbeck, B. Burkhardt, T. Hussain, M. Genet, R. Chabiniok, "Left ventricular torsion obtained using equilibrated warping in patients with repaired tetralogy of fallot", *Pediatr. Cardiol.* **42** (2021), p. 1275-1283.

[25] M. I. Miller, G. E. Christensen, Y. Amit, U. Grenander, "Mathematical textbook of deformable neuroanatomies", *Proc. Natl. Acad. Sci. USA* **90** (1993), no. 24, p. 11944-11948.

[26] J. Réthoré, S. Roux, F. Hild, "An extended and integrated digital image correlation technique applied to the analysis of fractured samples: The equilibrium gap method as a mechanical filter", *Eur. J. Comput. Mech.* **18** (2009), no. 3-4, p. 285-306.

[27] P. Ladevèze, J. P. Pelle, *Mastering Calculations in Linear and Nonlinear Mechanics*, Mechanical Engineering Series, Springer Science, New York, 2005.

[28] O. Zienkiewicz, R. Taylor, J. Zhu, "Ch. 15—Errors, recovery processes, and error estimates", in *The Finite Element Method: its Basis and Fundamentals*, Elsevier, Amsterdam, 2013, p. 493-543.

[29] A. Mendoza, J. Neggers, F. Hild, S. Roux, "Complete mechanical regularization applied to digital image and volume correlation", *Comput. Methods Appl. Mech. Eng.* **355** (2019), p. 27-43.

[30] P. Brandner, T. Jankuhn, S. Praetorius, A. Reuksen, A. Voigt, "Finite element discretization methods for velocity-pressure and stream function formulations of surface Stokes Equations", *SIAM J. Sci. Comput.* **44** (2022), no. 4, p. A1807-A1832.

[31] G. Allaire, *Conception Optimale Des Structures*, Mathématiques & Applications, vol. 58, Springer, Berlin, 2007.

[32] J. Passieux, R. Bouclier, "Classic and inverse compositional Gauss-Newton in global DIC", *Int. J. Numer. Methods Eng.* **119** (2019), no. 6, p. 453-468.

[33] R. W. Ogden, "Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids", *Proc. R. Soc. Lond. A. Math. Phys. Sci.* **328** (1972), no. 1575, p. 567-583.

[34] P. G. Ciarlet, G. Geymonat, "Sur Les Lois de Comportement En Élasticité Non-Linéaire Compressible", *C. R. Acad. Sci. Sér. II* **295** (1982), p. 423-426.

[35] P. Le Tallec, "Numerical methods for nonlinear elasticity", in *Handbook of Numerical Analysis*, vol. 3, Elsevier, 1994, p. 465-622.

[36] M. Genet, "A relaxed growth modeling framework for controlling growth-induced residual stresses", *Clin. Biomech.* **70** (2019), p. 270-277.
[37] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, *Numerical Recipes: The Art of Scientific Computing*, 3rd ed., Cambridge University Press, Cambridge, NY, 2007.

[38] R. Fedele, L. Galantucci, A. Ciani, "Global 2D digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation", *Int. J. Numer. Methods Eng.* **96** (2013), no. 12, p. 739-762.

[39] M. Genet, "Dolfin_warp", 2023, Zenodo.

[40] A. Logg, K.-A. Mardal, G. Wells (eds.), *Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book*, Lecture Notes in Computational Science and Engineering, vol. 84, Springer, Heidelberg, 2012.

[41] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells, "The FEniCS project version 1.5", in *Arch. Numer. Softw.*, vol. 3, 2015, p. 9-23.

[42] W. Schroeder, K. Martin, B. Lorensen, *The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics*, 4th ed., Kitware, Inc., Clifton Park, NY, 2006.

[43] M. Genet, "N-DEG-paper-demos", 2023, Zenodo.

[44] E. Berberoğlu, C. T. Stoeck, S. Kozerke, M. Genet, "Quantification of left ventricular strain and torsion by joint analysis of 3D tagging and cine MR images", *Med. Image Anal.* **82** (2022), article no. 102598.

[45] A. K. Rutz, S. Ryf, S. Plein, P. Boesiger, S. Kozerke, "Accelerated whole-heart 3D CSPAMM for myocardial motion quantification", *Magn. Reson. Med.* **59** (2008), no. 4, p. 755-763.

[46] M. Genet, L. C. Lee, R. Nguyen, H. Haraldsson, G. Acevedo-Bolton, Z. Zhang, L. Ge, K. Ordovas, S. Kozerke, J. M. Guccione, "Distribution of normal human left ventricular myofiber stress at end diastole and end systole: a target for in silico design of heart failure treatments", *J. Appl. Phys.* **117** (2014), p. 142-152.

[47] H. Finsberg, C. Xi, X. Zhao, J. L. Tan, M. Genet, J. Sundnes, L. C. Lee, L. Zhong, S. T. Wall, "Computational quantification of patient-specific changes in ventricular dynamics associated with pulmonary hypertension", *Am. J. Physiol. Heart Circ. Physiol.* **317** (2019), no. 6, p. H1363-H1375.

[48] C. Laville, C. Fetita, T. Gille, P.-Y. Brillet, H. Nunes, J.-F. Bernaudin, M. Genet, "Comparison of optimization parametrizations for regional lung compliance estimation using personalized pulmonary poromechanical modeling", *Biomech. Model. Mechanobiol.* **22** (2023), p. 1541-1554.

[49] F. Mathieu, H. Leclerc, F. Hild, S. Roux, "Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC", *Exp. Mech.* **55** (2015), no. 1, p. 105-119.

[50] S. Avril, M. Bonnet, A.-S. Bretelle, M. Grédiac, F. Hild, P. Lenny, F. Latourte, D. Lemosse, S. Pagano, E. Pagnacco, E. Pierron, "Overview of identification methods of mechanical parameters based on full-field measurements", *Exp. Mech.* **48** (2008), no. 4, p. 381-402.

[51] T. Leiner, D. Rueckert, A. Suinesiaputra, B. Baeßler, R. Nezafat, I. Iğum, A. A. Young, "Machine learning in cardiovascular magnetic resonance: basic concepts and applications", *J. Cardiovasc. Magn. Resonan.* **21** (2019), no. 1, p. 1-14.

[52] S. Friedrich, S. Groß, I. R. König, S. Engelhardt, M. Bahls, J. Heinz, C. Huber, L. Kaderali, M. Kelm, A. Leha, J. Rühl, J. Schaller, C. Scherer, M. Vollmer, T. Seidler, T. Friede, "Applications of artificial intelligence/machine learning approaches in cardiovascular medicine: a systematic review with recommendations", *Eur. Heart J. - Digital Health* **2** (2021), no. 3, p. 424-436.

[53] S. Koehler, T. Hussain, H. Hussain, D. Young, S. Sarikouch, T. Pickardt, G. Greil, S. Engelhardt, "Self-supervised motion descriptor for cardiac phase detection in 4D CMR based on discrete vector field estimations", in *Statistical Atlases and Computational Models of the Heart. Regular and CMRxMotion Challenge Papers* (O. Camara, E. Puyol-Antón, C. Qin, M. Sermesant, A. Suinesiaputra, S. Wang, A. Young, eds.), vol. 13593, Springer Nature, Switzerland, 2022, p. 65-78.

[54] P. Arratia López, H. Mella, S. Uribe, D. E. Hurtado, F. Sahli Costabal, "WarpPINN: Cine-MR image registration with physics-informed neural networks", *Med. Image Anal.* **89** (2023), article no. 102925.

[55] E. Berberoğlu, C. T. Stoeck, P. Moireau, S. Kozerke, M. Genet, "In-silico study of accuracy and precision of left-ventricular strain quantification from 3D tagged MRI", *PLoS One* **16** (2021), no. 11, article no. e0258965.

[56] K. Škardová, M. Rambausek, R. Chabiniok, M. Genet, "Mechanical and imaging models-based image registration", in *VipIMAGE 2019 (J. M. R. S. Tavares, R. M. Natal Jorge, eds.)*, vol. 34, Springer International Publishing, Cham, 2019, p. 77-85.
[37] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, *Numerical Recipes: The Art of Scientific Computing*, 3rd ed., Cambridge University Press, Cambridge, NY, 2007.

[38] R. Fedele, L. Galantucci, A. Ciani, "Global 2D digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation", *Int. J. Numer. Methods Eng.* **96** (2013), no. 12, p. 739-762.

[39] M. Genet, "Dolfin_warp", 2023, Zenodo.

[40] A. Logg, K.-A. Mardal, G. Wells (eds.), *Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book*, Lecture Notes in Computational Science and Engineering, vol. 84, Springer, Heidelberg, 2012.

[41] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells, "The FEniCS project version 1.5", in *Arch. Numer. Softw.*, vol. 3, 2015, p. 9-23.

[42] W. Schroeder, K. Martin, B. Lorensen, *The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics*, 4th ed., Kitware, Inc., Clifton Park, NY, 2006.

[43] M. Genet, "N-DEG-paper-demos", 2023, Zenodo.

[44] E. Berberoğlu, C. T. Stoeck, S. Kozerke, M. Genet, "Quantification of left ventricular strain and torsion by joint analysis of 3D tagging and cine MR images", *Med. Image Anal.* **82** (2022), article no. 102598.

[45] A. K. Rutz, S. Ryf, S. Plein, P. Boesiger, S. Kozerke, "Accelerated whole-heart 3D CSPAMM for myocardial motion quantification", *Magn. Reson. Med.* **59** (2008), no. 4, p. 755-763.

[46] M. Genet, L. C. Lee, R. Nguyen, H. Haraldsson, G. Acevedo-Bolton, Z. Zhang, L. Ge, K. Ordovas, S. Kozerke, J. M. Guccione, "Distribution of normal human left ventricular myofiber stress at end diastole and end systole: a target for in silico design of heart failure treatments", *J. Appl. Phys.* **117** (2014), p. 142-152.

[47] H. Finsberg, C. Xi, X. Zhao, J. L. Tan, M. Genet, J. Sundnes, L. C. Lee, L. Zhong, S. T. Wall, "Computational quantification of patient-specific changes in ventricular dynamics associated with pulmonary hypertension", *Am. J. Physiol. Heart Circ. Physiol.* **317** (2019), no. 6, p. H1363-H1375.

[48] C. Laville, C. Fetita, T. Gille, P.-Y. Brillet, H. Nunes, J.-F. Bernaudin, M. Genet, "Comparison of optimization parametrizations for regional lung compliance estimation using personalized pulmonary poromechanical modeling", *Biomech. Model. Mechanobiol.* **22** (2023), p. 1541-1554.

[49] F. Mathieu, H. Leclerc, F. Hild, S. Roux, "Estimation of elastoplastic parameters via weighted FEMU and integrated-DIC", *Exp. Mech.* **55** (2015), no. 1, p. 105-119.

[50] S. Avril, M. Bonnet, A.-S. Bretelle, M. Grédiac, F. Hild, P. Lenny, F. Latourte, D. Lemosse, S. Pagano, E. Pagnacco, E. Pierron, "Overview of identification methods of mechanical parameters based on full-field measurements", *Exp. Mech.* **48** (2008), no. 4, p. 381-402.

[51] T. Leiner, D. Rueckert, A. Suinesiaputra, B. Baeßler, R. Nezafat, I. Iğum, A. A. Young, "Machine learning in cardiovascular magnetic resonance: basic concepts and applications", *J. Cardiovasc. Magn. Resonan.* **21** (2019), no. 1, p. 1-14.

[52] S. Friedrich, S. Groß, I. R. König, S. Engelhardt, M. Bahls, J. Heinz, C. Huber, L. Kaderali, M. Kelm, A. Leha, J. Rühl, J. Schaller, C. Scherer, M. Vollmer, T. Seidler, T. Friede, "Applications of artificial intelligence/machine learning approaches in cardiovascular medicine: a systematic review with recommendations", *Eur. Heart J. - Digital Health* **2** (2021), no. 3, p. 424-436.

[53] S. Koehler, T. Hussain, H. Hussain, D. Young, S. Sarikouch, T. Pickardt, G. Greil, S. Engelhardt, "Self-supervised motion descriptor for cardiac phase detection in 4D CMR based on discrete vector field estimations", in *Statistical Atlases and Computational Models of the Heart. Regular and CMRxMotion Challenge Papers* (O. Camara, E. Puyol-Antón, C. Qin, M. Sermesant, A. Suinesiaputra, S. Wang, A. Young, eds.), vol. 13593, Springer Nature, Switzerland, 2022, p. 65-78.

[54] P. Arratia López, H. Mella, S. Uribe, D. E. Hurtado, F. Sahli Costabal, "WarpPINN: Cine-MR image registration with physics-informed neural networks", *Med. Image Anal.* **89** (2023), article no. 102925.

[55] E. Berberoğlu, C. T. Stoeck, P. Moireau, S. Kozerke, M. Genet, "In-silico study of accuracy and precision of left-ventricular strain quantification from 3D tagged MRI", *PLoS One* **16** (2021), no. 11, article no. e0258965.

[56] K. Škardová, M. Rambausek, R. Chabiniok, M. Genet, "Mechanical and imaging models-based image registration", in *VipIMAGE 2019 (J. M. R. S. Tavares, R. M. Natal Jorge, eds.)*, vol. 34, Springer International Publishing, Cham, 2019, p. 77-85.