\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,published]{cedram}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.935\linewidth}{#1}}
\usepackage{etoolbox}
\usepackage{upgreek}
\RequirePackage{caption}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\newcounter{runlevel}
\usepackage{multirow} 
\usepackage[boxed]{algorithm2e}
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20230521}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\0{\phantom{0}}
\def\ndash{\text{--}}
\def\rmmu{\upmu}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\def\beqimg#1#2{}
\def\eeqimg#1#2{}
\def\algOneREF{\ref}
\def\algTwoREF{\ref}
\def\algThreeREF{\ref}
\def\algFourREF{\ref}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\def\bcaption#1#2#3{\caption{\raggedright Caption continued on next page.}\end{figure}\setcounter{figure}{#2}\def\thefigure{6 (cont.)}\begin{figure}\caption{\space#1\space #3}}
%\def\bbPi{\boldsymbol{\Pi}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\input pdf-trans
\newbox\qbox
\def\usecolor#1{\csname\string\color@#1\endcsname\space}
\newcommand\mycolsplit[2]{\colorlet{tmpcolor}{#2}\edef\tmp{\usecolor{tmpcolor}}%
  \def\tmpB{}\expandafter\colsplithelp\tmp\relax%
  \ifnum0=#1\relax\edef\fillcol{\tmpB}\else\edef\bordercol{\tmpC}\fi}
\def\colsplithelp#1#2 #3\relax{%
  \edef\tmpB{\tmpB#1#2 }%
  \ifnum `#1>`9\relax\def\tmpC{#3}\else\colsplithelp#3\relax\fi
}
\newcommand\outline[1]{\leavevmode%
  \def\maltext{$#1$}%
  \setbox\qbox=\hbox{\maltext}%
  \boxgs{Q q 2 Tr \bbthickness\space w \fillcol\space \bordercol\space}{}%
  \copy\qbox%
}
\newcommand\bordercolor[1]{\mycolsplit{1}{#1}}
\newcommand\fillcolor[1]{\mycolsplit{0}{#1}}
\bordercolor{black}
\fillcolor{white}
\newcommand\bbthickness{.45}
\newcommand\mymathbb[1]{%\def\mydelim{$}
\outline{#1}} %$
\newcommand{\bbPi}{\mymathbb{\Pi}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\DOI{10.5802/crmeca.228}
\datereceived{2023-06-07}
\daterevised{2023-08-22}
\dateaccepted{2023-09-26}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage[T1]{fontenc} 
\usepackage{upgreek}
\def\CDRdoi{10.5802/crmeca.228}
\def\newline{\unskip\break}
\let\beqimg\relax
\let\eeqimg\relax
\def\bbPi{\mathbb{\Pi}}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
\def\algOneREF#1{1}
\def\algTwoREF#1{2}
\def\algThreeREF#1{3}
\def\algFourREF#1{4}
\def\bcaption#1#2#3{\caption{#1\space#3}}
}
\makeatother

\usepackage{moreverb}
\usepackage{upgreek}

\dateposted{2023-12-08}
\begin{document}

\begin{noXML}

\TopicEF{Research article}{Article de recherche}

\title{Finite strain formulation of the discrete equilibrium gap
principle: application to mechanically consistent regularization for
large motion tracking{\vspace*{-2pt}}}

\alttitle{Formulation en transformation finie du principe d'\'ecart
d'\'equilibre discret~: application \`a la r\'egularisation m\'ecanique
pour le suivi de mouvement{\vspace*{-2pt}}}

\author{\firstname{Martin} \lastname{Genet}\CDRorcid{0000-0003-2204-201X}}
\address{\mbox{Laboratoire de M\'ecanique des Solides, \'Ecole
Polytechnique/IPP/CNRS, Palaiseau,~France}}
\address{\'Equipe M$\Xi$DISIM, INRIA, Palaiseau, France}
\email{martin.genet@polytechnique.edu}

\keywords{\kwd{Motion tracking}\kwd{Mechanical
regularization}\kwd{Equilibrium gap principle}\kwd{Finite element
method}\kwd{Inverse problems}}

\altkeywords{\kwd{Suivi de mouvement}\kwd{R\'egularisation
m\'ecanique}\kwd{Principe de l'\'ecart d'\'equilibre}\kwd{M\'ethode des
\'el\'ements finis}\kwd{Probl\`emes inverses}}

\begin{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.
\end{abstract}

\begin{altabstract}
Le principe de l'\'ecart d'\'equilibre offre un bon compromis entre la
robustesse et la pr\'ecision pour la r\'egularisation du suivi du
mouvement, car il impose simplement que le mouvement suivi corresponde
\`a celui d'un corps se d\'eformant sous l'effet de charges
arbitraires. Cet article pr\'esente une extension du principe de
l'\'ecart d'\'equilibre dans le cadre des grandes d\'eformations, un
nouveau terme de r\'egularisation pour contr\^oler les tractions de
surface, les deux dans le contexte du suivi de mouvement par
\'el\'ements finis, et une reformulation {\og}~probl\`eme inverse~{\fg}
coh\'erente du probl\`eme de suivi de mouvement avec r\'egularisation
m\'ecanique. Les performances de suivi de la m\'ethode propos\'ee, avec
une r\'esolution de d\'eplacement allant jusqu'\`a la taille du pixel
de l'image, sont d\'emontr\'ees sur des images synth\'etiques
repr\'esentant divers mouvements avec diff\'erents rapports
signal-sur-bruit.
{\vspace*{-2pt}}
\end{altabstract}

\maketitle

\end{noXML}

\section{Introduction} \label{sec::intro}

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~\cite{Chu1985EM,
Hild2006S, Lenoir2007S, Tueni2020SR}. 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~\cite{Garot2000C, Zou2018FP}, or by
merging physical models and imaging data into so called digital
twins~\cite{Smith2011If, Patte2022JBE}.

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.~\cite{Bornert2009EM, Hild2012EM,
Sotiras2013ITMI, Tobon-Gomez2013MIA}. 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.~\cite{Sotiras2013ITMI, Tobon-Gomez2013MIA}.

One key question in motion tracking is the regularization, i.e., the
\textit{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~\cite{Passieux2012IJNME}, fluid-like mechanical
regularization~\cite{Christensen1996ITIPPISPS},
incompressibility~\cite{Mansi2011IJCV}, hyperelastic
energy~\cite{Veress2005JBE}, 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~\cite{Hild2006S,
Genet2018MIA}, 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~\cite{Claire2004IJNME}, and was later used in the context of motion
tracking notably in~\cite{Leclerc2010EM}. An extension to the non
linear setting was proposed, at the continuous level,
in~\cite{Genet2018MIA, Lee2019FIMHFa, Berberoglu2019PAMM}, with
multiple applications to biomedical images~\cite{Xi2016JBE,
Zou2020CMPB, Castellanos2021PC}. 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~\cite{Leclerc2010EM, Genet2018MIA}. To \mbox{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~\cite{Leclerc2010EM}. A
surface Laplacian of the displacement was used as an additional
regularization term in~\cite{Leclerc2010EM}, 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~\ref{sec::methods::problem}),
then we provide a short literature review on mechanical regularization
(Section~\ref{sec::methods::literature_review}), and we describe our
proposed regularization term, including body
(Section~\ref{sec::methods::nldeg::body}) {and} boundary
(Section~\ref{sec::methods::nldeg::boundary}) terms. Then, we give an
inverse problem formulation to the regularized motion tracking problem
(Section~\ref{sec::methods::inverse_problem}), and we describe our
numerical strategy for the resolution
(Section~\ref{sec::methods::resolution}). We finish the Methods section
with a description of the synthetic images that will be used for the
validation of our method (Section~\ref{sec::methods::data}), and then
analyse tracking Results in the case of rigid motion
(Section~\ref{sec::results::rigid}), non rigid but homogeneous motion
(Section~\ref{sec::results::nonrigid}) and non homogeneous motion
(Section~\ref{sec::results::heart}), as well as in the case of refined
meshes (Section~\ref{sec::results::refine}).

\section{Methods} \label{sec::methods}
\subsection{The motion tracking problem} \label{sec::methods::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$:
{\begin{equation} \label{eq1}
I_0:\begin{cases}
\square_0 \rightarrow \mathbb{R}\\
\underline{X} \mapsto I_0(\underline{X}),
\end{cases} \quad I:\begin{cases}
\square \rightarrow \mathbb{R}\\
\underline{x} \mapsto I(\underline{x}),
\end{cases}
\end{equation}}\noindent
where $\text{\ensuremath{\square}}_0$ {and}
$\text{\ensuremath{\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
$\underline{\Phi}$, or equivalently the smooth displacement field
$\underline{U}$, between material points of the reference and deformed
domains:
{\begin{equation} \label{eq2}
\underline{\Phi}:\begin{cases}
\Omega_0 \rightarrow \omega\\
\underline{X} \mapsto \underline{x} = \underline{\Phi}(\underline{X}),
\end{cases} \quad \underline{U}:\begin{cases}
\Omega_0 \rightarrow \mathbb{R}^3\\
\underline{X} \mapsto \underline{U}(\underline{X}) := \underline{\Phi}
(\underline{X}) - \underline{X},
\end{cases}
\end{equation}}\noindent
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:
{\begin{equation}
\text{Find}\ \underline{U^{\mathrm{sol}}}:=\text{argmin}_{\{\underline{U}\}}
\left\{J({\underline{U}}) := ({1-\beta}) \frac{J^{\mathrm{ima}}
({\underline{U}})}{J^{\mathrm{ima}}_0} + \beta \frac{J^{\mathrm{reg}}
({\underline{U}})}{J^{\mathrm{reg}}_0}\right\},
\label{eq::problem_cont}
\end{equation}}\noindent
where $J^{\mathrm{ima}}$ is the image similarity metric, or
``correlation energy'', $J^{\mathrm{reg}}$ is the regularization
energy, $J^{\mathrm{ima}}_0$ {and} $J^{\mathrm{reg}}_0$ are
normalization terms, and $\beta$ defines the regularization strength.
The normalization terms $J^{\mathrm{ima}}_0$ {and} $J^{\mathrm{reg}}_0$
allow for the consistent addition of ``energies'' with very different
physical units, and are typically taken equal to the value of
$J^{\mathrm{ima}}$ {and} $J^{\mathrm{reg}}$ for a chosen displacement
field~\cite{Leclerc2010EM}, 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:
{\begin{equation}
J^{\mathrm{ima}}({\underline{U}}) := \frac{1}{2} \int_{\Omega_0}
({I({\underline{X}+\underline{U}(\underline{X})})-I_0(\underline{X})})^2
\,\mathrm{d}\Omega_0. \label{eq::J_ima_cont}
\end{equation}}\noindent
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{\underline{\mathbb{N}}}}(\underline{X}) \cdot
\underline{\mathbb{U}}$ with $\underline{\underline{\mathbb{N}}}$ the
array of shape functions. Thus, the problem becomes a finite
dimensional problem:
{\begin{equation}
\text{Find}\ \underline{\mathbb{U}^{\mathrm{sol}}} :=
\text{argmin}_{\{\underline{\mathbb{U}}\}}\left\{J({\underline{\mathbb{U}}})
:= ({1-\beta}) \frac{J^{\mathrm{ima}}({\underline{\mathbb{U}}})}
{J^{\mathrm{ima}}_0} + \beta\frac{J^{\mathrm{reg}}
({\underline{\mathbb{U}}})}{J^{\mathrm{reg}}_0}\right\},
\label{eq::problem_disc}
\end{equation}}\noindent
where
{\begin{equation}
J^{\mathrm{ima}}({\underline{\mathbb{U}}}):=\frac{1}{2}\int_{\Omega_0}
\left({I\left({\underline{X} + {}^{t}{\underline{\underline{\mathbb{N}}}}
(\underline{X}) \cdot \underline{\mathbb{U}}}\right)- I_0(\underline{X})}
\right)^2 \,\mathrm{d}\Omega_0. \label{eq::J_ima_disc}
\end{equation}}\noindent
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.

\subsection{Short literature review on mechanical regularization}
\label{sec::methods::literature_review}

Many regularization terms have been proposed in the literature, see for
instance~\cite{Sotiras2013ITMI}. 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~\cite{Miller1993PNASU} {and} hyperelastic~\cite{Veress2005JBE}
regularizations have been proposed, which consist in penalizing the
strain energy of the body:
{\begin{equation}
J^{\mathrm{reg},\mathrm{el}}({\underline{U}}) := \int_{\Omega_0}
\rho_0 \Psi({\underline{U}})\,\mathrm{d}\Omega_0,
\label{eq::J_reg_el_cont}
\end{equation}}\noindent
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~\cite{Claire2004IJNME}. It
was formulated in the linear setting, and directly at the discrete
level, in~\cite{Rethore2009EJoCM, Leclerc2010EM}.\vspace*{2pt}
We first define $\underline{\underline{\mathbb{K}}} := \int_{\Omega_0}
\underline{\underline{\underline{\mathbb{B}}}} :
\underline{\underline{\underline{\underline{K}}}} :
{}^t{\underline{\underline{\underline{\mathbb{B}}}}}\,
\mathrm{d}\Omega_0$ the system stiffness matrix, with
$\underline{\underline{\underline{\underline{K}}}}$ the material
stiffness tensor (like for elastic regularization, a generic isotropic
Hooke law is usually used, such that
$\underline{\underline{\underline{\underline{K}}}} = \lambda
\underline{\underline{1}} \otimes \underline{\underline{1}} + 2\mu
\underline{\underline{\underline{\underline{1}}}}$ with $\lambda$ and
$\mu$ the Lam\'e constants) and
$\underline{\underline{\underline{\mathbb{B}}}}$ the array of shape
functions symmetric gradients (i.e., such that
$\underline{\underline{\varepsilon}}(\underline{X}) :=
{(\underline{\underline{\text{Grad}}}({\underline{U}})
(\underline{X}))_{\mathrm{sym}}} = ({1}/{2})
({{}^{t}{\underline{\underline{\text{Grad}}}({\underline{U}})}
(\underline{X}) + \underline{\underline{\text{Grad}}}
({\underline{U}})(\underline{X})})$ $\approx$ ${}^{t}
\underline{\underline{\underline{\mathbb{B}}}} (\underline{X}) \cdot
\underline{\mathbb{U}}$). Then the regularization is expressed as
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq},\mathrm{lin}}({\underline{\mathbb{U}}}):=
\tfrac{1}{2} {}^{t}{\underline{\mathbb{U}}} \cdot {}^{t}
{\underline{\underline{\mathbb{K}^*}}} \cdot
\underline{\underline{\mathbb{K}^*}} \cdot \underline{\mathbb{U}},
\label{eq::J_reg_eq_lin}
\end{equation}}\noindent
where $\underline{\underline{\mathbb{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 \mbox{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~\cite{Leclerc2010EM}) 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{\break} ``energy''.

A first attempt toward extending this principle to the large
deformation nonlinear setting has been made in~\cite{Genet2018MIA}.
{The main idea is to directly penalize the non verification of the
internal linear momentum equilibrium equation ($\underline{\text{Div}}
({\underline{\underline{P}}}) = 0$, where $\underline{\underline{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:}
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq},\mathrm{cont}}({\underline{U}}):=\frac{1}{2}
\sum_{{E}}\int_{{E}}{\|\underline{\text{Div}}
({\underline{\underline{P}}})\|}^2
\,\mathrm{d}E + \frac{1}{2h} \sum_F \int_F [[{\underline{\underline{P}}
\cdot \underline{N}}]]^2\, \mathrm{d}F \label{eq::J_reg_eq_cont}
\end{equation}}\noindent
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 (${}^{t}{\underline{\underline{P}}} \cdot
{}^{t}{\underline{\underline{F}}}^{-1} = \underline{\underline{F}}^{-1}
\cdot \underline{\underline{P}}$, where $\underline{\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
(\ref{eq::J_reg_eq_cont}) 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 (\ref{eq::J_reg_eq_lin}),
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 (\ref{eq::J_reg_eq_cont}) 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
(\ref{eq::J_reg_eq_cont}), 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.

\subsection{The nonlinear discrete equilibrium gap regularization}
\label{sec::methods::nldeg}
\subsubsection{Body stresses regularization}
\label{sec::methods::nldeg::body}

In order to avoid the limitation of the regularization term
(\ref{eq::J_reg_eq_cont}), 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}}
({\underline{\underline{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{\Pi_b}$} onto a
space of square integrable functions, typically a standard continuous
finite element space denoted $V^h$:
{\begin{equation}
\underline{\Pi_b} \in V^h\; \left|\; \int_{\Omega_0} \underline{\Pi_b}
\cdot \underline{\Pi_b^*} = -\int_{\Omega_0} \underline{\text{Div}}
({\underline{\underline{P}}}) \cdot \underline{\Pi_b^*}\right.\quad
\forall \underline{\Pi_b^*} \in V^h_0,
\end{equation}}\noindent
where $V^h_0$ denotes the space of functions of $V^h$ that vanish at
the boundary, i.e., $\underline{\Pi_b^*}({\partial\Omega_0}) = 0$.
Indeed, for now we are only interested in body equilibrium gap. Thus,
after integration by parts, we obtain:
{\begin{equation}
\underline{\Pi_b} \in V^h\; \left|\; \int_{\Omega_0} \underline{\Pi_b}
\cdot \underline{\Pi_b^*} = \int_{\Omega_0}{\underline{\underline{P}}}:
{\underline{\underline{\text{Grad}}}}({\Pi_b^*})\right.\quad \forall
\underline{\Pi_b^*} \in V^h_0. \label{eq::projection_Pi_b}
\end{equation}}\noindent
Similar projections, though often performed element-by-element, are
used in \textit{a posteriori} error estimation methods, to quantify the
distance between the finite element solution and an actual equilibrium
solution, i.e., the discretization error~\cite{Ladeveze2005,
Zienkiewicz2013TFEMiBaFk}. 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:
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq}}_b({\underline{U}}) := \frac{1}{2}
\int_{\Omega_0} \underline{\Pi_b} \cdot \underline{\Pi_b}\,
\mathrm{d}\Omega_0. \label{eq::J_reg_eq_disc_b}
\end{equation}}\noindent
A key question is the calculation of this term in the finite element
context. Projection (\ref{eq::projection_Pi_b}) simply leads to the
following linear system:
{\begin{equation}
\underline{\underline{\mathbb{M}}} \cdot \underline{{\bbPi}_b} =
\underline{\mathbb{R}_b}
\end{equation}}\noindent
with $\underline{\underline{\mathbb{M}}} := \int_{\Omega_0}
\underline{\underline{\mathbb{N}}} \cdot
{}^{t}{\underline{\underline{\mathbb{N}}}}\, \mathrm{d}\Omega_0$ the
mass matrix, $\underline{{\bbPi}_b}$ such that $\underline{\Pi_b} =
{}^{t}{\underline{\underline{\mathbb{N}}}} \cdot
\underline{{\bbPi}_b}$, and
{\begin{equation}
({\mathbb{R}_b})_i := \begin{cases}
\displaystyle\int_{\Omega_0} \underline{\underline{P}} :
\underline{\underline{\text{Grad}}}({\underline{\mathbb{N}_i}})\,
\mathrm{d}\Omega_0 & \text{if}\ i\ \text{body d.o.f.}\\
0 & \text{if}\ i\ \text{boundary d.o.f.},
\end{cases}
\end{equation}}\noindent
where $\underline{\mathbb{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{\underline{\mathbb{N}}}$ array. Thus, the norm
(\ref{eq::J_reg_eq_disc_b}) can be expressed as
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq}}_b({\underline{\mathbb{U}}}) = \tfrac{1}{2}\,
{}^{t}{\underline{{\bbPi}_b}} \cdot \underline{\underline{\mathbb{M}}}
\cdot \underline{{\bbPi}_b} = \tfrac{1}{2} {}^{t}
{\underline{\mathbb{R}_b}}\cdot\underline{\underline{\mathbb{M}}}^{-1}
\cdot {\underline{\mathbb{R}_b}}. \label{eq::J_reg_eq_disc_b_final}
\end{equation}}\noindent
Note that if we linearize this expression, we obtain
$\underline{\mathbb{R}_b} \approx \underline{\underline{\mathbb{K}^*}}
\cdot {\underline{\mathbb{U}}}$ and thus $J^{\mathrm{reg},\mathrm{eq}}_b
\approx ({1}/{2}){}^{t}{\underline{\mathbb{U}}} \cdot
{{}^{t}{\underline{\underline{\mathbb{K}^*}}}} \cdot
{\underline{\underline{\mathbb{M}}}^{-1}} \cdot
{\underline{\underline{\mathbb{K}^*}}} \cdot {\underline{\mathbb{U}}}$,
i.e., an expression similar to (\ref{eq::J_reg_eq_lin}), the mass
matrix allowing to make the term{\break} consistent when refining the mesh.

\subsubsection{Boundary tractions regularization}
\label{sec::methods::nldeg::boundary}

The regularization term proposed in
Section~\ref{sec::methods::nldeg::body} 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~\cite{Leclerc2010EM}, 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~\cite{newbib}.

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}}) := {}^{t}{\underline{N}} \cdot
\underline{\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$:
{\begin{equation}
\underline{\Pi_n} \in \partial V^h\;\left|\;\int_{\partial\Omega_0}\right.
\underline{\Pi_n} \cdot \underline{\Pi_n^*} = - \int_{\partial\Omega_0}
\underline{\text{Grad}_s}{(F_n)} \cdot \underline{\Pi_n^*}\quad
\forall \underline{\Pi_n^*} \in \partial V^h,
\end{equation}}\noindent
where $\underline{\text{Grad}_s}{(F_n)} = \underline{\underline{\Pi}}
\cdot \underline{\text{Grad}}{(F_n)}$ with $\underline{\underline{\Pi}}
:= \underline{\underline{1}} - \underline{N} \otimes \underline{N}$ the
projection operator onto the domain boundary~\cite{Brandner2021}. After
integration by parts of the right hand side we obtain:
{\begin{equation}
\underline{\Pi_n}\in\partial V^h\;\left|\;\int_{\partial\Omega_0}\right.
\underline{\Pi_n} \cdot \underline{\Pi_n^*}=\int_{\partial\Omega_0}
F_n \cdot \text{Div}_s{(\underline{\Pi_n^*})}\quad \forall
\underline{\Pi_n^*} \in \partial V^h,
\end{equation}}\noindent
where $\text{Div}_s{(\underline{\Pi_n^*})} = \text{tr}
({\underline{\underline{\Pi}} \cdot \underline{\underline{\text{Grad}}}
{(\underline{\Pi_n^*})} \cdot
\underline{\underline{\Pi}}})$~\cite{Brandner2021}. Then, the
regularization term is actually defined as:
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq}}_n({\underline{U}}):= \frac{1}{2}
\int_{\partial\Omega_0} \underline{\Pi_n} \cdot \underline{\Pi_n}\,
\mathrm{d}\partial\Omega_0. \label{eq::J_reg_eq_disc_n}
\end{equation}}\noindent
The discretization procedure is similar to the bulk term
(\ref{eq::J_reg_eq_disc_b_final}):
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq}}_n({\underline{\mathbb{U}}})=\tfrac{1}{2}
{}^{t}{\underline{\mathbb{R}_n}} \cdot
\underline{\underline{\mathbb{M}_{\partial\Omega_0}^{-1}}} \cdot
\underline{\mathbb{R}_n} \label{eq::J_reg_eq_disc_n_final}
\end{equation}}\noindent
where $\underline{\underline{\mathbb{M}_{\partial\Omega_0}}} :=
\int_{\partial\Omega_0} \underline{\underline{\mathbb{N}}} \cdot
{}^{t}{\underline{\underline{\mathbb{N}}}}\, \mathrm{d}\partial
\Omega_0$ is the ``mass'' matrix of the domain boundary, and
{\begin{equation}
({\mathbb{R}_n})_i := \begin{cases}
0 & \text{if}\ i\ \text{body d.o.f.}\\
\displaystyle\int_{\partial\Omega_0} F_n \cdot
\text{Div}_s{(\underline{\mathbb{N}_i})}\, \mathrm{d}\partial
\Omega_0 & \text{if}\ i\ \text{boundary d.o.f.}
\end{cases}
\end{equation}}\noindent
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{T}} \cdot
\underline{\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~(\ref{eq::J_reg_eq_disc_n_final}):
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq}}_{t,2D}({\underline{\mathbb{U}}})=
\tfrac{1}{2} {}^{t}{\underline{\mathbb{R}_{t,2D}}} \cdot
\underline{\underline{\mathbb{M}_{\partial\Omega_0}^{-1}}} \cdot
\underline{\mathbb{R}_{t,2D}}
\end{equation}}\noindent
where
{\begin{equation} \label{eq::J_reg_eq_disc_t_2D}
({\mathbb{R}_{t,2D}})_i := \begin{cases}
0 & \text{if}\ i\ \text{body d.o.f.}\\
\displaystyle \int_{\partial\Omega_0} F_{t,2D} \cdot
\text{Div}_s{(\underline{\mathbb{N}_i})}\, \mathrm{d}\partial\Omega_0
& \text{if}\ i\ \text{boundary d.o.f.}
\end{cases}
\end{equation}}\noindent
In 3D the tangential force is a vector, defined as
$\underline{F_{t,3D}}({\underline{U}}) := \underline{\underline{\Pi}}
\cdot ({\underline{\underline{P}}({\underline{U}}) \cdot
\underline{N}})$, and we propose to simply penalize the gradient of its
norm, which leads to:
{\begin{equation}
J^{\mathrm{reg},\mathrm{eq}}_{t,3D}({\underline{\mathbb{U}}})=
\tfrac{1}{2} {}^{t}{\underline{\mathbb{R}_{t,3D}}} \cdot
\underline{\underline{\mathbb{M}_{\partial\Omega_0}^{-1}}} \cdot
\underline{\mathbb{R}_{t,3D}}
\end{equation}}\noindent
where
{\begin{equation}
({\mathbb{R}_{t,3D}})_i := \begin{cases}
0 & \text{if}\ i\ \text{body d.o.f.}\\
\displaystyle \int_{\partial\Omega_0} \|{\underline{F_{t,3D}}}\|\cdot
\text{Div}_s{(\underline{\mathbb{N}_i})}\, \mathrm{d}\partial\Omega_0
& \text{if}\ i\ \text{boundary d.o.f.}
\end{cases} \label{eq::J_reg_eq_disc_t_3D}
\end{equation}}\noindent

\subsection{Inverse problem reformulation}
\label{sec::methods::inverse_problem}

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:
{\begin{equation}
\text{Find}\ \underline{T^{\mathrm{sol}}} :=
\text{argmin}_{\{\underline{T}\}}\left\{J({\underline{T}}) :=
({1-\beta}) \|{I \circ \underline{\Phi}({\underline{T}}) - I_0}
\|_{\Omega_0}^2 + \beta \|{\underline{T}}\|_{\partial
\Omega_0}^2\right\},
\end{equation}}\noindent
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~\cite{Allaire2007a}. Here we propose to simply use
penalization, such that the problem becomes:
{\begin{eqnarray}
\text{Find}\ ({\underline{T^{\mathrm{sol}}},
\underline{U^{\mathrm{sol}}}}) &:=& \text{argmin}_{\{\underline{T},
\underline{U}\}}\left\{J({\underline{T}, \underline{U}}):=({1-\beta})
\|{I \circ \underline{\Phi} - I_0}\|_{\Omega_0}^2\right.\nonumber\\
&& +\left. \beta \|{\underline{T}}\|_{\partial\Omega_0}^2 + \gamma
\left({\|{\underline{\text{Div}}({\underline{\underline{P}}})}
\|_{\Omega_0}^2 + \|{{}^{t}{\underline{\underline{P}}} \cdot
{}^{t}{\underline{\underline{F}}}^{-1}-\underline{\underline{F}}^{-1}
\cdot \underline{\underline{P}}}\|_{\Omega_0}^2 +
\|{\underline{\underline{P}} \cdot \underline{N} - \underline{T}}
\|_{\partial\Omega_0}^2}\right)\right\},\qquad
\end{eqnarray}}\noindent
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{\underline{F}}$ {and} $\underline{\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{\underline{P}}({\underline{U}}) \cdot \underline{N}$,
leading to the following formulation:
{\begin{equation}
\text{Find}\ \underline{U^{\mathrm{sol}}} :=
\text{argmin}_{\{\underline{U}\}}\left\{J({\underline{U}}) :=
({1-\beta}) \|{I \circ \underline{\Phi} - I_0}\|_{\Omega_0}^2
+ \beta \|{\underline{\underline{P}} \cdot \underline{N}}
\|_{\partial\Omega_0}^2 + \gamma \|{\underline{\text{Div}}
{(\underline{\underline{P}})}}\|_{\Omega_0}^2\right\}.
\end{equation}}\noindent
This formulation corresponds formally (i.e., using the right norms) to
the original formulation~(\ref{eq::problem_cont}), when using
$J_b^{\mathrm{reg},\mathrm{eq}}$ for the bulk regularization and
$J_{n/t}^{\mathrm{reg},\mathrm{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
(\ref{eq::problem_cont}) with elastic regularization term
(\ref{eq::J_reg_el_cont}) 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.

\subsection{Numerical resolution} \label{sec::methods::resolution}

After describing existing and new regularization terms, we recall that
the general motion tracking problem formulation was given by
Equation~(\ref{eq::problem_disc}), with the image correlation term
$J^{\mathrm{ima}}$ given by Equation~(\ref{eq::J_ima_disc}), and where
the regularization term $J^{\mathrm{reg}}$ may take the form
$J^{\mathrm{reg},\mathrm{el}}$ (Equation~(\ref{eq::J_reg_el_cont})),
$J^{\mathrm{reg},\mathrm{eq},\mathrm{cont}}$
(Equation~(\ref{eq::J_reg_eq_cont})), $J_b^{\mathrm{reg},\mathrm{eq}}$
(Equation~(\ref{eq::J_reg_eq_disc_b})),
$J_n^{\mathrm{reg},\mathrm{eq}}$
(Equation~(\ref{eq::J_reg_eq_disc_n})),
$J_t^{\mathrm{reg},\mathrm{eq}}$
(Equation~(\ref{eq::J_reg_eq_disc_t_2D}) in 2D,
Equation~(\ref{eq::J_reg_eq_disc_t_3D}) 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~\cite{Leclerc2010EM}; 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 (\ref{eq::problem_disc}) is actually formulated as
a root finding problem:
{\begin{equation}
\underline{\mathbb{U}} \mid \underline{\nabla \mathbb{J}}
({\underline{\mathbb{U}}}) = ({1-\beta}) \underline{\nabla
\mathbb{J}^{\mathrm{ima}}}({\underline{\mathbb{U}}}) + \beta
\underline{\nabla \mathbb{J}^{\mathrm{reg}}}({\underline{\mathbb{U}}})
= 0 \label{eq::root_finding_problem}
\end{equation}}\noindent
where the gradient of the image term is simply:
{\begin{equation}
\underline{\nabla \mathbb{J}^{\mathrm{ima}}}({\underline{\mathbb{U}}})
:= \int_{\Omega_0} \left({I\left({\underline{X} +
{}^{t}{\underline{\underline{\mathbb{N}}}} \cdot
\underline{\mathbb{U}}}\right) - I_0}\right) \frac{\partial I}
{\partial\underline{x}}\left({\underline{X} +
{}^{t}{\underline{\underline{\mathbb{N}}}} \cdot
\underline{\mathbb{U}}}\right)\cdot\underline{\underline{\mathbb{N}}}
\,\mathrm{d}\Omega_0,
\end{equation}}\noindent
and the gradients of the proposed regularization terms are:
{\begin{equation}
\underline{\nabla \mathbb{J}^{\mathrm{reg},\mathrm{eq}}_b}
({\underline{\mathbb{U}}}) :=
{}^{t}{\underline{\underline{\mathrm{d}\mathbb{R}_b}}}
({\underline{\mathbb{U}}}) \cdot
\underline{\underline{\mathbb{M}}}^{-1} \cdot \underline{\mathbb{R}_b}
({\underline{\mathbb{U}}})
\end{equation}}\noindent
with
{\begin{equation}
({\mathrm{d}\mathbb{R}_b})_{ij} := \begin{cases}
\displaystyle \int_{\Omega_0} \underline{\underline{\text{Grad}}}
({\underline{\mathbb{N}_i}}) : \frac{\partial\underline{\underline{P}}}
{\partial\underline{\underline{F}}} : \underline{\underline{\text{Grad}}}
({\underline{\mathbb{N}_j}})\, \mathrm{d}\Omega_0 &
\text{if}\ i\ \text{body d.o.f.}\\
0 & \text{if}\ i\ \text{boundary d.o.f.},
\end{cases}
\end{equation}}\noindent
and
{\begin{equation}
\underline{\nabla \mathbb{J}^{\mathrm{reg},\mathrm{eq}}_{n/t}}
({\underline{\mathbb{U}}}) :=
{}^{t}{\underline{\underline{\mathrm{d}\mathbb{R}_{n/t}}}}
({\underline{\mathbb{U}}}) \cdot
\underline{\underline{\mathbb{M}_{\partial\Omega_0}^{-1}}} \cdot
\underline{\mathbb{R}_{n/t}}({\underline{\mathbb{U}}})
\end{equation}}\noindent
with
{\begin{equation}
({\mathrm{d}\mathbb{R}_{n/t}})_{ij} := \begin{cases}
0 & \text{if}\ i\ \text{body d.o.f.}\\
\displaystyle \int_{\partial\Omega_0} \text{Div}_s
{(\underline{\mathbb{N}_i})} \frac{\partial F_{n/t}}
{\partial\underline{\underline{F}}} :
\underline{\underline{\text{Grad}}}{(\underline{\mathbb{N}_j})}\,
\mathrm{d}\partial\Omega_0 & \text{if}\ i\ \text{boundary d.o.f.}
\end{cases}
\end{equation}}\noindent
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~\cite{Passieux2019IJNME}.
Details were given in~\cite{Genet2018MIA} in the context of the
regularization term $J^{\mathrm{reg},\mathrm{eq},\mathrm{cont}}$
(Equation~(\ref{eq::J_reg_eq_cont})). 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:
{\begin{equation} \label{eq::linear_problem}
\underline{\Delta \mathbb{U}} \mid \left({({1-\beta})
\underline{\underline{\nabla\nabla \mathbb{J}^{\mathrm{ima}}}} + \beta
\underline{\underline{\nabla\nabla \mathbb{J}^{\mathrm{reg}}}}}\right) \cdot
\underline{\Delta \mathbb{U}} = -\left({({1-\beta})
\underline{\nabla \mathbb{J}^{\mathrm{ima}}} + \beta
\underline{\nabla \mathbb{J}^{\mathrm{reg}}}}\right),
\end{equation}}\noindent
where
{\begin{equation}
\underline{\underline{\nabla\nabla \mathbb{J}^{\mathrm{ima}}}} :=
\int_{\Omega_0} \underline{\underline{\mathbb{N}}} \cdot
\left({\frac{\partial I}{\partial \underline{x}} \cdot 
{\frac{\mbox{}^{t}\partial I}{\partial \underline{x}}}}\right) \cdot
{}^{t}{\underline{\underline{\mathbb{N}}}}\, \mathrm{d}\Omega_0,
\end{equation}}\noindent
which corresponds to a mass matrix weighted by the deformed image
gradients, and
{\begin{equation}
\underline{\underline{\nabla\nabla \mathbb{J}^{\mathrm{reg},
\mathrm{eq}}_{b/n/t}}} :=
{}^{t}{\underline{\underline{\mathrm{d}\mathbb{R}_{b/n/t}}}} \cdot
\underline{\underline{\mathbb{M}_{\Omega_0/\partial\Omega_0}^{-1}}}
\cdot \underline{\underline{\mathrm{d}\mathbb{R}_{b/n/t}}}.
\end{equation}}\noindent
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:
{\begin{equation}
\text{err} := \frac{\|{\underline{\Delta\mathbb{U}}}\|}
{\|{\underline{\mathbb{U}}}\|} < \text{tol}.
\end{equation}}\noindent
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~\cite{Ogden1972PRSLA, Ciarlet1982CRASSI}), 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~\cite{LeTallec1994HoNA, Genet2019CB}. 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~\cite{Press2007} that prevents the passing of
energy barriers.

Finally, the image series integrator is described in
Algorithm~\algOneREF{alg::image_integrator}, the nonlinear solver in
Algorithm~\algTwoREF{alg::nonlinear_solver}, and the line search in
Algorithm~\algThreeREF{alg::line_search}.

\beqimg\begin{eqimg*}
\begin{algorithm}[htp]
 \SetKwBlock{Initialisation}{Initialisation}{}
 \Initialisation
 {
 Read initial image $I_0$ \\
 Instantiate finite element solution: $\underline{\mathbb{U}}
 \leftarrow 0$ \\
 }
 \ForEach{\upshape frame index $t = 1, 2, \ldots$}
 {
 Read current image $I_t$ \\
 Compute current displacement: $\underline{\mathbb{U}}
 \leftarrow \texttt{nonlinear\_solver}({I_0, I_t})$
 }
\caption{Image series integrator. The \texttt{nonlinear\_solver} is
detailed Algorithm~\algTwoREF{alg::nonlinear_solver}. 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~\algFourREF{alg::image_integrator_with_refine}.}
\label{alg::image_integrator}
\end{algorithm}
\eeqimg\end{eqimg*}

\beqimg\begin{eqimg*}
\begin{algorithm}[htp]
 \While{$\mathrm{err} > \mathrm{tol}$}
 {
 Assemble residual, Jacobian: $\underline{\nabla \mathbb{J}}$,
 $\underline{\underline{\nabla\nabla \mathbb{J}}}$\\
 Compute solution increment: $\underline{\Delta\mathbb{U}}
 \leftarrow \texttt{linear\_solver}({\underline{\nabla \mathbb{J}},
 \underline{\underline{\nabla\nabla \mathbb{J}}}})$\\
 Update error: $\text{err} \leftarrow
 \frac{\|{\underline{\Delta\mathbb{U}}}\|}
 {\|{\underline{\mathbb{U}}}\|}$\\
 Compute relaxation: $\alpha \leftarrow \texttt{line\_search}
 ({\underline{\mathbb{U}}, \underline{\Delta\mathbb{U}}})$\\
 Update solution: $\underline{\mathbb{U}} \leftarrow
 \underline{\mathbb{U}} + \alpha~\underline{\Delta\mathbb{U}}$
 }
\caption{Nonlinear solver. In practice we use a tolerance of
$\text{tol} = 1\%$.} \label{alg::nonlinear_solver}
\end{algorithm}
\eeqimg\end{eqimg*}

\beqimg\begin{eqimg*}
\begin{algorithm}[htp]
 \SetKwBlock{Initialisation}{Initialisation}{}
 \Initialisation
 {
 Compute initial energy: $J_0\leftarrow J({\underline{\mathbb{U}}})$\\
 Initialise relaxation counter: $k \leftarrow 0$
 }
 \While{$J > J_0$}
 {
 Define current relaxation: $\alpha \leftarrow f^k$\\
 Compute current energy: $J \leftarrow J({\underline{\mathbb{U}}
 + \alpha~\underline{\Delta\mathbb{U}}})$\\
 Update relaxation counter: $k \leftarrow k+1$
 }
\caption{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~(\ref{eq::problem_disc}).}} \label{alg::line_search}
\end{algorithm}
\eeqimg\end{eqimg*}

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 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~\cite{Bornert2009EM, Leclerc2010EM}, 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~\algOneREF{alg::image_integrator}. 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~\algFourREF{alg::image_integrator_with_refine}.
{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~\cite{Fedele2013IJNME}; however, for the sake of
simplicity, here we use the full images at all mesh refinement levels.}

\beqimg\begin{eqimg*}
\begin{algorithm}[t]
 \SetKwBlock{Initialisation}{Initialisation}{}
 \ForEach{\upshape refinement level $k = 1, 2, \ldots$}
 {
 \Initialisation
 {
 Read initial image $I_0$\\
 Instantiate finite element solution $\underline{\mathbb{U}^k}$\\
 }
 \ForEach{\upshape frame index $t = 1, 2, \ldots$}
 {
 Read current image $I_t$\\
 \uIf{$k > 1$}
 {
 Initialize displacement: $\underline{\mathbb{U}^k}
 \leftarrow \underline{\mathbb{U}_t^{k-1}}$
 }
 Compute current displacement: $\underline{\mathbb{U}^k}
 \leftarrow \texttt{nonlinear\_solver}({I_0, I_t})$\\
 Save displacement: $\underline{\mathbb{U}_t^k} \leftarrow
 \underline{\mathbb{U}^{k}}$
 }
 }
\caption{Image series integrator with multi-resolution. Compared to the
single-refinement-level Algorithm~\algOneREF{alg::image_integrator}, 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$.}
\label{alg::image_integrator_with_refine}
\end{algorithm}
\eeqimg\end{eqimg*}

These algorithms have been implemented in an open-source
library~\cite{dolfin_warp} written in python and based on the
FEniCS~\cite{Logg2012, Alnaes2015ANS} and VTK~\cite{Schroeder2006}
libraries. It is currently freely available online at
\url{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~\cite{N-DEG-paper-demos}: static versions are given in the
Supplementary material of the paper while interactive versions are currently
available online at
\url{https://mgenet.gitlabpages.inria.fr/N-DEG-paper-demos/index.html}.

\subsection{Synthetic data} \label{sec::methods::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~\cite{Berberoglu2022MIA}), and textured (with used a
tagged-MRI-like pattern, see~\cite{Rutz2008MRM}) images. The images
occupy the spatial domain ${0\underline{;}1}^2$ (arbitrary unit), the
temporal domain ${0\underline{;}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\underline{;}\ 0.7} \times
{0.2\underline{;}\ 0.8}$ for the translation case and $\Omega_0^{r,c,s}
:= {0.2\underline{;}\ 0.8}^2$ for the rotation, compression {and} shear
cases. The motion {models} are given by
{\begin{equation}
\begin{cases}
\underline{x^t}(\underline{X},t) := \underline{X}+t~\underline{D} &
\text{with}\ \underline{D} := \left(\begin{smallmatrix}
0.2\\ 0
\end{smallmatrix}\right)\\
\underline{x^r}(\underline{X},t) := \underline{X_0} +
\underline{\underline{R}} \cdot ({\underline{X}-\underline{X_0}}) &
\text{with}\ \underline{X_0} := \left(\begin{smallmatrix}
0.5 \\ 0.5
\end{smallmatrix}\right),\ \underline{\underline{R}} :=
\left(\begin{smallmatrix}
{+}\text{cos}({\theta}) & {-}\text{sin}({\theta})\\
{+}\text{sin}({\theta}) & {+}\text{cos}({\theta})
\end{smallmatrix}\right), \ \theta = \frac{t\uppi}{4}
\vspace*{.2pc}\\
\underline{x^c}(\underline{X},t) := \underline{X_0} +
\underline{\underline{F^c}} \cdot ({\underline{X}-\underline{X_0}}) &
\text{with}\ \underline{X_0} := \left(\begin{smallmatrix}
0.5 \\ 0.5
\end{smallmatrix}\right), \ \underline{\underline{F^c}} :=
\left(\begin{smallmatrix}
\sqrt{1-2t 0.2} & 0\\
0 & 1
\end{smallmatrix}\right)\\
\underline{x^s}(\underline{X},t) := \underline{X_0} +
\underline{\underline{F^s}} \cdot ({\underline{X}-\underline{X_0}}) &
\text{with}\ \underline{X_0} := \left(\begin{smallmatrix}
0.5 \\ 0.5
\end{smallmatrix}\right), \ \underline{\underline{F^s}} :=
\left(\begin{smallmatrix}
1 & 0.2 \\
0 & 1
\end{smallmatrix}\right),
\end{cases} \label{eq::mappings}
\end{equation}}\noindent
which represent a 0.2 (arbitrary unit) translation, a ${\uppi}/{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
$\underline{x}$), the initial (i.e., first time frame) position of the
pixel $(\underline{X}(\underline{x},t))$ was determined through
inversion of the mappings (\ref{eq::mappings}) and the intensity was
given by
{\begin{equation}
I(\underline{x},t) := I_0({\underline{X}(\underline{x},t)})\quad
\text{with}\ I_0(\underline{X}) := \sqrt{\left|{\sin\left({\frac{\uppi
X}{s}}\right)}\right|\ \left|{\sin\left({\frac{\uppi Y}{s}}\right)}
\right|}. \label{eq::texture_model}
\end{equation}}\noindent
Note that a more complex imaging model~\cite{Berberoglu2022MIA} could
have been considered as well.

For the cardiac-like case, a ring (center $\underline{X_0} :=
\left(\begin{smallmatrix}0.5 \\ 0.5\end{smallmatrix}\right)$, 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 $-{\uppi}/{4}$~rad) and external
(inward displacement of 0.05 and rotation of $-{\uppi}/{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~$\underline{x}$) the
reference position of the pixel was computed as
$\underline{X}(\underline{x},t) = \underline{x} - \underline{U}
(\underline{x}, t)$. Finally, the texture model
(\ref{eq::texture_model}) 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~\cite{Berberoglu2022MIA}, could have been considered.

\section{Results {and} discussion} \label{sec::results}

To establish the tracking performance of the proposed method (described
in Section~\ref{sec::methods::nldeg}), we now present tracking results
on various synthetic images (described in
Section~\ref{sec::methods::data}), for various regularization terms,
namely elastic (Equation~(\ref{eq::J_reg_el_cont}), continuous version
of the equilibrium gap (Equation~(\ref{eq::J_reg_eq_cont})), discrete
version of the equilibrium gap (Equation~(\ref{eq::J_reg_eq_disc_b}))
and discrete version of the nonlinear equilibrium gap including surface
traction regularization terms (Equations~(\ref{eq::J_reg_eq_disc_b}),
(\ref{eq::J_reg_eq_disc_n}), (\ref{eq::J_reg_eq_disc_t_2D}) and/or
(\ref{eq::J_reg_eq_disc_t_3D}))). 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
{\begin{equation}
\text{err} := \frac{\sqrt{\displaystyle\frac{1}{T}\int_0^T \frac{1}
{|{\Omega_0}|} \int_{\Omega_0} \|{\underline{U}-
\underline{U^{\mathrm{ex}}}}\|^2}}{\displaystyle\sqrt{\frac{1}{T}
\int_0^T \frac{1}{|{\Omega_0}|} \int_{\Omega_0}
\|{\underline{U^{\mathrm{ex}}}}\|^2}}, \label{eq::error}
\end{equation}}\noindent
where $\underline{U}$ is the tracked displacement {and}
$\underline{U^{\mathrm{ex}}}$ is the ground truth used to generate the
images.

\subsection{Rigid body motion} \label{sec::results::rigid}

We start with simple rigid motions, namely pure translation and pure
rotation. Figures~\ref{fig1} and~\ref{fig2} show normalized tracking
error as a function of regularization strength $\beta \in
{0\underline{;}1}$, for various levels of image SNR and for various
regularization terms. In the plots, dots represent noise realizations,
and lines represent their average.

{
\begin{figure}
{\vspace*{-2pt}}
\includegraphics{fig01a}
{\vspace*{-.5pc}}
\caption{\label{fig1}{\fontsize{9.3}{10}\selectfont
Translation case. (a)~Normalized displacement
error (\ref{eq::error}) as a function of regularization strength
(parameter~$\beta$), for various levels of image noise, for various
regularization terms (first row: $J^{\mathrm{reg}, \mathrm{el}}$
(\ref{eq::J_reg_el_cont}); second row: $J^{\mathrm{reg}, \mathrm{eq},
\mathrm{cont}}$ (\ref{eq::J_reg_eq_cont}); third row: $J^{\mathrm{reg},
\mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}); fourth row:
$J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}),
$J^{\mathrm{reg}, \mathrm{eq}}_n$ (\ref{eq::J_reg_eq_disc_n}) and
$J^{\mathrm{reg}, \mathrm{eq}}_{t,2D}$ (\ref{eq::J_reg_eq_disc_t_2D})),
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.)}
\end{figure}
\def\thefigure{1 (cont.)}
\begin{figure}
\ContinuedFloat
{\vspace*{.3pc}}
\includegraphics{fig01}
{\vspace*{.3pc}}
\caption{ 
%%%
(b)~Tracking solutions for the
discrete version of the nonlinear equilibrium gap
(\ref{eq::J_reg_eq_disc_b}) including surface traction regularization
(\ref{eq::J_reg_eq_disc_n}) and (\ref{eq::J_reg_eq_disc_t_2D}), 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.}
{\vspace*{.3pc}}
\end{figure}
}

{
\begin{figure}
{\vspace*{-.3pc}}
\includegraphics{fig02a}
{\vspace*{-6pt}}
\caption{\label{fig2}\fontsize{9.8}{11.6}\selectfont Rotation case.
(a)~Normalized displacement error (\ref{eq::error}) as a function of
regularization strength (parameter~$\beta$), for various levels of
image noise, for various regularization terms (first row:
$J^{\mathrm{reg}, \mathrm{el}}$ (\ref{eq::J_reg_el_cont}); second row:
$J^{\mathrm{reg}, \mathrm{eq}, \mathrm{cont}}$
(\ref{eq::J_reg_eq_cont}); third row: $J^{\mathrm{reg}, \mathrm{eq}}_b$
(\ref{eq::J_reg_eq_disc_b}); fourth row: $J^{\mathrm{reg},
\mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}), $J^{\mathrm{reg},
\mathrm{eq}}_n$ (\ref{eq::J_reg_eq_disc_n}) and $J^{\mathrm{reg},
\mathrm{eq}}_{t,2D}$ (\ref{eq::J_reg_eq_disc_t_2D})), 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.)}
\end{figure}
\def\thefigure{2 (cont.)}
\begin{figure}
\ContinuedFloat
{\vspace*{.5pc}}
\includegraphics{fig02b}
{\vspace*{.5pc}}
\caption{ 
%%%
(b)~Tracking
solutions for the discrete version of the nonlinear equilibrium gap
(\ref{eq::J_reg_eq_disc_b}) including surface traction regularization
(\ref{eq::J_reg_eq_disc_n}) and (\ref{eq::J_reg_eq_disc_t_2D}), 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.)}
{\vspace*{.5pc}}
\end{figure}
\def\thefigure{2 (cont.)}
\begin{figure}
\ContinuedFloat
{\vspace*{.3pc}}
\includegraphics{fig02}
{\vspace*{.3pc}}
\caption{ 
%%%
(c)~Tracking solutions for the linear elastic regularization
term (\ref{eq::J_reg_el_cont}), 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.}
{\vspace*{.3pc}}
\end{figure}}

Let us first discuss the pure translation case, i.e.,
Figure~\ref{fig1}. 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~\ref{fig1}a) and
continuous equilibrium gap (second row in Figure~\ref{fig1}a) terms
behave better than the discrete equilibrium gap terms (third and fourth
rows in Figure~\ref{fig1}a). 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~\ref{fig1}a), 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, 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~\ref{fig1}a), 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~\ref{fig1}a).

Let us now discuss the pure rotation case, i.e., Figure~\ref{fig2}. For
the nonlinear elastic (upper right in Figure~\ref{fig2}a) and all
equilibrated (second, third and fourth rows in Figure~\ref{fig2}a)
regularization terms, the conclusions are the same as for the pure
translation case. For the linear elastic regularization term (upper
left in Figure~\ref{fig2}a), 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~\ref{fig2}c). This is not the case for
hyperelastic regularization (upper right in Figure~\ref{fig2}a), 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~\ref{fig2}a), 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.

\subsection{Nonrigid homogeneous deformation}
\label{sec::results::nonrigid}

We now consider nonrigid---though still homogeneous---deformations,
namely compression (Figure~\ref{fig3}) and shear (Figure~\ref{fig4}).
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~\ref{fig3}a and~\ref{fig4}a),
which perform very well. Regularization terms based on the elastic
energy (i.e.,~first row in Figures~\ref{fig3}a and~\ref{fig4}a) 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.

{
\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig03a}
{\vspace*{-6pt}}
\caption{\label{fig3}\fontsize{9.5}{11.6}\selectfont 
Compression case. (a)~Normalized
displacement error (\ref{eq::error}) as a function of regularization
strength (parameter~$\beta$), for various levels of image noise, for
various regularization terms (first row: $J^{\mathrm{reg},
\mathrm{el}}$ (\ref{eq::J_reg_el_cont}); second row: $J^{\mathrm{reg},
\mathrm{eq}, \mathrm{cont}}$ (\ref{eq::J_reg_eq_cont}); third row:
$J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}); fourth
row: $J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}) and
$J^{\mathrm{reg}, \mathrm{eq}}_{t,2D}$ (\ref{eq::J_reg_eq_disc_t_2D})),
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.
%%%
\break (Continued on next page.)}
\end{figure}
\def\thefigure{3 (cont.)}
\begin{figure}
\ContinuedFloat
{\vspace*{6pt}}
\includegraphics{fig03b}
{\vspace*{6pt}}
\caption{ 
%%%
(b)~Tracking
solutions for the discrete version of the nonlinear equilibrium gap
(\ref{eq::J_reg_eq_disc_b}) including surface traction regularization
(\ref{eq::J_reg_eq_disc_t_2D}), 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.)}
{\vspace*{8pt}}
\end{figure}
\def\thefigure{3 (cont.)}
\begin{figure}
\ContinuedFloat
\includegraphics{fig03}
\caption{ 
%%%
(c)~Tracking solutions for the
elastic regularization term (\ref{eq::J_reg_el_cont}), 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 tends to prevent the mesh from deforming, which
degrades the tracking.}
{\vspace*{6pt}}
\end{figure}}

{
\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig04a}
{\vspace*{-4pt}}
\caption{\label{fig4}\fontsize{9.6}{11.6}\selectfont 
Shear case. (a)~Normalized
displacement error (\ref{eq::error}) as a function of regularization
strength (parameter~$\beta$), for various levels of image noise, for
various regularization terms (first row: $J^{\mathrm{reg},
\mathrm{el}}$ (\ref{eq::J_reg_el_cont}); second row: $J^{\mathrm{reg},
\mathrm{eq}, \mathrm{cont}}$ (\ref{eq::J_reg_eq_cont}); third row:
$J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}); fourth
row: $J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}) and
$J^{\mathrm{reg}, \mathrm{eq}}_n$ (\ref{eq::J_reg_eq_disc_n})), 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.)}
\end{figure}
\def\thefigure{4 (cont.)}
\begin{figure}
\ContinuedFloat
\includegraphics{fig04b}
\caption{ 
%%%
(b)~Tracking solutions for the discrete version of
the nonlinear equilibrium gap (\ref{eq::J_reg_eq_disc_b}) including
surface traction regularization (\ref{eq::J_reg_eq_disc_n}), 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.)}
\end{figure}
\def\thefigure{4 (cont.)}
\begin{figure}
\ContinuedFloat
{\vspace*{6pt}}
\includegraphics{fig04}
{\vspace*{6pt}}
\caption{ 
%%%
(c)~Tracking solutions for the elastic regularization term
(\ref{eq::J_reg_el_cont}), 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 tends to prevent the mesh
from deforming, which degrades the tracking.}
{\vspace*{6pt}}
\end{figure}}

One can also notice that the continuous formulation of the equilibrium
gap regularization (second row in Figures~\ref{fig3}a and~\ref{fig4}a)
seems to perform better than the discrete formulation (third and fourth
rows in Figures~\ref{fig3}a and~\ref{fig4}a); as we will see later,
this is actually due to the fact that $J^{\mathrm{reg}, \mathrm{eq},
\mathrm{cont}}$ (\ref{eq::J_reg_eq_cont}) 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\break cases.

It is important to notice also that for the compression case, only the
tangential term (\ref{eq::J_reg_eq_disc_t_2D}) was included in the
boundary traction terms (fourth row in Figure~\ref{fig3}a). 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
(\ref{eq::J_reg_eq_disc_n}) was included in the boundary tractions
terms (fourth row in Figure~\ref{fig4}a).

\subsection{Nonrigid heterogeneous deformation}
\label{sec::results::heart}

Let us now consider more realistic images, mimicking cardiac tagged
magnetic resonance imaging slices. Tracking results are presented in
Figure~\ref{fig5}. Elastic regularization (first row of
Figure~\ref{fig5}a) 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~\ref{fig5}a) 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~\ref{sec::results::rigid} 
and~\ref{sec::results::nonrigid}, but not in more realistic cases. The
discrete version of the equilibrium gap principle (third and fourth
rows of Figure~\ref{fig5}a), 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 stain framework ant the
Hooke law). Also, the surface traction regularization terms (fourth row
in Figure~\ref{fig5}a) help {reducing} the tracking error mean and
dispersion compared to the tracking with the bulk term only (third row
in\break Figure~\ref{fig5}a).

{
\begin{figure}
{\vspace*{-2pt}}
\includegraphics{fig05a}
{\vspace*{-2pt}}
\caption{\label{fig5}\fontsize{9.4}{11.6}\selectfont 
Cardiac-like case. (a)~Normalized
displacement error (\ref{eq::error}) as a function of regularization
strength (parameter $\beta$), for various levels of image noise, for
various regularization terms (first row: $J^{\mathrm{reg},
\mathrm{el}}$ (\ref{eq::J_reg_el_cont}); second row: $J^{\mathrm{reg},
\mathrm{eq}, \mathrm{cont}}$ (\ref{eq::J_reg_eq_cont}); third row:
$J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}); fourth
row: $J^{\mathrm{reg}, \mathrm{eq}}_b$ (\ref{eq::J_reg_eq_disc_b}),
$J^{\mathrm{reg}, \mathrm{eq}}_n$ (\ref{eq::J_reg_eq_disc_n}) and
$J^{\mathrm{reg}, \mathrm{eq}}_{t,2D}$ (\ref{eq::J_reg_eq_disc_t_2D})),
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.)}
\end{figure}
\def\thefigure{5 (cont.)}
\begin{figure}
\ContinuedFloat
{\vspace*{6pt}}
\includegraphics{fig05b}
{\vspace*{6pt}}
\caption{ 
%%%
(b)~Tracking
solutions for the discrete version of the nonlinear equilibrium gap
(\ref{eq::J_reg_eq_disc_b}) including surface traction regularization
terms (\ref{eq::J_reg_eq_disc_n}) and (\ref{eq::J_reg_eq_disc_t_2D}),
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.)}
\end{figure}
\def\thefigure{5 (cont.)}
\begin{figure}
\ContinuedFloat
\includegraphics{fig05}
\caption{ 
%%%
(c)~Tracking solutions for the continuous version of the
nonlinear equilibrium gap regularization term
(\ref{eq::J_reg_eq_cont}), 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.}
{\vspace*{6pt}}
\end{figure}}

\subsection{Impact of mesh size} \label{sec::results::refine}

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~\ref{fig6} 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~\ref{fig6}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~\ref{fig6}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{\break} 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.

\begin{figure}
\includegraphics{fig06}
\bcaption{\label{fig6}}{5}{Cardiac-like case, impact of mesh size.
(a)~Normalized displacement error (\ref{eq::error}) as a function of
mesh size, for various levels of image noise, for the discrete version
of the nonlinear equilibrium gap (\ref{eq::J_reg_eq_disc_b}) including
surface traction regularization terms (\ref{eq::J_reg_eq_disc_n}) and
(\ref{eq::J_reg_eq_disc_t_2D}), with various regularization strength.
(b)~Tracking solutions for the discrete version of the nonlinear
equilibrium gap (\ref{eq::J_reg_eq_disc_b}) including surface traction
regularization terms (\ref{eq::J_reg_eq_disc_n}) and
(\ref{eq::J_reg_eq_disc_t_2D}), with regularization strength $\beta =
0.1$.}
\end{figure}

\section{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 heterogenous
motion, which is the most representative of actual applications of
motion tracking~\cite{Genet2014JAP, Finsberg2019AJPCP}, 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 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~\cite{dolfin_warp} (currently at
\url{https://gitlab.inria.fr/mgenet/dolfin\_warp}), and the specific
code to reproduce the results of this paper is also freely
available~\cite{N-DEG-paper-demos} (in the supplementary file and
currently at
\url{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~\cite{Patte2022JBE, Laville2023BMM}. 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~\cite{Mathieu2015EM}. 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~\cite{Claire2004IJNME}, but in large 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~\cite{Avril2008EM}.

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~\cite{Leiner2019JCMR,
Friedrich2021EHJDH, Koehler2022SAaCMotHRaCCP, Lopez2022}.

An important perspective of our motion tracking algorithm is toward low
resolution images, which have been proven to drastically impact the
tracking quality~\cite{Berberoglu2021PO}. One option is to combine
images from multiple modalities and combine the strength of both
images~\cite{Berberoglu2022MIA}. Another option to consider is to
introduce a model of the imaging modality (e.g.,
MRI~\cite{Skardova2019V2}), 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.

\section*{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.

\section*{Acknowledgements}

MG would like to thank Philippe Moireau, Patrick Le Tallec {and} Pierre
Kerfriden for fruitful discussions.
\vspace*{-.4pc}

\back{}

\section*{Supplementary data}

Supporting information for this article is available on the journal's
website under \printDOI\ or from the author. 

\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{crmeca-228-suppl.pdf}}

\bibliographystyle{crunsrt}
\bibliography{crmeca20230521}
\refinput{crmeca20230521-reference.tex}

\end{document}
