A Multiscale Deformation Representation*

Noémie Debroux†, Carole Le Guyader‡, and Luminita A. Vese §

**Abstract.** Motivated by Tadmor et al.'s work ([31]) dedicated to multiscale image representation using hierarchical ($BV, L^2$) decompositions, we propose transposing their approach to the case of registration, task which consists in determining a smooth deformation aligning the salient constituents visible in an image into their counterpart in another. The underlying goal is to obtain a hierarchical decomposition of the deformation in the form of a composition of intermediate deformations: the coarser one, computed from versions of the two images capturing the essential features, encodes the main structural/geometrical deformation, while iterating the procedure and refining the versions of the two images yields more accurate deformations that map faithfully small-scale features. The proposed model falls within the framework of variational methods and hyperelasticity by viewing the shapes to be matched as Ogden materials. The material behaviour is described by means of a specifically tailored strain energy density function, complemented by $L^\infty$-penalisations ensuring that the computed deformation is a bi-Lipschitz homeomorphism. Theoretical results emphasising the mathematical soundness of the model are provided, among which the existence of minimisers/asymptotic results, and a suitable numerical algorithm is supplied, along with numerical simulations demonstrating the ability of the model to produce accurate hierarchical representations of deformations.

A very preliminary version of this work has been accepted for publication in the Eighth International Conference on Scale Space and Variational Methods in Computer Vision, 2021 ([14]) but it does not include all the theoretical results, nor the detailed related proofs. A more complete and detailed analysis of the numerical experiments is also provided. The theoretical analysis of the numerical algorithm (introduced in Section 3 and which is a result in itself) will be the subject of a separate article in preparation ([13]).

**Key words.** Multiscale representation, hyperelasticity, Ogden materials, bi-Lipschitz homeomorphisms, asymptotic results

**AMS subject classifications.** 68U10, 49, 65D18, 28

# 1. Introduction.

**1.1. Motivations.** The grey-level representation of a real scene, assumed to be an $L^2$ observation, encompasses scale-varying and noticeable objects, whether it be edges —well-identified within the small subclass of functions of bounded variation ($BV$) —, homogeneous regions, or oscillating patterns/texture, these latter features requiring the introduction of more involved intermediate spaces. Medical images for instance exemplify this multiscale structure: they often comprise structural organs irrigated by finer blood

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>**Funding:** L.A. Vese acknowledges support from the National Science Foundation under Grant # 2012868. This project was co-financed by the European Union with the European regional development fund (ERDF, 18P03390/18E01750/18P02733), by the Haute-Normandie Régional Council via the M2SINUM project and by the French Research National Agency ANR via AAP CE23 MEDISEG ANR project.</footnote>

<footnote>†Université Clermont Auvergne, CNRS, SIGMA Clermont, Institut Pascal (noemie.debroux@uca.fr).</footnote>

<footnote>‡Normandie Univ, Institut National des Sciences Appliquées de Rouen, Laboratory of Mathematics (carole.le-guyader@insa-rouen.fr)</footnote>

<footnote>§Department of Mathematics, University of California Los Angeles (lvese@math.ucla.edu).</footnote>
vessels, and possibly a noise degradation. The purpose of multiscale representation is thus to quantify accurately these subclasses lying in between the rougher space $L^2$ and the smaller space $BV$.

A special instance of such an algorithm is [31]. This latter is the foundation of our method as will be seen later. Taking $X$ as the larger functional space $L^2$, and $Y$, the smaller space $BV$, the authors in [31] assess how accurately an $L^2$-object can be approximated by its $BV$-characteristics, this being quantified by means of the family of functionals

$$J(f, \lambda) = \inf_{u+v=f} \{\lambda \|v\|_{L^2}^2 + \|u\|_{BV}\},$$

with increasing $\lambda$'s. (To understand better how TV regularisation and more precisely, how the amount of regularisation applied to the image, encodes scale of individual features, we refer the reader to [30]. Note however that the focus is shifted in [30] since the tuning parameter weights the TV semi-norm). Component $u$ of the decomposition captures the geometrical/structural features of the observation $f$, while component $v$ encodes textures and oscillatory patterns, the degree of detail fineness in component $v$ (or symmetrically, the level of coarseness of $u$) being dictated by scale parameter $\lambda$: the larger this parameter is, the fewer details the $v$ component contains. This latter observation reflects the fact that the discrimination between these two components is scale-dependent: what is viewed as texture at a fixed scale, will be part of the structural component at a more refined scale (higher $\lambda$).

An iterative dyadic refinement scheme is applied, with $\lambda_0$ a given initial scale, and reads as:

$$f = u_0 + v_0, \quad [u_0, v_0] = \arg \min_{u+v=f} J(f, \lambda_0),$$

$$v_j = u_{j+1} + v_{j+1}, \quad [u_{j+1}, v_{j+1}] = \arg \min_{u+v=v_j} J(v_j, \lambda_0 2^{j+1}), \quad j = 0, 1, \dots,$$

producing at the end of the $k^{th}$ step, the following hierarchical decomposition:

$$f = u_0 + v_0 = u_0 + u_1 + v_1 = \dots = u_0 + u_1 + \dots + u_k + v_k,$$

the dyadic blocks $u_j = u_j(f)$'s encoding different scales and resolving finer edges, while $v_k$ being thought of as a residual in the approximation of $f$ by $\sum_{j=0}^k u_j$. With an additional slight amount of smoothness on $f$, a strong $L^2$-convergence result of $\sum_{j=0}^k u_j$ towards $f$ can be established ([31, Theorem 2.2]).

## 1.2. Contributions.
Equipped with this material, we now focus on the core of the contribution which aims to transpose this idea of multiscale representation of an image to the multiscale representation of a deformation pairing two images. The underlying goal is twofold:

* (i) obtaining a hierarchical expression of the sought deformation in the form of a composition $\varphi_0 \circ \dots \circ \varphi_k \circ \dots \circ \varphi_n$, $\varphi_0$ encoding the main geometry-driven/structural
deformation, while the $\varphi_0 \circ \cdots \circ \varphi_k$'s capture more refined deformations —in comparison to [31], the composition of deformations is now a substitute for the sum of the scale-varying constituents $u_j$'s —. This enables one to dissociate the main deformation from the more localised displacements, and to handle more accurately the different levels of granularity of the deformation. In others words, it allows to separate the global deformations of the main features forming the image from the more refined and local deformations of smaller items appearing at higher scales. In a medical context for instance, the method disassociates the organ matching from the vessel one, while isolating the noise, and one can thus answer to questions like: do blood vessels have an inherent movement besides the movement induced by the organ in which they are lying?

• (ii) opening the way to a posteriori analyses, like disclosing the hidden structure of a deformation, or deriving some statistics such as mean deformations or prevailing dynamics. Precisely, a line of research in medical image analysis could consist of constructing an atlas, i.e. a mean representative of a collection of images, to estimate variability of shapes inside a population and to understand how structural changes may affect health (refer to [11] for instance). This involves identifying significant shape constituents of the set of images (possibly neglecting very small scale details that are not necessary) and mapping this group of images to an unknown mean image with the desired level of details, which could be achieved using our multiscale registration method (with the mean image as unknown in addition to the deformation mappings of each image of the cohort to this average image). The added value of the proposed model is that one can control more finely the level of detail that should be encoded in the reconstructed mean representative. A statistical analysis on the obtained deformations in order to retrieve the main modes of variations in terms of geometric distortions in the initial set of images could then be conducted at each scale : the first step would yield a mean version of the global/structural deformation, while with increasing $k$, mean representatives of more refined deformations could be generated.

To this end, an iterative refinement scheme, straightforwardly connected to [31] is introduced, incorporating finer details on both images at each step to increase the recovered deformation accuracy.

As in [31], the multiscale model is built on a functional, indifferently called parent functional or generating functional, from which the successive minimisation subproblems are designed. This parent functional is tailored to comply with some prescribed conditions (smoothness, orientation preservation, physical interpretability, etc.). The outline of the article follows this progression from the single generating functional to the multiscale scheme, and provides theoretical mathematical results ensuring the well-posed character of each step as well as intuitive interpretations of them:

• The structure of Section 2 reflects this linear progression. Subsection 2.2 is devoted to the design of the parent functional $\mathcal{F}$ related to the basic registration problem and
on which the multiscale approach relies. Arguments from the theory of mechanics motivate the way this functional is built, the objects to be matched being viewed as bodies subjected to external forces. A first theoretical result (Theorem 2.3) ensures that given a pair of images $(R, T) \in BV(\Omega) \times BV(\Omega)$, this primary registration model admits at least one minimiser exhibiting fine smoothness properties. In particular, it is a bi-Lipschitz homeomorphism and self-penetration of the matter does not occur, which is mechanically mandatory. Note that the detailed proof was not provided in [14];

• Subsection 2.3 constitutes the core of our contribution, both methodologically and theoretically, since it introduces the multiscale model. Its construction is patterned after the multiscale image representation [31] insofar as the stage $k$ depends on the recovered deformations at the previous stages. However, the composition is now a substitute for the addition which was the natural operation in the context of hierarchical decompositions of images. Given $R$ (resp. $T$) the Reference (resp. Template) image and its related hierarchical decomposition $\sum_{j=0}^{k} R_j$ (resp. $\sum_{j=0}^{k} T_j$) obtained from [31], the algorithm reads as

$$ \varphi_0 = \arg \min_{\varphi} \mathcal{F}(\varphi; R_0, T_0), $$

$$ (D) \qquad \vdots $$

$$ \varphi_k = \arg \min_{\varphi} \mathcal{F}(\varphi_0 \circ \dots \circ \varphi_{k-1} \circ \varphi; \sum_{j=0}^{k} R_j, \sum_{j=0}^{k} T_j). $$

The deformation $\varphi_0$ thus maps the coarser version $R_0$ of $R$ to the coarser version $T_0$ of $T$, while iterating the procedure yields more refined deformations of the type $\varphi_0 \circ \dots \circ \varphi_i$ pairing versions of the original images encoding finer details ($\sum_{j=0}^{i} R_j$ and $\sum_{j=0}^{i} T_j$).

A first result (first part of Theorem 2.4) shows that each subproblem of the hierarchical scheme admits at least one minimiser on a suitable functional space using an induction process, while a second result (second part of Theorem 2.4) which proves to be an asymptotic result, emphasises that for $k$ large enough, the recovered deformation $\varphi_0 \circ \dots \circ \varphi_k$ constitutes a good approximation of the deformation that maps $R$ and $T$. Note that the detailed proof was not provided in [14]; The added value of the model is that it allows the control and analysis of the granularity of the obtained deformation by selecting carefully the tuning parameter $\lambda$ of the hierarchical decomposition of images.

• If the model exhibits desirable properties (well-posedness, exhaustiveness of the family of generated deformations, physical interpretation) linked in particular to its non-linear character, it falls within the non-convex and non-differentiable class of optimisation problems which is the hardest one to solve numerically. Section 3 aims to split the initial problem into subproblems encoding only a part of the numerical difficulty by means of auxiliary variables, the underlying goal being to alleviate the computational burden. These auxiliary variables are related to the variables they are supposed to simulate by $L^p$-penalisations weighted by a parameter $\gamma$ (doomed to be large. The larger parameter $\gamma$, the closer the auxiliary variable is to the
quantity it simulates in $L^p$-norm). As previously mentioned, an asymptotic result
that is a result by itself and that will be the topic of a paper on its own ([13])
can be proved. It says in substance that for a sufficiently large parameter $\gamma$, the
decoupled problem is a suitable approximation of the original problem.

• Section 4 is dedicated to numerical experiments and to their analysis. The question of evaluating the potential of the proposed model encompasses several levels of discussion: ability of the method to both discriminate between main global tendency and more localised displacements and model deformations capturing increasingly fine details (in particular, synthetic images have been created, exhibiting features of different scales), stability of the algorithm regarding the choice of parameters and with respect to noise on the data, variety of the panel of generated deformations (large deformations can occur), quantitative evaluation of the registration accuracy with several metrics supporting the theoretical asymptotic result, comparisons with related methods etc.. Such a discussion was not provided in [14].

To summarise, our contributions are of different kinds: (i) first, of a methodological nature, by providing a nested algorithm capable of representing a deformation pairing two images in a multiscale fashion; (ii) second, of a more theoretical nature, by devising several theoretical results supporting the soundness of the model and which can be interpreted intuitively; (iii) at last, of a more applied nature, with multi-factorial evaluations that sustain the theoretical results and the intended objectives.

**1.3. Prior works.** Before depicting in depth our model and for the sake of completeness, we review some prior related works and highlight the main differences with the proposed work.

Prior related works ([24, 25, 26], [23]) suggest fostering the use of this multiscale representation of images —separation of the coarse and fine scales —in the context of registration. The work [24] (and then its extensions to landmark-driven registration in a B-spline setting from the one hand ([25]), and non-rigid deformations from the other hand ([26])) focuses on spatial alignment of medical images degraded by significant levels of noise, the underlying goal being to highlight the real differences due to actual variations of the objects, while removing artificial deviation. It is achieved by mapping the truncated hierarchical representations of both images. Two main differences can be noticed compared to our model: (i) first, unlike our model, there is an independent treatment at each hierarchical level. No connection is made between the deformations $\phi_0, \dots, \phi_{k-1}$ computed at the previous iterations and the current deformation $\phi_k$; (ii) second, in [24], the final optimal deformation meant to map the two images is computed as a weighted average of the form $\frac{1}{m} \sum_{l=0}^{m-1} b_l \phi_l$ with suitable weights $b_l$'s, while we promote composition of deformations as it is the most natural and geometrically meaningful operation the space of non-parametric spatial transformations can be endowed with. Indeed, mapping a point through a first transformation and then through a second one amounts to mapping the point through the composition of these two spatial transformations. On the contrary, except for the case of small deformations where linearisation is applied, addition of spatial deformations has no geometrical meaning.

A work closer to ours can be found in [23] in the sense that analogous hierarchical expan-
sions of diffeomorphisms as composition of maps are constructed. This model can be viewed
as a sequence of Large Deformation Diffeomorphic Metric Mapping (LDDMM, [5])-based
steps fine-tuned with suitable weighting parameters. Here, the adjective multiscale applies
to the setting of these parameters, i.e., on how strong the penalisations on the similarity
measure and on the deviation from the identity mapping are. If our method and [23] share
this idea of composing deformations to refine the registration process, some differentiating
points can be highlighted in addition to the one mentioned above: among them, the point
of view we adopt to describe the framework in which the objects to be registered are viewed.
If physical assumptions and more precisely hyperelasticity arguments promoting large and
nonlinear deformations rule the design of our model, [23] is purely built on mathematical
considerations and strays to some extent to the physics of the problem.

Starting again from the observation that a deformation is a combination of local and
global deformations of different scales and locations, the work [19] proposes decomposing
an orientation-preserving deformation into different components —each one inheriting this
property of sense preservation —, based on the theory of quasiconformal mappings ([2],
[20]) which are mainly mappings of bounded distortion. Several facts motivate the ap-
proach: (i) first, a homeomorphism $f$ is $K$-quasiconformal if and only if $f$ is an $L^2$-solution
of an equation of the type $\partial f = \mu \partial f$, where $\mu$, named complex dilation or Beltrami Coef-
ficient (BC), satisfies $|\mu(z)| \le \frac{K-1}{K+1} < 1$ for almost every $z$ ([20, Theorem 4.1]) —note that
$\mu$ is a measure of non conformality: it quantifies to what extent a deformation deviates
from a conformal map and $\|\mu\|_{\infty} < 1$ implies that $f$ is sense-preserving —; (ii) second, by
shifting the focus ([20, Theorem 4.4]): given a measurable function $\mu$ in a domain $A$ with
$\|\mu\|_{\infty} < 1$, there exists a quasiconformal mapping of $A$ whose complex dilation agrees with
$\mu$ almost everywhere; (iii) third, adjusting a mapping by working with its complex dilation
is easier than handling its coordinate functions. Thus based on these elements and on
the fact that a deformation is entirely described by its associated BC, once an orientation-
mapping $f$ is extracted and its complex dilation $\mu(z) = \mu_f(z) = \left(\frac{\partial f}{\partial z}/\frac{\partial f}{\partial z}\right)$ is computed, the
authors suggest applying a wavelet transform to $\mu$, yielding a decomposition into distinct
components of different frequencies compactly supported on different sub-domains. The
multiscale components of the deformation are then recovered by converting the successive
normalised —in order to ensure that the supreme norm is strictly less than 1 —truncations
of the wavelet transform into their associated quasiconformal map, yielding a sequence of
deformations encoding finer and finer details. This is achieved by solving elliptic PDE’s
derived from Beltrami equations.

If their method and ours agree on this latter point, there are however, beyond the math-
ematical formalism, dissimilarities. Mainly, on the structuring of the algorithms: in [19],
the original orientation-preserving deformation pairing the two images is an input, and the
multiscale decomposition of the deformation is computed only from the related complex
dilation by abstraction of the different levels of details encapsulated in the images. The
procedures of registration and deformation decomposition are thus independent with each
others. On the contrary, in our method, the deformation allowing to match the two images
is the expected output, at least from a theoretical point of view since it is viewed as the
asymptotic behaviour of the intermediate deformation composition, and is computed tak-
ing into account the level of granularity of the image constituents. Our model thus sticks more to the information contained in the images and reflects more faithfully the features of the underlying deformation involved in the registration process. Also, in [19], deformation analysis can be carried out locally but requires to introduce a mask function on the wavelet coefficients, while it is more straightforward in our approach as will be seen in section 2. Finally and for the sake of completeness, we refer the reader to [3], [17], [27] and [28] for alternative approaches.

We now turn to the mathematical foundations of our physics-based multiscale registration model. We would like to emphasise that the focus of the paper is on the mathematical analysis of the proposed model including well-posedness of the original minimisation problem, asymptotic behaviour (meaning that the deformation obtained at step $k$ converges to the deformation matching the two images as $k$ increases to $+\infty$), suitable algorithm, etc. The model is restricted to the two-dimensional case. Further work will be dedicated to higher dimensions (2D, 3D, 3D+t) and to the ability of the model to unveil the hidden structure of a deformation. The study will also be enriched by a theoretical analysis of the proposed numerical algorithm, which will be the subject of a forthcoming paper ([13]).

## 2. Mathematical modelling.

### 2.1. Motivations.
If image decomposition aims to partition a given image $f$ into the sum of a structural part encoding the main geometrical features and a texture component $v$ capturing the oscillatory patterns or noise, multiscale image representation goes beyond by reckoning the different levels of details of an image. In line with this idea of hierarchising the information carried by an object —in our case, a mapping —and relying on the multiscale image representation [31], we propose quantifying the noticeable characteristics of a deformation matching two images through the behaviour of a family of functionals. We first introduce the original minimisation problem based on the parent functional $\mathcal{F}$ from which the multiscale model will be derived.

### 2.2. Hyperelastic setting for the original minimisation problem.
Let $\Omega$ be a convex bounded open subset of $\mathbb{R}^2$ of class $C^1$ therefore satisfying the cone property. This latter requirement is for technical purposes to ensure that Ball's theorems ([4]) apply. It means that there exists a finite cone $C$ such that each point $x \in \Omega$ is the vertex of a finite cone $C_x$ contained in $\Omega$ and congruent to $C$. The moving Template image is represented by $T: \bar{\Omega} \to \mathbb{R}$, while the fixed Reference image is denoted by $R: \bar{\Omega} \to \mathbb{R}$. These are assumed to belong to the functional space $BV(\Omega)$. For theoretical purposes, we assume that $T$ is such that its essential support $\text{ess supp}(T)$ is included in $\Omega' \subset \subset \Omega$, $\Omega'$ being a bounded open set of $\Omega$. The mapping $\varphi: \bar{\Omega} \to \mathbb{R}^2$ is the sought non-parametric non-rigid deformation matching the two images. The deformation gradient is $\nabla\varphi: \Omega \to M_2(\mathbb{R})$, with $M_2(\mathbb{R})$ the set of $2 \times 2$ matrices. Mechanically [9], a deformation is a smooth mapping that is orientation-preserving and injective except possibly on $\partial\Omega$ where self-contact is authorised. This translates mathematically into the condition $\det \nabla\varphi > 0$ almost everywhere. This property should be included into the deformation model prescribing the nature of the allowed deformations if one aims to get physically meaningful and sense-preserving ones.
287

288

289

290

291

292

Remark 2.1. We acknowledge the fact that the deformation should be with values in $\bar{\Omega}$ in practice. However, from a mathematical point of view, if we work with such elements we lose the structure of vector space which is essential for the theoretical analysis. Furthermore, thanks to Ball’s results [4], we show in the sequel that our model generates deformations with values in $\bar{\Omega}$.

Since the registration problem is ill-posed, in a variational setting, the sought deformation is obtained by minimising a functional $\mathcal{F}$, that we call parent functional, comprising a fidelity term quantifying how close the deformed Template is to the Reference, and an additional regularisation acting as a deformation model. Since we focus on mono-modal registration in this work, we propose using the classical sum of squared difference metric to measure alignment:

$$ Fid(\varphi) = \|T \circ \varphi - R\|_{L^2(\Omega)}^2. $$

Several stances can be adopted to depict the deformation model the objects to be matched fall within. This deformation model must be a good compromise between computational efficiency and completeness of the generated family of deformations. As suggested in [29], the geometrical deformations can be classified into three categories: (i) those inspired by physical models and more precisely, by principles of mechanics: the objects contained in the images are viewed as bodies subjected to forces; (ii) those derived from interpolation and approximation theory and at last, (iii) those stemming from models including a priori knowledge such as biomechanical models whose design is dictated by specific anatomical/-physiological laws.

Our model falls within the former category and is more particularly part of the hyperelasticity setting, good compromise between computational performance and exhaustiveness of the panel of generated deformations, since including large deformations (please refer to [9, Part A, Chapter 4] for an introduction to hyperelasticity). Hyperelasticity provides a means of modeling the stress-strain behavior of certain highly deformable materials for which linear elasticity principles are inaccurate since too simplistic/reductive. A common example of this kind of material is rubber, whose stress-strain relationship can be defined as non-linearly elastic, isotropic and incompressible. Unlike linear elasticity defined explicitly by Hooke’s law for small deformations, the hyperelasticity framework postulates the existence of a stored energy density function whose derivatives with respect to the deformation in a given direction give the state of stress within the material in this same direction. Coming back to our model, the objects contained in the images are assumed to be isotropic (exhibiting the same mechanical properties in every direction), homogeneous (showing the same behaviour everywhere inside the material), and so hyperelastic (allowing large changes on shape while keeping a mechanical elastic behaviour) materials, and more precisely as Ogden ones (please refer to [9, Part B, Chapter 7]). Note that the Ogden material model is often used to describe the non-linear stress-strain behaviour of complex materials such as rubbers or biological tissues. This perspective drives the design of the regularisation on the deformations, which is thus based on the stored energy function of
an Ogden material, prescribing then a physically-meaningful nature. These elements are mathematically formalised next.

In 2D, for such a material [9, Part B, Chapter 7], the general expression of the stored energy with $F \in M_2(\mathbb{R})$ is :

$$W_O(F) = \sum_{i=1}^{K} a_i \|F\|^{\gamma_i} + \Gamma(\det F),$$

where $\forall i \in \{1, \dots, K\}$, $a_i > 0$, $\gamma_i \ge 0$ are material parameters, $\|\cdot\| = \|\cdot\|_F$ being the Frobe-

nius norm ($\|F\| = \sqrt{\det F^T F} = \sqrt{f_{11}^2 + f_{12}^2 + f_{21}^2 + f_{22}^2}$ if one sets $F = \begin{pmatrix} f_{11} & f_{12} \\ f_{21} & f_{22} \end{pmatrix}$), and

$\Gamma:]0, +\infty[ \to \mathbb{R}$ being any convex function satisfying $\lim_{\delta \to 0^+} \Gamma(\delta) = \lim_{\delta \to +\infty} \Gamma(\delta) = +\infty$.

Here $F$ denotes the deformation gradient $\nabla\varphi$. The first terms influence the changes in length, while the last one restricts the changes in area and ensures orientation preservation by preventing the Jacobian determinant from becoming negative. In this work, we introduce the following particular energy —$Op$ stands for Ogden particular —:

$$W_{Op}(F) = W_{Op}(F, \det F)$$

$$= \begin{cases} a_1 \|F\|^4 + a_2(\det F - 1)^2 + \frac{a_3}{(\det F)^{10}} - 4a_1 - a_3 & \text{if } \det F > 0 \\ +\infty & \text{otherwise} \end{cases},$$

which fulfils the previous assumptions and exhibits fine theoretical properties useful for the mathematical analysis conducted in the sequel. In particular, $W_{Op}$ is polyconvex since $W_{Op} := W_{Op}(F, \delta)$ is convex. Moreover, the choice of the power 4 in $\|F\|^4 - 4$ being strictly greater than 2, the dimension of the ambient domain —combined with the regularisation $R(F)$ below and the constraint $\det F > 0$ a.e. allows to recover deformations that are homeomorphisms as will be seen later). The first term controls the smoothness of the deformation, the second one restricts changes in area since promoting Jacobian determinant close to 1, while the third one prevents singularities and large contractions by penalising small values of the determinant. The last two constants are added to comply with the energy property $W_{Op}(I) = 0$, $I$ denoting the identity matrix, Jacobian of the identity mapping. We propose complementing this regularisation by the following term :

$$R(F) = 1_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(F) + 1_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}}(F^{-1}),$$

with $\alpha \ge 1$ and $\beta \ge 1$, $1_A$ being the convex characteristic function of a convex set $A$. This ensures that the obtained deformations are bi-Lipschitz homeomorphisms and subsequently, $T$ being an element of the space $BV(\Omega)$, that $T \circ \varphi$ remains in $BV(\Omega)$ according to [1, Theorem 3.16].

**Remark 2.2.** This additional constraint implicitly gives an upper and lower bound on the Jacobian determinant, controlling thus the amount of contraction and dilation allowed while preserving topology.
The proposed registration model in a variational setting therefore reads:

$$
\begin{equation}
\begin{aligned}
& \inf_{\varphi \in W} \left\{ \begin{aligned}[t]
& \mathcal{F}(\varphi) = \mathcal{F}(\varphi, T, R) = \frac{\lambda}{2} \mathrm{Fid}(\varphi) + \int_{\Omega} W_{Op}(\nabla \varphi) \, dx + R(\nabla \varphi), \\
& \mathcal{P} \qquad = \frac{\lambda}{2} \|T \circ \varphi - R\|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \varphi, \det \nabla \varphi) \, dx + 1_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \varphi) \\
& \qquad + 1_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}} ((\nabla \varphi)^{-1}) \end{aligned} \right\},
\end{aligned}
\tag{3.1}
\end{equation}
$$

with $\mathcal{W} = \{\psi \in \mathrm{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2) \mid \|\nabla\psi\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha, \|(\nabla\psi)^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta, \det \nabla\psi > 0 \text{ a.e. in } \Omega\}$, and $\lambda > 0$ a weighting parameter balancing the influence of the fidelity term with respect to the regularisation one. The first theoretical result claims that problem $(\mathcal{P})$ admits at least one minimiser. In particular, this result guarantees that the recovered deformations exhibit smoothness properties and that they are mechanically admissible with no self-intersection of matter.

Theorem 2.3. Problem $(\mathcal{P})$ admits at least one minimiser in $\mathcal{W}$.

*Proof.* The proof follows the arguments of the classical direct method of the calculus of variations. We first derive a coercivity inequality. Using the fact that $(a-b)^2 \ge \frac{1}{2}a^2 - b^2$, one has

$$
\begin{equation}
\begin{aligned}
\mathcal{F}(\varphi) \geq & a_1 \|\nabla\varphi\|_{L^4(\Omega, M_2(\mathbb{R}))}^4 + \frac{a_2}{2} \|\det \nabla\varphi\|_{L^2(\Omega)}^2 - a_2 \mathrm{meas}(\Omega) + \| \frac{a_3}{(\det \nabla\varphi)^{10}} \|_{L^1(\Omega)} \\
& - 4a_1 \mathrm{meas}(\Omega) - a_3 \mathrm{meas}(\Omega).
\end{aligned}
\tag{3.2}
\end{equation}
$$

The quantity $\mathcal{F}(\varphi)$ is thus bounded below by $-(4a_1 + a_2 + a_3)\mathrm{meas}(\Omega)$ and as for $\varphi = \mathrm{Id}$—and suitable $\alpha$ and $\beta$—, $\mathcal{F}(\varphi) = \frac{\lambda}{2}\|T - R\|_{L^2(\Omega)}^2$ is finite (due to the embedding $BV(\Omega) \subset L^2(\Omega)$ in the two-dimensional case), the infimum is finite.

Let then $(\varphi_k)_k \in \mathcal{W}$ be a minimising sequence—we may omit the index $k$ in the following when dealing with a sequence indexed by $k$, i.e., $\lim_{k \to +\infty} \mathcal{F}(\varphi_k) = \inf_{\Psi \in \mathcal{W}} \mathcal{F}(\Psi)$. Hence there exists $K \in \mathbb{N}$ such that $\forall k \in \mathbb{N}$, $(k \ge K \Rightarrow \mathcal{F}(\varphi_k) \le \inf_{\Psi \in \mathcal{W}} \mathcal{F}(\Psi) + 1)$. From now on, we assume that $k \ge K$. According to the coercivity inequality, one gets:

* $(\varphi_k)$ is uniformly bounded according to $k$ in $W^{1,4}(\Omega, \mathbb{R}^2)$, using the generalised Poincaré inequality ([15, pp. 106-107]) and the fact that $\varphi_k = \mathrm{Id}$ on $\partial\Omega$;

* $(\nabla\varphi_k)$ is uniformly bounded according to $k$ in $L^\infty(\Omega, M_2(\mathbb{R}))$;

* $(\nabla\varphi_k)^{-1}$ is uniformly bounded according to $k$ in $L^\infty(\Omega, M_2(\mathbb{R}))$;

* $(\det \nabla\varphi_k)$ is uniformly bounded according to $k$ in $L^2(\Omega)$.

Thus there exist a subsequence —still denoted by $(\varphi_k)$ —and $\bar{\varphi} \in W^{1,4}(\Omega, \mathbb{R}^2)$ such that

$$
\varphi_k \xrightarrow{k \to +\infty} \bar{\varphi} \text{ in } W^{1,4}(\Omega, \mathbb{R}^2).
$$

Moreover, there exist a subsequence (common with the previous one, which is always possible) —still denoted by $(\det \nabla \varphi_k)$ —and $\delta \in L^2(\Omega)$ such that

$$
\det \nabla \varphi_k \xrightarrow{k \to +\infty} \delta \text{ in } L^2(\Omega).
$$
By applying [10, Theorem 8.20], we deduce that $\delta = \det \nabla \bar{\varphi}$ and $\det \nabla \varphi_k \underset{k \to +\infty}{\to} \det \nabla \bar{\varphi}$ in $L^2(\Omega)$.

Now, recall that according to Gagliardo-Nirenberg inequalities ([6, p.195, Example 3], [7]), whenever $1 \le q \le p \le +\infty$ and $r > N$ (in the general case where $\Omega \subset \mathbb{R}^N$ is a bounded open set with smooth boundary),

$$ \|u\|_{L^p} \le C \|u\|_{L^q}^{1-a} \|u\|_{W^{1,r}}^a, \forall u \in W^{1,r}(\Omega), $$

with $a = \frac{\frac{1}{q} - \frac{1}{p}}{\frac{1}{q} + \frac{1}{N} - \frac{1}{r}}$.

In our case, $(\varphi_k)$ being uniformly bounded in $W^{1,4}(\Omega, \mathbb{R}^2)$ and taking $\begin{cases} p = +\infty \\ q = 4 \\ r = 4 \end{cases}$ (yielding $a = \frac{1}{2}$),

$$ \|\varphi_k\|_{L^\infty(\Omega, \mathbb{R}^2)} \le C \|\varphi_k\|_{L^4(\Omega, \mathbb{R}^2)}^{\frac{1}{2}} \|\varphi_k\|_{W^{1,4}(\Omega, \mathbb{R}^2)}^{\frac{1}{2}}, $$

showing in the end that $(\varphi_k)$ is uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$. Thus there exist a sub-sequence (common with the previous one) —still denoted by $(\varphi_k)$ —and $\bar{\varphi} \in W^{1,\infty}(\Omega, \mathbb{R}^2)$ such that

$$ \varphi_k \underset{k \to +\infty}{\stackrel{*}{\rightharpoonup}} \bar{\varphi} \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2). $$

In particular, one has $\varphi_k \underset{k \to +\infty}{\to} \bar{\varphi}$ in $L^4(\Omega, \mathbb{R}^2)$ and $\varphi_k \underset{k \to +\infty}{\stackrel{*}{\rightharpoonup}} \bar{\varphi}$ in $L^\infty(\Omega, \mathbb{R}^2)$, so that owing to the property of uniqueness of the weak limit in $L^4(\Omega, \mathbb{R}^2)$, $\bar{\varphi} = \bar{\varphi}$ in $L^4(\Omega, \mathbb{R}^2)$. By definition of the functional space $W^{1,4}(\Omega, \mathbb{R}^2)$, $\forall i \in \{1, 2\}$, $\exists (g_{ij})_{j=1,2} \in L^4(\Omega)$ such that $\forall \Psi \in C_0^\infty(\Omega)$,

$$ \int_{\Omega} \bar{\varphi}_i \frac{\partial \Psi}{\partial x_j} dx = \int_{\Omega} \bar{\varphi}_i \frac{\partial \Psi}{\partial x_j} dx = - \int_{\Omega} g_{ij} \Psi dx, $$

leading to $\bar{\varphi} = \bar{\varphi} \in W^{1,\infty}(\Omega, \mathbb{R}^2)$. At last, by continuity of the trace operator ([6, Theorem III.9], [7]), we get that $\bar{\varphi} \in \text{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2)$.

The stored energy function $W_{Op}$ is continuous and convex. If $\psi_n \underset{n \to +\infty}{\to} \bar{\psi}$ in $W^{1,4}(\Omega, \mathbb{R}^2)$, thus $\nabla \psi_n \underset{n \to +\infty}{\to} \nabla \bar{\psi}$ in $L^4(\Omega, M_2(\mathbb{R}))$ and one can extract a subsequence still denoted by $(\nabla \psi_n)$ such that $\nabla \psi_n \underset{n \to +\infty}{\to} \nabla \bar{\psi}$ almost everywhere in $\Omega$. Similarly, if $\kappa_n \underset{n \to +\infty}{\to} \bar{\kappa}$ in $L^2(\Omega)$, then one can extract a subsequence (common subsequence) still denoted by $(\kappa_n)$ such that $\kappa_n \underset{n \to +\infty}{\to} \kappa$ almost everywhere in $\Omega$. Then, by continuity of $W_{Op}$, one gets that $W_{Op}(\nabla \psi_n, \kappa_n) \underset{n \to +\infty}{\to} W_{Op}(\nabla \bar{\psi}, \bar{\kappa})$ almost everywhere in $\Omega$. Applying Fatou's lemma then yields

$$ \int_{\Omega} W_{Op}(\nabla \bar{\psi}, \bar{\kappa}), dx \le \liminf_{n \to +\infty} \int_{\Omega} W_{Op}(\nabla \psi_n, \kappa_n) dx. $$
As $W_{Op}$ is convex, so is $\int_{\Omega} W_{Op}(\xi, \kappa) dx$ and invoking [6, Corollaire III.8], [7] leads to:

$$ \int_{\Omega} W_{Op}(\nabla \bar{\varphi}, \det \nabla \bar{\varphi}) dx \leq \liminf_{k \to +\infty} \int_{\Omega} W_{Op}(\nabla \varphi_k, \det \nabla \varphi_k) dx < +\infty. $$

Since $W_{Op}(\nabla \bar{\varphi}(x), \det \nabla \bar{\varphi}(x)) = +\infty$ when $\det \nabla \bar{\varphi}(x) \le 0$, the set on which it occurs is necessarily of null measure otherwise we would have $F(\bar{\varphi}) = +\infty$. So $\det \nabla \bar{\varphi} > 0$ almost everywhere in $\Omega$. Besides, for all $q > 2$ and all $k \ge K$,

$$ \int_{\Omega} \|(\nabla \varphi_k)^{-1}\|_F^q \det \nabla \varphi_k dx \le C, $$

$C = C(\alpha, \beta, q, \Omega) > 0$ being a constant depending only on $\alpha, \beta$ and $\Omega$.

The assumptions of Ball's theorems ([4, Theorems 1 and 2]) thus hold yielding that $\varphi_k$ is a homeomorphism of $\bar{\Omega}$ onto $\bar{\Omega}$ and $\varphi_k^{-1} \in W^{1,q}(\Omega, \mathbb{R}^2)$. The mapping $\varphi_k$ is one-to-one almost everywhere, maps measurable sets in $\bar{\Omega}$ to measurable sets in $\bar{\Omega}$, and the change of variables formula

$$ \int_A f(\varphi(x)) \det \nabla \varphi(x) dx = \int_{\varphi(A)} f(v) dv $$

holds for any measurable $A \subset \bar{\Omega}$ and any measurable function $f : \mathbb{R}^2 \to \mathbb{R}$, provided only that one of the integrals exists. The matrix of weak derivatives of $\varphi_k^{-1}$ is given by $\nabla(\varphi_k)^{-1} = (\nabla \varphi_k)^{-1}(\varphi_k^{-1})$ almost everywhere in $\Omega$.

Let $\mathcal{N}_k \subset \Omega$ be such that $\mathrm{meas}(\mathcal{N}_k) = 0$ and for all $x \in \Omega \setminus \mathcal{N}_k$,

$$ |(\nabla \varphi_k)^{-1}(x)| \le \|(\nabla \varphi_k)^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}. $$

Let now $\mathcal{N}'_k$ be such that $\mathcal{N}'_k = \varphi_k(\mathcal{N}_k)$. Then $\mathrm{meas}(\mathcal{N}'_k) = 0$ since $\varphi_k$ is a Lipschitz map (and thus for every measurable set $E$, $\mathrm{meas}(\varphi_k(E)) \le C'$meas$(E)$, $C'$ being a constant depending only on the dimension and on the Lipschitz constant of $\varphi_k$ itself uniformly bounded with respect to $k$) and for every $y \notin \mathcal{N}'_k$, $|(\nabla \varphi_k)^{-1}(\varphi_k^{-1}(y))| \le \|(\nabla \varphi_k)^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}$. Indeed, if $y \notin \mathcal{N}'_k$, then $y = \varphi_k(x)$ with $x \notin \mathcal{N}_k$, resulting in $(\nabla \varphi_k)^{-1}(\varphi_k^{-1}(y)) = (\nabla \varphi_k)^{-1}(x)$.

Consequently,

$$ \|\nabla(\varphi_k^{-1})\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} = \|(\nabla \varphi_k)^{-1}(\varphi_k^{-1})\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \|(\nabla \varphi_k)^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}. $$

Invoking again the generalised Poincaré inequality and Gagliardo-Nirenberg inequalities allows to conclude that $\varphi_k^{-1}$ is uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$.

Applying the same reasoning as the one we did previously, that is, denoting by $\mathcal{N}_k \subset \Omega$, set of null measure such that $\forall x \in \Omega \setminus \mathcal{N}_k$,

$$ |(\nabla \varphi_k)^{-1}(\varphi_k^{-1})(x)| \le \|(\nabla \varphi_k)^{-1}(\varphi_k^{-1})\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}. $$

Let $\mathcal{N}'_k$ be such that $\mathcal{N}'_k = \varphi_k^{-1}(\mathcal{N}_k)$. Then $\mathrm{meas}(\mathcal{N}'_k) = 0$ since $\varphi_k^{-1}$ is a Lipschitz map. For every $y \notin \mathcal{N}'_k$, $|(\nabla \varphi_k)^{-1}(\varphi_k^{-1} \circ \varphi_k(y))| \le \|(\nabla \varphi_k)^{-1}(\varphi_k^{-1})\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}$. Consequently,

$$ \| (\nabla \varphi_k)^{-1} \|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \| (\nabla \varphi_k)^{-1} (\varphi_k^{-1}) \|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} = \| \nabla (\varphi_k^{-1}) \|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}. $$
422 By gathering the two previous results, it follows that

$$ \| (\nabla \varphi_k)^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \| \nabla (\varphi_k^{-1}) \|_{L^\infty(\Omega, M_2(\mathbb{R}))}. \quad 423 $$

425 By the weak-* lower semi-continuity of $\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))}$, we deduce that $\|\nabla \bar{\varphi}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le$

$$ \liminf_{k \to +\infty} \|\nabla \varphi_k\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha \text{ and} \quad 426 $$

$$ 427 \quad 1\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\} (\nabla \bar{\varphi}) \le \liminf_{k \to +\infty} 1\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\} (\nabla \varphi_k). \text{ Also, for all } q \in [2, 11], $$

$$ 428 \quad \int_{\Omega} \|(\nabla \bar{\varphi})^{-1}\|_{F}^{q} \det \nabla \bar{\varphi} \, dx = \int_{\Omega} \frac{1}{(\det \nabla \bar{\varphi})^{q}} \|\nabla \bar{\varphi}\|_{F}^{q} \det \nabla \bar{\varphi} \, dx, $$

$$ 429 \quad = \int_{\Omega} \|\nabla \bar{\varphi}\|_{F}^{q} (\det \nabla \bar{\varphi})^{1-q} \, dx, $$

$$ 430 \quad \le \alpha^q \frac{1}{\det \nabla \bar{\varphi}} \|\nabla \bar{\varphi}\|_{L^{q-1}(\Omega)}^{q-1} < +\infty, \quad 431 $$

432 since $\det \nabla \bar{\varphi} > 0$ almost everywhere and $\frac{1}{\det \nabla \bar{\varphi}} \in L^{10}(\Omega)$.

433 Ball's theorems ([4, Theorems 1 and 2]) allow to conclude that $\bar{\varphi}$ is a homeomorphism from 434 $\bar{\Omega}$ to $\bar{\Omega}$, $\bar{\varphi}^{-1} \in W^{1,q}(\Omega, \mathbb{R}^2) \subset L^\infty(\Omega, \mathbb{R}^2)$ (continuous embedding) and an upper bound of 435 $\|\bar{\varphi}\|_{L^\infty(\Omega, \mathbb{R}^2)}$ with respect to $\|\bar{\varphi}\|_{W^{1,q}(\Omega, \mathbb{R}^2)}$ is easily obtained.

436 The sequence $(\varphi_k^{-1})$ being uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$, there exists a subsequence 437 still denoted by $\varphi_k^{-1}$ and $\bar{u} \in W^{1,\infty}(\Omega, \mathbb{R}^2)$ such that

$$ 438 \quad \varphi_k^{-1} \underset{k \to +\infty}{\stackrel{*}{\longrightarrow}} \bar{u} \quad 439 $$

440 in $W^{1,\infty}(\Omega, \mathbb{R}^2)$. Let us now prove that $\bar{u} = \bar{\varphi}^{-1}$. Due to Rellich-Kondrachov theorem, the 441 compact injection $W^{1,\infty}(\Omega, \mathbb{R}^2) \subset C^0(\bar{\Omega}, \mathbb{R}^2)$ holds so that $(\varphi_k)$ uniformly converges to $\bar{\varphi}$ 442 on $\bar{\Omega}$ while $(\varphi_k^{-1})$ uniformly converges to $\bar{u}$ on $\bar{\Omega}$. Also,

$$ 443 \quad \|\varphi_k^{-1} \circ \bar{\varphi} - \varphi_k^{-1} \circ \varphi_k\|_{C^0(\bar{\Omega}, \mathbb{R}^2)} \le \beta \|\bar{\varphi} - \varphi_k\|_{C^0(\bar{\Omega}, \mathbb{R}^2)} \xrightarrow{k \to +\infty} 0, \quad 444 $$

445 leading to $\varphi_k^{-1} \circ \bar{\varphi} \xrightarrow{k \to +\infty} \text{Id}$ pointwise everywhere on $\bar{\Omega}$. But as $(\varphi_k^{-1})$ uniformly converges 446 to $\bar{u}$ on $\bar{\Omega}$, for all $x \in \bar{\Omega}$,

$$ 447 \quad \varphi_k^{-1} \circ \bar{\varphi}(x) \xrightarrow{k \to +\infty} \bar{u} \circ \bar{\varphi}(x). \quad 448 $$

449 By uniqueness of the pointwise limit, $\bar{u} \circ \bar{\varphi} = \text{Id}$ on $\bar{\Omega}$, resulting in $\bar{u} = \bar{\varphi}^{-1}$ everywhere 450 on $\bar{\Omega}$ and $\bar{\varphi}^{-1} \in W^{1,\infty}(\Omega, \mathbb{R}^2)$. The mapping $\bar{\varphi}$ is thus a bi-Lipschitz homeomorphism.

451 Invoking again the weak-* lower semi-continuity of $\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))}$ and arguing with the 452 same arguments as before yields

$$ 453 \quad \|\nabla (\bar{\varphi}^{-1})\|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \| (\nabla \bar{\varphi})^{-1} \circ \bar{\varphi}^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \| (\nabla \bar{\varphi})^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))} $$

$$ 454 \quad \le \liminf_{k \to +\infty} \|\nabla (\varphi_k^{-1})\|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \liminf_{k \to +\infty} \|(\nabla \varphi_k)^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta, \quad 455 $$
so that $1_{\{\|\cdot\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \beta\}}((\nabla \bar{\varphi})^{-1}) \le \liminf_{k \to +\infty} 1_{\{\|\cdot\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \beta\}}((\nabla \varphi_k)^{-1})$.

Since $T \in BV(\Omega)$ and all $\varphi_k$ and $\bar{\varphi}$ are bi-Lipschitz homeomorphisms, we get that $T \circ \varphi_k \in BV(\Omega) \subset L^2(\Omega)$ for all $k \in \mathbb{N}$ and $T \circ \bar{\varphi} \in BV(\Omega)$ from [1, Theorem 3.16]. We first prove that $\varphi_k \circ \bar{\varphi}^{-1} \xrightarrow{k \to +\infty} \text{Id}$ in $C^{0,\alpha}(\bar{\Omega}, \mathbb{R}^2)$ with $\alpha < 1$. Recall that (see [10, Definition 12.5]) with $0 < \alpha \le 1$, $C^{0,\alpha}(\bar{\Omega}, \mathbb{R}^2)$ is the set of functions $u \in C^0(\bar{\Omega}, \mathbb{R}^2)$ such that

$$ [u]_{\alpha, \bar{\Omega}} := \sup_{(x,y) \in \bar{\Omega} \times \bar{\Omega} \atop x \neq y} \frac{|u(x) - u(y)|}{|x - y|^{\alpha}} < +\infty. $$

It is equipped with the norm

$$ \|u\|_{C^{0,\alpha}(\bar{\Omega}, \mathbb{R}^2)} := \|u\|_{C^0(\bar{\Omega}, \mathbb{R}^2)} + [u]_{\alpha, \bar{\Omega}}. $$

Additionally, $W^{1,\infty}(\Omega, \mathbb{R}^2) \subset C^{0,\lambda}(\bar{\Omega}, \mathbb{R}^2)$ ([10, Sobolev embedding theorem, Theorem 12.11]) for every $\lambda \in [0, 1]$ and the embedding is compact for every $0 \le \lambda < 1$ ([10, Rellich-Kondrachov theorem, Theorem 12.12]). Straightforward computations thus give

$$ \begin{aligned} & \|\varphi_k \circ \bar{\varphi}^{-1} - \bar{\varphi} \circ \bar{\varphi}^{-1}\|_{C^{0,\alpha}(\bar{\Omega}, \mathbb{R}^2)} \le \sup_{x \in \bar{\Omega}} |\varphi_k(x) - \bar{\varphi}(x)| \\ & \qquad + \sup_{(x,y) \in \bar{\Omega} \times \bar{\Omega} \atop x \neq y} \frac{|\varphi_k \circ \bar{\varphi}^{-1}(x) - \varphi_k \circ \bar{\varphi}^{-1}(y) - \bar{\varphi} \circ \bar{\varphi}^{-1}(x) + \bar{\varphi} \circ \bar{\varphi}^{-1}(y)|}{|\bar{\varphi}^{-1}(x) - \bar{\varphi}^{-1}(y)|^{\alpha}} \frac{|\bar{\varphi}^{-1}(x) - \bar{\varphi}^{-1}(y)|^{\alpha}}{|x - y|^{\alpha}}, \\ & \le \sup_{x \in \bar{\Omega}} |\varphi_k(x) - \bar{\varphi}(x)| + \sup_{(x,y) \in \bar{\Omega} \times \bar{\Omega} \atop x \neq y} \frac{|\varphi_k(x) - \varphi_k(y) - \bar{\varphi}(x) + \bar{\varphi}(y)|}{|x - y|^{\alpha}} \sup_{(x,y) \in \bar{\Omega} \times \bar{\Omega} \atop x \neq y} \frac{|\bar{\varphi}^{-1}(x) - \bar{\varphi}^{-1}(y)|^{\alpha}}{|x - y|^{\alpha}}, \\ & \le (1 + \|\bar{\varphi}^{-1}\|_{C^{0,1}(\bar{\Omega}, \mathbb{R}^2)}^{\alpha}) \|\varphi_k - \bar{\varphi}\|_{C^{0,\alpha}(\bar{\Omega}, \mathbb{R}^2)} \xrightarrow{k \to +\infty} 0. \end{aligned} $$

The conclusion is immediate. It remains to prove that $T \circ \varphi_k \xrightarrow{k \to +\infty} T \circ \bar{\varphi}$ in $L^2(\Omega)$.

Let $\varepsilon > 0$ be fixed. Let $(T_n)_{n \in \mathbb{N}}$ be a sequence of functions of $D(\Omega)$ such that $T_n \xrightarrow{n \to +\infty} T$ in $L^2(\Omega)$. Let $N(\varepsilon) \in \mathbb{N}$ be such that $\forall n \in \mathbb{N}$, $(n \ge N(\varepsilon) \Rightarrow \|T_n - T\|_{L^2(\Omega)}^2 \le \frac{\varepsilon}{3})$. Since $\varphi_k \circ \bar{\varphi}^{-1} \xrightarrow{k \to +\infty} \text{Id}$ in $C^{0,\alpha}(\bar{\Omega}, \mathbb{R}^2)$, it strongly converges to $\text{Id}$ in $L^2(\Omega, \mathbb{R}^2)$, so that there exists $K = K(\varepsilon) \in \mathbb{N}$ such that $\forall k \in \mathbb{N}$, $(k \ge K(\varepsilon) \Rightarrow \|\varphi_k \circ \bar{\varphi}^{-1} - \text{Id}\|_{L^2(\Omega, \mathbb{R}^2)}^2 \le \frac{\varepsilon}{3L_{N(\varepsilon)}^2}$, with $L_{N(\varepsilon)}$ the Lipschitz constant of $T_{N(\varepsilon)}$. Here the constant $C = C(\alpha, \beta, \Omega) > 0$ may change line to line. According to Ball's theorems ([4, Theorems 1 and 2]), the following
change of variable formula holds:

$$
\begin{align*}
\int_{\Omega} |T \circ \varphi_k - T \circ \bar{\varphi}|^2 dx &= \int_{\Omega} |T \circ \varphi_k \circ \bar{\varphi}^{-1} - T|^2 \det \nabla (\bar{\varphi}^{-1}) dx, \\
&\le \frac{1}{2} \|\nabla(\varphi)^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))}^2 \int_{\Omega} |T \circ \varphi_k \circ \bar{\varphi}^{-1} - T|^2 dx, \\
&\le C \left( \int_{\Omega} |T \circ \varphi_k \circ \bar{\varphi}^{-1} - T_{N(\varepsilon)} \circ \varphi_k \circ \bar{\varphi}^{-1}|^2 dx \right. \\
&\qquad \left. + \int_{\Omega} |T_{N(\varepsilon)} \circ \varphi_k \circ \bar{\varphi}^{-1} - T_{N(\varepsilon)}|^2 dx + \int_{\Omega} |T_{N(\varepsilon)} - T|^2 dx \right), \\
&\le C \left( \|T - T_{N(\varepsilon)}\|_{L^2(\Omega)}^2 + L_{N(\varepsilon)}^2 \|\varphi_k \circ \bar{\varphi}^{-1} - \mathrm{Id}\|_{L^2(\Omega, \mathbb{R}^2)}^2 + \frac{\varepsilon}{3} \right), \\
&\le C \varepsilon.
\end{align*}
$$

We thus have proved that $\forall k \in \mathbb{N}$, $(k \ge K(\varepsilon) \Rightarrow \|T \circ \varphi_k - T \circ \bar{\varphi}\|_{L^2(\Omega)}^2) \le C\varepsilon$. By gathering all the results, we get

$$
\mathcal{F}(\bar{\varphi}) \le \liminf_{k \to +\infty} \mathcal{F}(\varphi_k) = \inf_{\varphi \in \mathcal{W}} \mathcal{F}(\varphi) < +\infty,
$$

with $\bar{\varphi} \in \mathcal{W}$, which completes the proof.

With the basic framework in place, we now see how the multiscale model is constructed
from this parent functional.

2.3. Towards a multiscale representation of the deformation . Equipped with this
original minimisation problem, we now derive a multiscale representation of the deforma-
tion, relying on the hierarchical decomposition of both the Reference and the Template
into the sum of scale-varying components [31]. Let then (T_j)_j ∈ BV(Ω) ⊂ L^2(Ω) and
(R_j)_j ∈ BV(Ω) ⊂ L^2(Ω) be the sequence of scale-varying structural features of respec-
tively T and R computed from the following problems —S standing for either R or T
below —:

$$
\left\{
\begin{array}{ll}
(S_0, v_0) & = \displaystyle \operatorname*{arg\,min}_{(u,v) \in BV(\Omega) \times L^2(\Omega) \mid S=u+v} \{\lambda_0 \|v\|_2^2 + TV(u)\}, \\
(S_{j+1}, v_{j+1}) & = \displaystyle \operatorname*{arg\,min}_{(u,v) \in BV(\Omega) \times L^2(\Omega) \mid v_j=u+v} \{2^{j+1} \lambda_0 \|v\|_2^2 + TV(u)\}, \quad j=1, \dots,
\end{array}
\right.
$$

$\lambda_0$ being an initial scale parameter provided by the user. It is assumed that $T$ and $R$ have
similar scale structures and that each level of the following hierarchical decomposition of
$T$, $\left(\sum_{j=0}^{k} T_{j}\right)$, can be matched to the corresponding level of hierarchical decomposition of
$R$, $\left(\sum_{j=0}^{k} R_{j}\right)$. The related hierarchical expansion of the deformation, starting from main
structural deformations to more localised ones, and based on the composition operator —a
more natural and physically meaningful operator than addition —is derived and reads as
follows, formulation in which one recognises the parent functional :

$$
512 \quad (\mathcal{P}_0) \qquad \varphi_0 = \arg \min_{\varphi \in \mathcal{W}} \{\mathcal{F}(\varphi, T_0, R_0)\},
$$

$$
\begin{array}{ll}
513 & (\mathcal{P}_k) \qquad \varphi_k = \arg \min_{\varphi \in \mathcal{X}_k} \{ \mathcal{F}(\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \varphi, \sum_{j=0}^k T_j, \sum_{j=0}^k R_j) \}, \\
514 & \text{with} \\
515 & \phantom{\varphi_k =} \\
516 & \mathcal{F}(\varphi, T, R) = \frac{\lambda}{2} \|T \circ \varphi - R\|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \varphi, \det(\nabla \varphi)) \, dx \\
517 & \phantom{\mathcal{F}(\varphi, T, R) =} + \mathbb{1}\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\} (\nabla \varphi) + \mathbb{1}\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\} ((\nabla \varphi)^{-1}), \\
518 & \text{and} \\
519 & \mathcal{X}_k = \left\{ \psi \mid \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi \in \text{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2), \det \nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi) > 0 \text{ a.e.,} \right. \\
520 & \phantom{\mathcal{X}_k = \left\{ \right.} \left. \|(\nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi))\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha, \|(\nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi))^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta \right\}.
\end{array}
$$

$$
\begin{array}{ll}
516 & \mathcal{F}(\varphi, T, R) = \frac{\lambda}{2} \|T \circ \varphi - R\|_{L^2(\Omega)}^2 + \displaystyle\int_{\Omega} W_{Op}(\nabla \varphi, \det(\nabla \varphi)) \, dx \\
517 & \phantom{\mathcal{F}(\varphi, T, R) =} + \mathbb{1}\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\} (\nabla \varphi) + \mathbb{1}\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\} ((\nabla \varphi)^{-1}), \\
518 & \text{and} \\
519 & \mathcal{X}_k = \left\{ \psi \mid \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi \in \text{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2), \det \nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi) > 0 \text{ a.e.,} \right. \\
520 & \phantom{\mathcal{X}_k = \left\{ \right.} \left. \|(\nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi))\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha, \|(\nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi))^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta \right\}.
\end{array}
$$

$$
\begin{array}{ll}
520 & \mathcal{X}_k = \left\{ \psi \mid \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi \in \mathrm{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2), \det \nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi) > 0 \text{ a.e.,} \right. \\
521 & \phantom{\mathcal{X}_k = \left\{ \right.} \left. \|(\nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi))\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha, \|(\nabla (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \psi))^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta \right\}. \\
522 & \text{and} \\
523 & \text{The next theorem contains two results: a first one that ensures that at each step } k, \text{ the minimisation problem admits at least one solution and that this solution exhibits fine properties (smoothness, topology-preserving feature, etc.). The second one is an asymptotic result: it highlights the fact that the recovered deformation } \varphi_0 \circ \dots \circ \varphi_k \text{ constitutes a good approximation of the deformation that would map } R \text{ and } T. \text{ The added value, in comparison to a standard registration algorithm, is that our proposed algorithm yields a good multiscale approximation of this original deformation.}
\end{array}
$$

The next theorem contains two results: a first one that ensures that at each step $k$, the minimisation problem admits at least one solution and that this solution exhibits fine properties (smoothness, topology-preserving feature, etc.). The second one is an asymptotic result: it highlights the fact that the recovered deformation $\varphi_0 \circ \dots \circ \varphi_k$ constitutes a good approximation of the deformation that would map $R$ and $T$. The added value, in comparison to a standard registration algorithm, is that our proposed algorithm yields a good multiscale approximation of this original deformation.

Theorem 2.4. Problem $(\mathcal{P}_k)$ admits at least one minimiser. Additionally, $\varphi_k$ denoting a minimiser of $(\mathcal{P}_k)$ and setting $\phi_k := \varphi_0 \circ \dots \circ \varphi_k$, one has

$$
\begin{array}{ll}
532 & \lim_{k \to +\infty} \mathcal{F}(\phi_k, \sum_{j=0}^{k} T_j, \sum_{j=0}^{k} R_j) = \mathcal{F}(\bar{\phi}, T, R) := \mathcal{F}(\bar{\phi}) = \inf_{\varphi \in \mathcal{W}} \mathcal{F}(\varphi), \\
533 & \text{with } \phi_k \underset{k \to +\infty}{\stackrel{*}{\rightharpoonup}} \bar{\phi} \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2).
\end{array}
$$

$$
534 \quad \text{with } \phi_k \underset{k \to +\infty}{\stackrel{*}{\rightharpoonup}} \bar{\phi} \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2).
$$

Before giving the proof, we recall a preliminary result.

Proposition 2.5. Taken from [1, Proposition 2.13, p.46]

Let $\Omega \subset \mathbb{R}^N$ be a bounded, convex, open set, and $u : \Omega \to \mathbb{R}$. Then $u \in W^{1,\infty}(\Omega)$ if and only if $\mathrm{Lip}(u, \Omega) < +\infty$ and $\|\nabla u\|_{\mathrm{L}^\infty(\Omega)} = \mathrm{Lip}(u, \Omega)$, with

$$
539 & \mathrm{Lip}(u, \Omega) = \sup \left\{ \frac{|u(x) - u(y)|}{|x - y|} \mid x \neq y, x, y \in \Omega \right\}. \\
540 & \text{This manuscript is for review purposes only.}
$$
Proof. The proof relies on an induction principle: at stage k, $\varphi_0 \circ \cdots \circ \varphi_{k-1}$ is assumed to be a bi-Lipschitz homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$ such that $\|\nabla(\varphi_0 \circ \cdots \circ \varphi_{k-1})\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha$, $\| (\nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1}))^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta$, and $\det \nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1}) > 0$ a.e.. One proves that a minimiser of $(\mathcal{P}_k)$ exists and is also a bi-Lipschitz homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$. It rests on arguments similar to those previously used, among them the fact that bi-Lipschitz orientation-preserving homeomorphisms form a group stable for the composition. At last, it is composed of two parts: the first one is devoted to the existence of minimisers to problem $(\mathcal{P}_k)$ for fixed $k$, while the second one focuses on the asymptotic result.

• Existence of minimisers to problem $(\mathcal{P}_k)$ for fixed $k$

For $\varphi = \text{Id}$, we have that $\mathcal{F}_k(\varphi) := \mathcal{F}(\varphi_0 \circ \cdots \circ \varphi_{k-1} \circ \varphi, \sum_{j=0}^k T_j, \sum_{j=0}^k R_j) < +\infty$ (the uniform bound in $n$ may depend on $k$), and an inequality of coercivity holds, showing that the infimum is finite.

Let $(\varphi_{k,n})_{n \in \mathbb{N}}$ be a minimising sequence. We then set $\varphi_{k,n} = \varphi_0 \circ \varphi_1 \circ \cdots \circ \varphi_{k-1} \circ \varphi_{k,n}$ so that for $n$ large enough

$$ 
\begin{aligned}
& \frac{\lambda}{2} \left\| \left( \sum_{j=0}^{k} T_j \right) \circ \phi_{k,n} - \sum_{j=0}^{k} R_j \right\|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \phi_{k,n}, \det \nabla \phi_{k,n}) dx \\
& + 1 \{ \| \cdot \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha \} (\nabla \phi_{k,n}) + 1 \{ \| \cdot \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta \} ((\nabla \phi_{k,n})^{-1}) \\
& \le \inf_{\varphi \in \mathcal{X}_k} \mathcal{F}_k(\varphi) + 1 < +\infty.
\end{aligned}
$$

Applying again the generalised Poincaré inequality and Gagliardo-Nirenberg interpolation inequalities, it follows that $(\phi_{k,n})_{n \in \mathbb{N}}$ is uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$ and thus

$$ \phi_{k,n} \underset{n \to +\infty}{\stackrel{*}{\longrightarrow}} \bar{\phi}_k \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2). $$

Reasoning as in the proof of Theorem 2.3, one shows that $\phi_{k,n}$ is a homeomorphism of $\bar{\Omega}$ onto $\bar{\Omega}$ and $\bar{\phi}_{k,n}^{-1} \in W^{1,q}(\Omega, \mathbb{R}^2)$. The mapping $\phi_{k,n}$ is one-to-one almost everywhere, maps measurable sets in $\bar{\Omega}$ to measurable sets in $\bar{\Omega}$, and the change of variables formula holds. Also, $\|\nabla (\bar{\phi}_{k,n}^{-1})\|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \| (\nabla \phi_{k,n})^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))}$. Focusing now on $\bar{\phi}_k$ and arguing as before, $\det \nabla \bar{\phi}_k > 0$ a.e., $\frac{1}{\det \nabla \bar{\phi}_k} \in L^1(\Omega)$ and $\bar{\phi}_k^{-1} \in W^{1,q}(\Omega, \mathbb{R}^2)$ with $q \in [2, 11]$. At last, Rellich-Kondrachov theorem enables one to conclude that in fact $\bar{\phi}_k^{-1} \in W^{1,\infty}(\Omega, \mathbb{R}^2)$ with $\|\nabla (\bar{\phi}_k^{-1})\|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \| (\nabla \bar{\phi}_k)^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta$. Weak lower semi-continuity arguments among others yield

$$ 
\begin{aligned}
& \frac{\lambda}{2} \left\| \sum_{j=0}^{k} T_j \circ \bar{\phi}_k - \sum_{j=0}^{k} R_j \right\|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \bar{\phi}_k, \det \nabla \bar{\phi}_k) dx \\
& + 1 \{ \| \cdot \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha \} (\nabla \bar{\phi}_k) + 1 \{ \| \cdot \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta \} ((\nabla \bar{\phi}_k)^{-1}) = \inf_{\varphi \in \mathcal{X}_k} \mathcal{F}_k(\varphi).
\end{aligned}
$$

The mapping $\varphi_0 \circ \cdots \circ \varphi_{k-1}$ being a bi-Lipschitz orientation-preserving homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$ with $\text{Lip}(\varphi_0 \circ \cdots \circ \varphi_{k-1}, \Omega) \le \alpha$ and $\text{Lip}((\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1}, \Omega) \le \beta$.
578

579

$\beta$ from the previous stage of the induction process, setting
$\varphi_{k,n} = (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} \circ \phi_{k,n}$, one has,

580

581

582

583

584

585

586

587

588

589

590

591

592

$$
\begin{align*}
\mathrm{Lip}(\varphi_{k,n}, \Omega) &\le \mathrm{Lip}((\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1}, \Omega) \mathrm{Lip}(\phi_{k,n}, \Omega), \\
&\le \alpha\beta.
\end{align*}
$$

$\varphi_{k,n}$ is thus uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$ and there exists $\bar{\varphi}_k \in W^{1,\infty}(\Omega, \mathbb{R}^2)$ (in fact, $\mathrm{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2)$) such that $\varphi_{k,n} \xrightarrow[n \to +\infty]{} \bar{\varphi}_k$ in $W^{1,\infty}(\Omega, \mathbb{R}^2)$ and uniformly on $\bar{\Omega}$ up to a subsequence. At the same time, up to a subsequence (common subsequence), $(\phi_{k,n})_{n \in \mathbb{N}}$ uniformly converges to $\bar{\phi}_k$ on $\bar{\Omega}$ and thus $(\varphi_{k,n} = (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} \circ \phi_{k,n})_{n \in \mathbb{N}}$ uniformly converges to $(\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} \circ \bar{\phi}_k$ on $\bar{\Omega}$. By uniqueness of the limit, it follows that $\bar{\varphi}_k = (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} \circ \bar{\phi}_k$. Similar arguments to those previously used enable one to get a bound on $\mathrm{Lip}(\bar{\varphi}_k, \Omega)$, respectively $\mathrm{Lip}(\bar{\varphi}_k^{-1}, \Omega)$. Also, Corollary 2.21 of [15, p. 64] states that the derivatives of $\bar{\varphi}_k$ in the sense of distributions are given by the usual derivation formulas for composed functions. Thus

593

594

$$
\nabla \bar{\varphi}_k = \nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} (\bar{\phi}_k) \nabla \bar{\phi}_k \quad \text{a.e.,}
$$

595

yielding

596

597

$$
\det \nabla \bar{\varphi}_k = \det \nabla \bar{\phi}_k \det \left( \nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} (\bar{\phi}_k) \right).
$$

598

599

600

601

From the above, $\det \nabla \bar{\phi}_k > 0$ a.e. and $\det (\nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1}) > 0$ a.e.. Besides, let $\mathcal{N}_k \subset \bar{\Omega}$ be such that $\mathrm{meas}(\mathcal{N}_k) = 0$ and for all $x \in \Omega \setminus \mathcal{N}_k$, $\det (\nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1}) > 0$. Now let $\mathcal{N}_k' = \bar{\phi}_k^{-1}(\mathcal{N}_k)$. Then $\mathrm{meas}(\mathcal{N}_k') = 0$ and for every $y \notin \mathcal{N}_k'$,

602

603

$$
\det \left( \nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} (\bar{\phi}_k(y)) \right) = \det \left( \nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} (x) \right)
$$

604

605

606

607

608

609

610

611

612

with $x \in \Omega \setminus \mathcal{N}_k$, yielding $\det (\nabla (\varphi_0 \circ \cdots \circ \varphi_{k-1})^{-1} (\bar{\phi}_k)) > 0$ a.e. and subsequently $\det \nabla \bar{\varphi}_k > 0$ a.e..

• Asymptotic analysis

Let us now consider the sequence $(\varphi_k)_{k \in \mathbb{N}}$ of minimisers —$\varphi_k$ a minimiser of $(\mathcal{P}_k)$— and let $\bar{\varphi} \in \mathcal{W}$ be a minimiser of $\mathcal{F}(\cdot)$. Let us set $\phi_k = \varphi_0 \circ \cdots \circ \varphi_k$.

Since $\forall k \in \mathbb{N}$, $\|\nabla \phi_k\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \alpha$, invoking the generalised Poincaré inequality and Gagliardo-Nirenberg interpolation inequalities yields that $(\phi_k)_{k \in \mathbb{N}}$ is uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$. Thus there exist a subsequence still denoted by $(\phi_k)_{k \in \mathbb{N}}$ and $\bar{\phi} \in W^{1,\infty}(\Omega, \mathbb{R}^2)$ such that

613

614

$$
\phi_k \xrightarrow{k \to +\infty} \bar{\phi} \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2).
$$
By continuity of the trace operator and since $\phi_k \in \text{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2)$ by construction, one gets that $\bar{\phi} \in \text{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2)$.

From [10, Theorem 8.20], we have that $\det \nabla \phi_k \underset{k \to +\infty}{\to} \det \nabla \bar{\phi}$ in $L^2(\Omega)$ for all $p \in ]1, \infty[$. We also know that $W_{Op}$ is continuous and convex. If $\psi_n \underset{n \to +\infty}{\to} \bar{\psi}$ in $W^{1,4}(\Omega, \mathbb{R}^2)$, thus $\nabla \psi_n \underset{n \to +\infty}{\to} \nabla \bar{\psi}$ in $L^4(\Omega, M_2(\mathbb{R}))$ and we can extract a subsequence such that $\nabla \psi_n \underset{n \to +\infty}{\to} \nabla \bar{\psi}$ almost everywhere in $\Omega$. If $\kappa_n \underset{n \to +\infty}{\to} \bar{\kappa}$ in $L^2(\Omega)$, then one can extract a subsequence still denoted by $(\kappa_n)$ such that $\kappa_n \underset{n \to +\infty}{\to} \bar{\kappa}$ almost everywhere in $\Omega$. Then by applying Fatou's lemma, we get

$$ \liminf_{n \to +\infty} \int_{\Omega} W_{Op}(\nabla \psi_n, \kappa_n) \, dx \geq \int_{\Omega} W_{Op}(\nabla \bar{\psi}, \bar{\kappa}) \, dx. $$

As $W_{Op}$ is convex, so is $(\xi, \kappa) \mapsto \int_{\Omega} W_{Op}(\xi, \kappa) \, dx$ and we can apply [6, Corollaire III.8], [7] so that $\int_{\Omega} W_{Op}(\xi, \kappa) \, dx$ is also weakly lower semi-continuous in $L^4(\Omega, M_2(\mathbb{R})) \times L^2(\Omega)$ yielding

$$ \int_{\Omega} W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) \, dx \leq \liminf_{k \to +\infty} \int_{\Omega} W_{Op}(\nabla \phi_k, \det \nabla \phi_k) \, dx. $$

Now, again, due in particular to the property of stability by composition, $\forall k \in \mathbb{N}^*$, $\varphi_{k-1}^{-1} \circ \dots \circ \varphi_0^{-1} \circ \bar{\varphi} \in \mathcal{X}_k$ and by definition of $\varphi_k$,

$$ \begin{aligned} \mathcal{F}(\varphi_0 \circ \dots \circ \varphi_k, \sum_{j=0}^k T_j, \sum_{j=0}^k R_j) &\le \mathcal{F}(\bar{\varphi}, \sum_{j=0}^k T_j, \sum_{j=0}^k R_j), \\ &= \frac{\lambda}{2} \left| \left(\sum_{j=0}^k T_j\right) \circ \bar{\varphi} - \left(\sum_{j=0}^k R_j\right) \right|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \bar{\varphi}, \det \nabla \bar{\varphi}) \, dx \\ &\quad + \mathbb{1}\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\} (\nabla \bar{\varphi}) + \mathbb{1}\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\} ((\nabla \bar{\varphi})^{-1}). \end{aligned} $$

Since $(\sum_{j=0}^k T_j) \underset{k \to +\infty}{\to} T$ and $(\sum_{j=0}^k R_j) \underset{k \to +\infty}{\to} R$ in $L^2(\Omega)$ and $\bar{\varphi}$ is a minimiser of $\mathcal{F}$ —meaning that $\bar{\varphi}$ is a bi-Lipschitz homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$ with $\det \nabla \bar{\varphi} > 0$ a.e., $\|\nabla \bar{\varphi}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha$ and $\|(\nabla \bar{\varphi})^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta$, we get, applying the classical change of variable and $C = C(\beta) > 0$ denoting a constant depending only on $\beta$,

$$ \begin{aligned} & \|(\sum_{j=0}^k T_j) \circ \bar{\varphi} - (\sum_{j=0}^k R_j) - T \circ \bar{\varphi} + R\|_{L^2(\Omega)} \\ & \le \|(\sum_{j=0}^k T_j) \circ \bar{\varphi} - T \circ \bar{\varphi}\|_{L^2(\Omega)} + \|(\sum_{j=0}^k R_j) - R\|_{L^2(\Omega)}, \\ & \le C \|(\sum_{j=0}^k T_j) - T\|_{L^2(\Omega)} + \|(\sum_{j=0}^k R_j) - R\|_{L^2(\Omega)} \underset{k \to +\infty}{\to} 0. \end{aligned} $$
645

This results in:

646

$$
\begin{align*}
\liminf_{k \to +\infty} \mathcal{F}(\phi_k, \sum_{j=0}^k T_j, \sum_{j=0}^k R_j) & \\
&\le \liminf_{k \to +\infty} \frac{\lambda}{2} \| (\sum_{j=0}^k T_j) \circ \bar{\phi} - (\sum_{j=0}^k R_j) \|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) \, dx \\
&\quad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \bar{\phi}) + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}} ((\nabla \bar{\phi})^{-1}), \\
&= \lim_{k \to +\infty} \frac{\lambda}{2} \| (\sum_{j=0}^k T_j) \circ \bar{\phi} - (\sum_{j=0}^k R_j) \|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) \, dx \\
&\quad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \bar{\phi}) + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}} ((\nabla \bar{\phi})^{-1}), \\
&= \frac{\lambda}{2} \| T \circ \bar{\phi} - R \|_{L^2(\Omega)}^2 + \int_{\Omega} W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) \, dx \\
&\quad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \bar{\phi}) + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}} ((\nabla \bar{\phi})^{-1}), \\
&= \inf_{\varphi \in \mathcal{W}} \mathcal{F}(\varphi) < +\infty,
\end{align*}
$$

647

648

649

650

651

652

653

654

655

this latter quantity being independent of $k$. We thus deduce that

656

$$
\begin{align*}
\int_{\Omega} W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) \, dx &\le \liminf_{k \to +\infty} \int_{\Omega} W_{Op}(\nabla \phi_k, \det \nabla \phi_k) \, dx \\
&\le \liminf_{k \to +\infty} \mathcal{F}(\phi_k, \sum_{j=0}^k T_j, \sum_{j=0}^k R_j), \\
&\le \inf_{\varphi \in \mathcal{W}} \mathcal{F}(\varphi) < +\infty.
\end{align*}
$$

657

658

659

660

661

662

663

664

665

Since $W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) = +\infty$ where $\det \nabla \bar{\phi} \le 0$, the set on which it occurs must be of null measure, otherwise we would have $\int_{\Omega} W_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) \, dx = +\infty$. Consequently, $\det \nabla \bar{\phi} > 0$ almost everywhere in $\Omega$. Also, by the weak-$\ast$ lower semicontinuity of $\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))}$, $\|\nabla \bar{\phi}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \liminf_{k \to +\infty} \|\nabla \phi_k\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha$, $\mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \bar{\phi}) \le \liminf_{k \to +\infty} \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \phi_k)$, and for any $q \in ]2, 11]$,

666

$$
\begin{align*}
\int_{\Omega} \|(\nabla \bar{\phi})^{-1}\|^q \det \nabla \bar{\phi} \, dx &\le \int_{\Omega} \frac{1}{(\det \nabla \bar{\phi})^q} \|\nabla \bar{\phi}\|_F^q \det \nabla \bar{\phi} \, dx \\
&\le \alpha^q \frac{1}{\det \nabla \bar{\phi}} \|\nabla \bar{\phi}\|_{L^{q-1}(\Omega)}^{q-1}, \\
&< +\infty,
\end{align*}
$$

667

668

669

670

671

672

owing to the fact that $\frac{1}{\det \nabla \bar{\phi}} \in L^{10}(\Omega)$. Ball's theorems ([4, Theorems 1 and 2]) enable one to conclude that $\bar{\phi}$ is an homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$ with $\bar{\phi}^{-1} \in W^{1,q}(\Omega, \mathbb{R}^2)$. (Recall that $\bar{\phi} \in \text{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2)$). Again $((\nabla \phi_k)^{-1})_{k \in \mathbb{N}}$ is uniformly bounded with respect to $k$ in $L^\infty(\Omega, M_2(\mathbb{R}))$, which, combined with the generalised
Poincaré inequality, Gagliardo-Nirenberg interpolation inequalities and Ball's theorems, shows that $\phi_k$ is a bi-Lipschitz orientation-preserving homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$, $\phi_k^{-1}$ being uniformly bounded in $W^{1,\infty}(\Omega, \mathbb{R}^2)$. One can thus extract a subsequence, still denoted by $(\phi_k^{-1})_{k \in \mathbb{N}}$ such that

$$ \phi_k^{-1} \underset{k \to +\infty}{\overset{*}{\to}} \bar{u} \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2), $$

and, up to a subsequence, $(\phi_k^{-1})_{k \in \mathbb{N}}$ uniformly converges to $\bar{u}$ in $\bar{\Omega}$. Arguing as before, we prove that $\bar{u} = \bar{\phi}^{-1}$ by demonstrating that $\bar{u} \circ \bar{\phi} = \text{Id}$ everywhere on $\bar{\Omega}$. The mapping $\bar{\phi}$ is thus a bi-Lipschitz orientation-preserving homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$, and again, applying the same arguments as before, $\|(\nabla \bar{\phi})^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} = \|\nabla \bar{\phi}^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \liminf_{k \to +\infty} \|\nabla \phi_k^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} = \liminf_{k \to +\infty} \|(\nabla \phi_k)^{-1}\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \beta$ and $\mathbb{1}\{\|\cdot\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \beta\}((\nabla \bar{\phi})^{-1}) \le \liminf_{k \to +\infty} \mathbb{1}\{\|\cdot\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \beta\}((\nabla \phi_k)^{-1})$.

We now prove that $\|\sum_{j=0}^k T_j \circ \phi_k - \sum_{j=0}^k R_j\|_{L^2(\Omega)} \xrightarrow{k \to +\infty} \|T \circ \bar{\phi} - R\|_{L^2(\Omega)}$, using the first triangle inequality and the classical change of variable. In that purpose, one has

$$ \begin{aligned} & \|(\sum_{j=0}^k T_j) \circ \phi_k - \sum_{j=0}^k R_j - T \circ \bar{\phi} + R\|_{L^2(\Omega)} \\ & \le \|(\sum_{j=0}^k T_j) \circ \phi_k - T \circ \bar{\phi}\|_{L^2(\Omega)} + \|\sum_{j=0}^k R_j - R\|_{L^2(\Omega)}, \\ & \le \|(\sum_{j=0}^k T_j) \circ \phi_k - T \circ \phi_k\|_{L^2(\Omega)} + \|T \circ \phi_k - T \circ \bar{\phi}\|_{L^2(\Omega)} + \|\sum_{j=0}^k R_j - R\|_{L^2(\Omega)}, \\ & \le C(\beta) \|(\sum_{j=0}^k T_j) - T\|_{L^2(\Omega)} + \|T \circ \phi_k - T \circ \bar{\phi}\|_{L^2(\Omega)} + \|\sum_{j=0}^k R_j - R\|_{L^2(\Omega)}. \end{aligned} $$

As $\sum_{j=0}^k T_j \xrightarrow{k \to +\infty} T$ and $\sum_{j=0}^k R_j \xrightarrow{k \to +\infty} R$ in $L^2(\Omega)$ ([31]), it suffices to show that $\|T \circ \phi_k - T \circ \bar{\phi}\|_{L^2(\Omega)}$ converges to 0.

Let $(T_n) \in C_0^\infty(\Omega)$ be a sequence such that $T_n \xrightarrow{n \to +\infty} T$ in $L^2(\Omega)$ (property of density of $C_0^\infty(\Omega)$ in $L^2(\Omega)$). Let $\epsilon > 0$ be fixed. Let $N = N(\epsilon) \in \mathbb{N}$ be such that $\|T_N - T\|_{L^2(\Omega)} \le \epsilon/3$ and let $L_{N(\epsilon)}$, be a Lipschitz constant associated to $T_N$. According to the Sobolev embedding theorem, there exists $K \in \mathbb{N}$ such that $\forall k \in \mathbb{N}$,
698

$$ (k \ge K \Rightarrow \|\phi_k - \bar{\phi}\|_{L^2(\Omega, \mathbb{R}^2)} \le \frac{\epsilon}{3L_{N(\epsilon)}}). \text{ Let us take } k \ge K \text{ so that} $$

699

$$ \|T \circ \phi_k - T \circ \bar{\phi}\|_{L^2(\Omega)} \le \|T \circ \phi_k - T_N \circ \phi_k\|_{L^2(\Omega)} + \|T_N \circ \phi_k - T_N \circ \bar{\phi}\|_{L^2(\Omega)} \\ + \|T_N \circ \bar{\phi} - T \circ \bar{\phi}\|_{L^2(\Omega)}, $$

700

701

$$ \le 2C(\beta) \|T - T_N\|_{L^2(\Omega)} + L_N \|\phi_k - \bar{\phi}\|_{L^2(\Omega)}, $$

$$ \le \frac{2C(\beta)\epsilon}{3} + \frac{\epsilon}{3}. $$

702

703

We thus have proved that $\exists K = K(\epsilon) \in \mathbb{N}$ such that $\forall k \in \mathbb{N}$, one has

704

705

$$ (k \ge K(\epsilon) \Rightarrow \|T \circ \phi_k - T \circ \bar{\phi}\|_{L^2(\Omega)} \le \frac{2C(\beta)\epsilon}{3} + \frac{\epsilon}{3}). \text{ Then} $$

706

707

$$ \|(\sum_{j=0}^{k} T_j) \circ \phi_k - \sum_{j=0}^{k} R_j - T \circ \bar{\phi} + R\|_{L^2(\Omega)} \xrightarrow{k \to +\infty} 0. $$

708

By gathering all the results, we finally get that

709

710

$$ \mathcal{F}(\bar{\phi}) = \frac{\lambda}{2} \|T \circ \bar{\phi} - R\|_{L^2(\Omega)}^2 + \int_{\Omega} \mathcal{W}_{Op}(\nabla \bar{\phi}, \det \nabla \bar{\phi}) dx + 1_{\{\|\cdot\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \alpha\}} (\nabla \bar{\phi}) \\ + 1_{\{\|\cdot\|_{L^{\infty}(\Omega, M_2(\mathbb{R}))} \le \beta\}} ((\nabla \bar{\phi})^{-1}), $$

711

712

$$ \le \liminf_{k \to +\infty} \mathcal{F}(\phi_k, \sum_{j=0}^{k} T_j, \sum_{j=0}^{k} R_j) \le \mathcal{F}(\bar{\phi}) = \inf_{\varphi \in \mathcal{W}} \mathcal{F}(\varphi), $$

713

714

with $\bar{\phi} \in \mathcal{W}$ and $\phi_k \xrightarrow{k \to +\infty} \bar{\phi}$ in $W^{1,\infty}(\Omega, \mathbb{R}^2)$. Thus $\bar{\phi}$ is a minimiser of the initial problem formulated in $T$ and $R$, and the result is proved. ■

715

716

**3. Numerical Resolution.** We assume that the $(k-1)$-th stage is reached and we aim to numerically solve problem $(\mathcal{P}_k)$ for a fixed $k$.

717

718

719

720

721

722

723

724

725

726

**3.1. Motivations.** Problem $(\mathcal{P}_k)$ falls within the non-convex and non-differentiable class of optimisation problems which is the hardest one to solve numerically, due to the nonlinearity on both the deformation and its Jacobian as well as the $L^\infty$ penalties on the Jacobian deformation and its inverse. We therefore adopt a common strategy in nonlinear elasticity which consists in introducing auxiliary variables to lift the nonlinearity from Jacobian deformation to a new variable and to move the nonconvexity from the Jacobian deformation and its inverse to new variables. We adjust it to the registration setting following [12]. The underlying idea is to obtain either a non-convex differentiable problem or a convex non-differentiable problem in each variable which are more tractable from a computational point of view.

727

728

729

730

**3.2. Decoupled problem.** We therefore introduce multiple auxiliary variables: (i) $\phi$ simulates the composition of deformations at scale $k$, i.e. $\phi \approx \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1} \circ \varphi$ to deal with the nonlinearity coming from the composition with the Template image at scale $k$, (ii) $\psi$ mimics the inverse of $\phi$ to facilitate the handling of the inverse Jacobian deformation,
(iii) $V$ approximates the Jacobian of the composition of deformations that is $V \approx \nabla\phi$ to deal with the regularisation, (iv) $W$ reproduces the Jacobian of $\psi$.

Remark 3.1. To handle the $L^\infty$ penalty on $(\nabla\varphi_0\circ\varphi_1\circ\dots\circ\varphi_{k-1}\circ\varphi)^{-1}$ we take advantage of the following property: if $u$ is a homeomorphism from $\Omega$ into $\Omega$, and the inverse function $u^{-1}$ belongs to $W^{1,q}(\Omega, \mathbb{R}^2)$, the matrix of weak derivatives reads $\nabla(u^{-1}) = (\nabla u)^{-1}(u^{-1})$ ([4]) and the property proved in previous computations that for bi-Lipschitz homeomorphisms $u$, $\|\nabla(u^{-1})\|_{L^\infty(\Omega, M_2(\mathbb{R}))} = \|(\nabla u)^{-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))}$.

Let $(\gamma_{k,i})_{i\in\mathbb{N}}$ be an increasing sequence of positive real numbers such that $\lim_{i\to+\infty} \gamma_{k,i} = +\infty$ for a fixed $k$. We then derive a decoupled problem ($\mathcal{DP}_{k,i}$) using $L^p$-type penalties:

$$ 
\begin{aligned}
\inf_{\varphi,\phi,\psi,V,W} \{ & \mathcal{F}_{k,i}(\varphi,\phi,\psi,V,W) = \frac{\lambda}{2} \int_{\Omega} \left( \sum_{j=0}^{k} T_j - \sum_{j=0}^{k} R_j \circ \psi \right)^2 \det \nabla \psi \, dx + \int_{\Omega} W_{Op}(V, \det V) \, dx \\
& \qquad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(V) + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}}(W) \quad (\mathcal{DP}_{k,i}) \\
& \qquad + \frac{\gamma_{k,i}}{4} \|V - \nabla\phi\|_{L^4(\Omega, M_2(\mathbb{R}))}^4 + \frac{\gamma_{k,i}}{4} \|W - \nabla\psi\|_{L^4(\Omega, M_2(\mathbb{R}))}^4 \\
& \qquad + \frac{\gamma_{k,i}}{2} \|\zeta_{k-1}^{-1} \circ \phi - \varphi\|_{L^2(\Omega, \mathbb{R}^2)}^2 + \frac{\gamma_{k,i}}{2} \|\psi \circ \phi - \mathrm{Id}\|_{L^2(\Omega, \mathbb{R}^2)}^2 \},
\end{aligned}
$$

where we have set $\zeta_{k-1} = \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_{k-1}$, bi-Lipschitz homeomorphism from $\bar{\Omega}$ to $\bar{\Omega}$ with $\zeta_{k-1} \in \mathrm{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^2)$, $\det \zeta_{k-1} > 0$ a.e., $\|\nabla \zeta_{k-1}\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha$ and $\| (\nabla \zeta_{k-1})^{-1} \|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta$, and with $\varphi \in L^2(\Omega, \mathbb{R}^2)$, $\phi \in \{u \in \mathrm{Id} + W_0^{1,4}(\Omega, \mathbb{R}^2)\}$, $\psi \in \{u \in \mathrm{Id} + W_0^{1,4}(\Omega, \mathbb{R}^2), \det \nabla u > 0 \text{ a.e.}\}$, $V \in \{u \in L^4(\Omega, M_2(\mathbb{R}))$, $(\det u)^{-1} \in L^{10}(\Omega)$, $\det u > 0 \text{ a.e.}$, $\|u\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}$ and $W \in \{u \in L^\infty(\Omega, M_2(\mathbb{R})), \|u\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}$. In a companion paper [13] in preparation, we prove that for fixed $k$ and for $i$ large enough, problem ($\mathcal{DP}_{k,i}$) constitutes a good approximation of problem ($\mathcal{P}_k$).

Remark 3.2. Little more regularity is assumed on $\bar{T}_k := \sum_{j=0}^k T_j$, namely $\bar{T}_k \in L^4(\Omega)$. With the prescribed functional space for $\Psi$,

$$ \int_{\Omega} (\bar{T}_k - \bar{R}_k \circ \Psi)^2 \det \nabla \Psi \, dx \le 2 \int_{\Omega} |\bar{T}_k|^2 \det \nabla \Psi \, dx + 2 \int_{\Omega} (\bar{R}_k \circ \Psi)^2 \det \nabla \Psi \, dx. $$

While Cauchy-Schwarz inequality guarantees that the first term is finite, Theorem 1 of [4] holds, ensuring that the classical change of variable formula applies to the second term which is well-defined.

Remark 3.3. The $L^4$-penalties on both $V$ and $W$ are used for theoretical purposes to ensure the Jacobian of $\phi$ and the Jacobian of $\psi$ are both in $L^4(\Omega, M_2(\mathbb{R}))$ which is needed for the asymptotic result derived hereafter. However in practice, the $L^4$-penalties are replaced by $L^2$ ones which is not too restrictive since the problem becomes discrete and all the norms turn to be equivalent.

**3.3. Numerical algorithm.** Since solving the decoupled problem ($\mathcal{DP}_{k,i}$) for fixed scale $k$ and large enough $\gamma_{k,i}$ gives a good approximation of a solution to the initial problem ($\mathcal{P}_k$) ([13]), we propose a numerical algorithm depicted in Algorithm 3.1 based on an alternating
minimisation scheme. That is, for each variable, we derive a more computationally tractable minimisation sub-problem by fixing the other ones. We now turn to the numerical details of each sub-problem. To make the reading more fluid, we remind the reader of the expression of the overall functional (the dependence on parameter $k$ is made explicit to enhance the fact that the resolution is done for each scale $k$, while that on index $i$ is omitted)

$$ 
\begin{aligned}
\inf_{\substack{\varphi_k, \phi_k, \\ \psi_k, V_k, W_k}} & \mathcal{F}(\varphi_k, \phi_k, \psi_k, V_k, W_k) \\
& = \frac{\lambda_k}{2} \int_{\Omega} (\bar{T}_k - \bar{R}_k \circ \psi_k)^2 \det \nabla \psi_k \, dx \\
& \quad + \int_{\Omega} W_{Op}(V_k, \det V_k) \, dx \\
& \quad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(V_k) \\
& \quad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}}(W_k) \\
& \quad + \frac{\gamma_{1,k}}{4} \|V_k - \nabla \phi_k\|_{L^4(\Omega, M_2(\mathbb{R}))}^4 \\
& \quad + \frac{\gamma_{2,k}}{4} \|W_k - \nabla \psi_k\|_{L^4(\Omega, M_2(\mathbb{R}))}^4 \\
& \quad + \frac{\gamma_{3,k}}{2} \|\zeta_{k-1}^{-1} \circ \phi_k - \varphi_k\|_{L^2(\Omega, \mathbb{R}^2)}^2 \\
& \quad + \frac{\gamma_{4,k}}{2} \|\psi_k \circ \phi_k - \mathrm{Id}\|_{L^2(\Omega, \mathbb{R}^2)}^2,
\end{aligned}
\quad (\mathcal{DP}_k) $$

and we explicitly state each resulting sub-problem by fixing all but one of the variables.

● **SUB-PROBLEM 1: OPTIMISATION OVER V.** For each scale $k$, the sub-problem in $V_k$ reads

$$ 
\begin{aligned}
\inf_{V_k} F(V_k) + Reg(V_k) = & \int_{\Omega} \left[ a_{1,k} \|V_k\|^4 + (\det V_k - 1)^2 + \frac{a_{3,k}}{(\det V_k)^{10}} \right] \, dx \\
& + \frac{\gamma_{1,k}}{2} \|V_k - \nabla \phi_k\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(V_k).
\end{aligned}
$$

**Remark 3.4.** Since in the discrete setting all the norms are equivalent, we have replaced here the $L^4$ penalty term by an $L^2$ one which is easier to handle from a numerical point of view.

This problem can be seen as the sum of a proper closed convex function

$$ Reg(\cdot) = \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(\cdot) $$

and a smooth function $F$, and as in [12], we use the simple iterative forward-backward splitting algorithm [18]:

$$ V_k^{n+1} = \mathrm{prox}_{\gamma_{1,k}Reg}(V_k^n - \gamma_{1,k}\nabla F(V_k^n)), $$

with $\mathrm{prox}_{\gamma_{1,k}Reg}(y) = \min_x \frac{1}{2}\|x - y\|_2^2 + \gamma_{1,k}Reg(y) = P_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(y)$, $P_C$ being the projection operator onto the convex set $C$.

● **SUB-PROBLEM 2: OPTIMISATION OVER W.** For each scale $k$, we solve the following minimisation problem

$$ 
\begin{aligned}
\inf_{W_k} \frac{\gamma_{2,k}}{2} \|W_k - \nabla \psi_k\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}}(W_k) = P_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(\nabla \psi_k).
\end{aligned}
$$
● **SUB-PROBLEM 3: OPTIMISATION OVER $\phi$**. For each scale $k$, the sub-problem in $\phi_k$ reads

$$ \inf_{\phi_k} \frac{\lambda_k}{2} \int_{\Omega} (\bar{T}_k \circ \phi_k - \bar{R}_k)^2 dx + \frac{\gamma_k}{2} \|V_k - \nabla \phi_k\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 + \frac{\gamma_{3,k}}{2} \|\zeta_{k-1}^{-1} \circ \phi_k - \varphi_k\|_{L^2(\Omega, \mathbb{R}^2)}^2 \\ + \frac{\gamma_{4,k}}{2} \|\psi_k \circ \phi_k - \mathrm{Id}\|_{L^2(\Omega, \mathbb{R}^2)}^2. $$

We then solve the associated Euler-Lagrange equation using an $L^2$-gradient flow scheme with an implicit Euler time stepping.

● **SUB-PROBLEM 4: OPTIMISATION OVER $\psi$**. For each scale $k$, the sub-problem in $\psi_k$ reads

$$ \inf_{\psi_k} \frac{\gamma_{2,k}}{2} \|W_k - \nabla \psi_k\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 + \frac{\gamma_{3,k}}{2} \|\psi_k \circ \phi_k - \mathrm{Id}\|_{L^2(\Omega, \mathbb{R}^2)}^2 $$

*Remark 3.5.* Numerically, to be tractable in practice, the fidelity term to the identity mapping is re-expressed by means of the change of variable formula.

We then solve the associated Euler-Lagrange equation using an $L^2$-gradient flow scheme with an implicit Euler time stepping.

● **SUB-PROBLEM 5: OPTIMISATION OVER $\varphi$**. For each scale $k$, the problem in $\varphi_k$ reads:

$$ \inf_{\varphi_k} \frac{\gamma_k}{2} \|\zeta_{k-1}^{-1} \circ \phi_k - \varphi_k\|_{L^2(\Omega, \mathbb{R}^2)}^2 $$

which has an explicit solution: $\varphi_k = \zeta_{k-1}^{-1} \circ \phi_k$. We emphasise that the $L^\infty$ penalties are applied componentwise in our algorithm.

We now test our method on both synthetic and real data from the medical imaging field. In all experiments, both the Template and Reference images are decomposed into 10 scales using Tadmor et al.'s algorithm ([31]) with parameters $\lambda_0 = 0.15$ and a prescribed number of iterations equal to 200, except for the toy example T-shape with texture where we consider only 8 scales.

**4. Numerical experiments.** This section is devoted to the analysis of numerical experiments: firstly, on a pair of 2 synthetic binary images 'device8-1' from the MPEG7 shape database (http://www.dabi.temple.edu) corresponding to T-shapes in order to (i) identify the mechanisms at work in the algorithm, (ii) assess the relevance of the results with regard to the intended objectives (in particular, the ability of the algorithm to model deformations capturing increasingly fine details while the scale grows). Then on real data stemming from the medical imaging domain: first on a slice of a 4DMRI sequence acquired during free
<caption>**Algorithm 3.1** Our Proposed Method ($L^\infty$ constraints applied componentwise)</caption>

1. Start from $\phi_{-1} \leftarrow \text{Id}, V_{11,-1} \leftarrow 1, V_{12,-1} \leftarrow 0, V_{21,-1} \leftarrow 0, V_{22,-1} \leftarrow 1,$
   $W_{11,-1} \leftarrow 1, W_{12,-1} \leftarrow 0, W_{21,-1} \leftarrow 0, W_{22,-1} \leftarrow 1, \psi_{-1} \leftarrow \text{Id}, \varphi_{-1} \leftarrow \text{Id},$
   and $\zeta_{-1} \leftarrow \text{Id}$;
2. Choose $N$, the number of scales.
3. Compute $(T_j)_{j=0,\dots,N}$ and $(R_j)_{j=0,\dots,N}$;
4. **for** $k = 0, \dots, N$:
5.     $\bar{T}_k \leftarrow \sum_{j=0}^k T_j$, and $\bar{R}_k \leftarrow \sum_{j=0}^k R_j$;
6.     $\phi_k \leftarrow \phi_{k-1}, V_{11,k} \leftarrow V_{11,k-1}, V_{12,k} \leftarrow V_{12,k-1}, V_{21,k} \leftarrow V_{21,k-1},$
        $V_{22,k} \leftarrow V_{22,k-1}, W_{11,k} \leftarrow W_{11,k-1}, W_{12,k} \leftarrow W_{12,k-1}, W_{21,k} \leftarrow W_{21,k-1},$
        $W_{22,k} \leftarrow W_{22,k-1}, \psi_k \leftarrow \psi_{k-1}, \varphi_k \leftarrow \text{Id}, \text{and } \zeta_{k-1} \leftarrow \phi_{k-1}$;
7.     **for** $l = 1, \dots, nbIter$:
8.         **for** each pixel:
9.             $V_{11,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \alpha\}}(V_{11,k} - (4\|V_k\|^2 V_{11,k} + 2a_2(\det V_k - 1))$
        10.             $V_{22,k} - \frac{10a_3V_{22,k}^2}{\det V_k^{11}} + \gamma_1(V_{11,k} - \frac{\partial \phi_{1,k}}{\partial x}))$;
11.             $V_{12,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \alpha\}}(V_{12,k} - (4\|V_k\|^2 V_{12,k} - 2a_2(\det V_k - 1))$
        12.             $V_{21,k} + \frac{10a_3V_{21,k}^2}{\det V_k^{11}} + \gamma_1(V_{12,k} - \frac{\partial \phi_{1,k}}{\partial y}))$;
13.             $V_{21,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \alpha\}}(V_{21,k} - (4\|V_k\|^2 V_{21,k} - 2a_2(\det V_k - 1))$
        14.             $V_{12,k} + \frac{10a_3V_{12,k}^2}{\det V_k^{11}} + \gamma_1(V_{21,k} - \frac{\partial \phi_{2,k}}{\partial x}))$;
15.             $V_{22,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \alpha\}}(V_{22,k} - (4\|V_k\|^2 V_{22,k} + 2a_2(\det V_k - 1))$
        16.             $V_{11,k} - \frac{10a_3V_{11,k}^2}{\det V_k^{11}} + \gamma_1(V_{22,k} - \frac{\partial \phi_{2,k}}{\partial y}))$;
17.             $W_{11,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \beta\}}(\frac{\partial \psi_{1,k}}{\partial x})$;
18.             $W_{12,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \beta\}}(\frac{\partial \psi_{1,k}}{\partial y})$;
19.             $W_{21,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \beta\}}(\frac{\partial \psi_{2,k}}{\partial x})$;
20.             $W_{22,k} \leftarrow \text{proj}_{\{\|\cdot\|_{L^\infty(\Omega)} \le \beta\}}(\frac{\partial \psi_{2,k}}{\partial y})$;
21.         **for** each pixel:
22.             Solve the Euler-Lagrange equation with respect to $\phi_k$ using an $L^2$ gradient flow with implicit Euler time stepping;
23.             Solve the Euler-Lagrange equation with respect to $\psi_k$ using an $L^2$ gradient flow with implicit Euler time stepping;
24.             $\varphi_k \leftarrow \zeta_{k-1}^{-1} \circ \phi_k$;
25.         **return** $\phi_k, \psi_k, V_{11,k}, V_{12,k}, V_{21,k}, V_{22,k}, W_{11,k}, W_{12,k}, W_{21,k}, W_{22,k}, \varphi_k, \bar{T}_k \circ \phi_k$;

835 breathing of the right lobe liver [32] http://vision.ee.ethz.ch/~organmot/chapter_download.

836 shtml. Second, on a slice of CINE cardiac sequence (courtesy of Caroline Petitjean).

837 We recall here that $\bar{T}_k$ represents the truncated decomposition of the Template image at

838 scale $k$, $\bar{R}_k$ stands for the truncated decomposition of the Reference image at scale $k$, $\phi_k$

839 denotes the composition of deformations at scale $k$, i.e. $\phi_k \approx \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k$, $\varphi_k$ is the

840 refined deformation obtained at scale $k$, and $\psi_k$ presents the inverse deformation at scale
$k$, i.e. $\psi_k \approx \phi_k^{-1} \approx (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k)^{-1}$.

Before presenting in depth the results of our multiscale model, we first discuss how to set
the parameters correctly.

**4.1. Parameter selection.** According to ($\mathcal{DP}_k$), 10 parameters are involved in the problem we numerically solve, and the chosen values for each experiment are reported in Table 1. Parameter $\lambda_k$ weighs the fidelity term at each scale $k$. When setting it, a trade-off must be met between accuracy of the alignment —requiring then high values of this parameter —and physically meaningful deformations —implying smaller values. The ranges are rather stable for each experiment as seen in Table 1 and go from 0.2 to 3. Parameters $a_{1,k}$, $a_{2,k}$, and $a_{3,k}$ involved in the Ogden stored energy function serving as part of the deformation regularisation to impose physical soundness affect respectively the average local change of length and the average local change of area at each scale, impacting subsequently the local rigidity of the deformations. The higher the $a_{i,k}$'s are, the more rigid the deformation is. These are rather stable for all experiments and with the scale growing, as one can see in Table 1. Parameters $\gamma_{1,k}$, $\gamma_{2,k}$, $\gamma_{3,k}$ and $\gamma_{4,k}$ are considered to be fixed for all scales and are chosen rather big as they ensure the closeness between the introduced auxiliary variables and those they are supposed to simulate as seen in the previous theoretical sections. $\alpha$ and $\beta$ are fixed for all experiments and all scales and ensure that the deformation Jacobian does not become too big. The choice for the number of scales $k$ considered follows from the discussion in [31] by using the following stopping criterion $\|u_k - u_{k+1}\|_{L^2(\Omega)} \le \delta$ with $\delta$ being a specified tolerance. Nevertheless, rather than the value of this parameter in itself, it seems to us that it is the combination of this parameter with the initial parameter $\lambda_0$ that is important since this latter dictates the level of detail contained in the images. To bound above the value of $k$, a data-driven preprocessing step based on [30] could be applied on both images to find the optimal value of $k$ or equivalently the optimal regularisation parameter $\lambda_0 2^{k+1}$ ensuring that features below a user-chosen threshold are removed.

**4.2. Evaluation protocol.** In order to quantitatively evaluate the accuracy and practicality of our model, in addition to a close and detailed visual inspection of the results, we consider the following metrics:

* the Dice coefficient [16] which measures set agreement (after binarising the images at each scale by thresholding). The closer it is to 1, the better the set agreement is and therefore the better the accuracy of the registration process is. A comparison of the Dice coefficient between $\bar{T}_k \circ \phi_k$ and $\bar{R}_k$, and the one between $\bar{T}_k$ and $\bar{R}_k$ allows us to quantitatively evaluate the quality of the registration at each scale. Then the Dice coefficient between $T$ and $R$ serves as a baseline to evaluate the improved quality of the registration process as the scale grows with the Dice coefficient between $T \circ \phi_k$ and $R$.

* $\min(\det(\nabla\phi_k))$ and $\max(\det(\nabla\phi_k))$ which ensure topology preservation of the global deformation at each scale. This range also indicates how far the deformation is from the volume preserving identity mapping and therefore quantifies the level of compression and dilation. The wider it is, the bigger the local compressions/dilations
<table><caption>**Table 1**
Selected parameters for the experiments.</caption><thead><tr><th>k</th><th>λ</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>γ<sub>3</sub></th><th>γ<sub>4</sub></th><th>α</th><th>β</th></tr></thead><tbody><tr><td colspan="11">Liver MRI</td></tr><tr><td>0</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>1</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>2</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>3</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>4</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>5</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>6</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>7</td><td>2</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>8</td><td>1.5</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>9</td><td>1.5</td><td>5</td><td>1000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td colspan="11">CINE cardiac MRI</td></tr><tr><td>0</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>1</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>2</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>3</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>4</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>5</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>6</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>7</td><td>2</td><td>8</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>8</td><td>2.5</td><td>8</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td>9</td><td>3</td><td>8</td><td>3500</td><td>4</td><td>80000</td><td>1</td><td>10</td><td>10</td><td>100</td><td>100</td></tr><tr><td colspan="11">T-shape</td></tr><tr><td>0</td><td>1</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>1</td><td>1</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>2</td><td>1</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>3</td><td>0.5</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>4</td><td>0.5</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>5</td><td>0.5</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>6</td><td>0.5</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>7</td><td>0.5</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>8</td><td>0.5</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>9</td><td>0.2</td><td>5</td><td>2000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td colspan="11">T-shape-texture</td></tr><tr><td>0</td><td>1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>1</td><td>1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>2</td><td>1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>3</td><td>1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>4</td><td>1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>5</td><td>0.1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>6</td><td>0.1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr><tr><td>7</td><td>0.1</td><td>5</td><td>3000</td><td>4</td><td>80000</td><td>1</td><td>1</td><td>1</td><td>100</td><td>100</td></tr></tbody></table>
are.

* $\text{Re-SSD}(\bar{T}_k \circ \phi_k, \bar{R}_k) = \frac{\|\bar{T}_k \circ \phi_k - \bar{R}_k\|^2}{\|\bar{T}_k - \bar{R}_k\|^2}$ and $\text{Re-SSD}(T \circ \phi_k, R) = \frac{\|T \circ \phi_k - R\|^2}{\|T - R\|^2}$ which measure the intensity alignment between the deformed truncated Template and the truncated Reference image at each scale and the initial Template image deformed by the transformation obtained at scale *k* and the initial Reference image. The closer it is to 0, the better the alignment is. The former one quantifies the quality of the registration process at each scale while the latter assesses the refined accuracy of the registration as the scale grows.

The results are reported for each experiment in Table 2, Table 3, Table 5 and Table 6 respectively. We now turn to the first synthetic numerical experiment.

**4.3. Toy example.** The proposed method is first evaluated on a synthetic example (Figure 2) to emphasise the ability of the model to generate large deformations and to handle noisy data. Indeed, white Gaussian noise is added to the Reference image to illustrate the benefit of multiscale image registration when dealing with noisy images (see Figure 1) while keeping a fidelity term based on intensity comparison. In Figure 2, we

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Template T-shape image on the left and noisy T-shape Reference image on the right.</figcaption></figure>

observe that on the first scales, the noise in the truncated Reference image is removed which allows our multiscale model to produce correct deformations without perturbations inherited from the noise. However, the first level is too blurry to generate accurate deformations, and more levels are needed to correct the displacements of the junctions between the vertical and horizontal bars of the T as seen in the zoom-in view at scale 9. This may also explain why the amplitude of the determinant range in Table 2 decreases as the scale increases: at the first scales, the hierarchical decomposition algorithm produces images with blurry and thus rough contours, which may entail larger deformations than the one required for pairing the original images. Additionally, it exemplifies the ‘corrective’ dimension of the proposed algorithm in the latter scales. We see that at the last level, the difference map between the deformed Template and the Reference is only composed of the noise, which supports the fact that this hierarchical decomposition of the deformations can help the registration process in case of noisy data. These observations are also supported by quantitative metrics in Table 2. In this case, the Dice coefficient is a much more reliable metric than the Re-SSD to evaluate the registration accuracy since Re-SSD compares the intensity values and is thus very sensitive to noise. We observe that the quantity Dice($T \circ \phi_k, R$) increases and gets closer to 1 as $k$ grows and therefore, several levels are needed to achieve the best registration accuracy.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Multiscale registration results on synthetic T-shape images with noise (size: 100×100, time: 7 minutes): each row represents a scale of the deformation; the first column displays the Template image at scale k, i.e. $\bar{T}_k$, the second column shows the Reference image at scale k, i.e. $\bar{R}_k$, the third one illustrates the deformed Template obtained at scale k, i.e. $\bar{T}_k \circ \phi_k$, the fourth one exhibits the absolute difference $|\bar{T}_k \circ \phi_k - \bar{R}_k|$ at scale k, the fifth column presents the inverse deformation at scale k, i.e. $\psi_k \approx \phi_k^{-1} \approx (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k)^{-1}$, the sixth column represents the composition of deformations at scale k, i.e. $\phi_k \approx \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k$, and finally the last column displays the deformation obtained at scale k, i.e. $\varphi_k$. A few scales have been removed to improve the readability of the figure.</figcaption></figure>
<table><caption>Table 2</caption><thead><tr><th>k</th><th>Dice(T<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T∘φ<sub>k</sub>,R)</th><th>min(det(∇φ<sub>k</sub>))</th><th>max(det(∇φ<sub>k</sub>))</th><th>Re-SSD(T∘φ<sub>k</sub>,R)</th><th>Re-SSD(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th></tr></thead><tbody><tr><td>0</td><td>0.4607</td><td>0.9621</td><td>0.9095</td><td>0.5991</td><td>1.4543</td><td>0.6773</td><td>0.6966</td></tr><tr><td>1</td><td>0.6168</td><td>0.9765</td><td>0.9176</td><td>0.6461</td><td>1.2582</td><td>0.6634</td><td>0.6828</td></tr><tr><td>2</td><td>0.6241</td><td>0.9649</td><td>0.9115</td><td>0.6024</td><td>1.2824</td><td>0.6825</td><td>0.6964</td></tr><tr><td>3</td><td>0.6076</td><td>0.9773</td><td>0.9205</td><td>0.7175</td><td>1.2023</td><td>0.6625</td><td>0.6810</td></tr><tr><td>4</td><td>0.5983</td><td>0.9820</td><td>0.9225</td><td>0.7434</td><td>1.2256</td><td>0.6591</td><td>0.6703</td></tr><tr><td>5</td><td>0.5888</td><td>0.9874</td><td>0.9267</td><td>0.7613</td><td>1.2562</td><td>0.6574</td><td>0.6622</td></tr><tr><td>6</td><td>0.5796</td><td>0.9896</td><td>0.9285</td><td>0.7834</td><td>1.2548</td><td>0.6546</td><td>0.6554</td></tr><tr><td>7</td><td>0.6203</td><td>0.9881</td><td>0.9293</td><td>0.7676</td><td>1.2991</td><td>0.6520</td><td>0.6540</td></tr><tr><td>8</td><td>0.6491</td><td>0.9810</td><td>0.9293</td><td>0.7727</td><td>1.3352</td><td>0.6505</td><td>0.6525</td></tr><tr><td>9</td><td>0.6037</td><td>0.9648</td><td>0.9303</td><td>0.8198</td><td>1.2173</td><td>0.6556</td><td>0.6569</td></tr></tbody></table>

Quantitative analysis of the multiscale registration model on synthetic T-shape images (for comparison, Dice(T,R) = 0.8227).

**4.4. Toy example T-shapes with texture.** In this example, two textures with different scales are added to the synthetic T-shape pair of images to illustrate the deformation decomposition mechanisms at work in our multiscale registration model. The first additional texture consists in fine vertical layers while the second one is composed of small black circles inside the T-shape. The number of layers and circles is the same in both images to comply with the topology preservation assumption of our model. However, the identity mapping on the boundary assumption is not fulfilled and we can see that even at the finer scale there is a small error in the alignment at the bottom of the image in Figure 3. We observe that on the first scales only the T-shape is kept and the deformation maps correctly match the two simplified images. From scale 2-3, the black circles are appearing and we see in the zoom-in views that the deformation $\phi_k$ obtained at this scale corresponds to the movement of these circles. The same phenomenon occurs at scale 4 with the apparition of the vertical lines. With the zoom-in view at scale 7, we see that after the apparition of all the textures, a corrective process is at hand to improve the accuracy of the registration process. This is further corroborated by the study of the quantitative metrics in Table 3. Indeed, both the Dice scores and the Re-SSD indicate that the matching is accurate at each scale and that as the scale grows the alignment between the original deformed Template and the original Reference improves significantly. We also notice a gap in the progression of these metrics at scale 3 when the circles appear and at scale 4 when the vertical lines emerge. It therefore shows the ability of our model to correct the deformations of small textures as the scale increases. As we add more deformations through the scales, we see that the range of the determinant Jacobian increases. Furthermore, it remains positive at all times and by looking at the deformation maps in Figure 3, we see that the produced deformations at each scale are topology preserving and invertible as requested.

We now present the numerical results on the first real medical dataset, that is the right lobe liver MRI.

**4.5. Liver.** We chose the images of size 195 × 166 corresponding to the liver in full inhalation and the liver in full inhalation. The goal is to illustrate the capability of our model to deal with large deformations and fine structures, and to refine deformations scale after scale. The results are illustrated in Figure 5. We notice that at each scale, the deformed truncated Template is well-aligned with the truncated Reference, showing the capability of our model to deal with large and complex deformations. At each level, the deformation
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Multiscale registration results on synthetic T-shape images with texture (size: 100×100, time: 4 minutes): each row represents a scale of the deformation; the first column displays the Template image at scale k, i.e. $\bar{T}_k$, the second column shows the Reference image at scale k, i.e. $\bar{R}_k$, the third one illustrates the deformed Template obtained at scale k, i.e. $\bar{T}_k \circ \phi_k$, the fourth one exhibits the absolute difference $|\bar{T}_k \circ \phi_k - \bar{R}_k|$ at scale k, the fifth column presents the inverse deformation at scale k, i.e. $\psi_k \approx \phi_k^{-1} \approx (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k)^{-1}$, the sixth column represents the composition of deformations at scale k, i.e. $\phi_k \approx \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k$, and finally the last column displays the deformation obtained at scale k, i.e. $\varphi_k$. A few scales have been removed to improve the readability of the figure.</figcaption></figure>

This manuscript is for review purposes only.
<figure><img src="image_4.png" /><figcaption>Figure 4. Multiscale registration results from [8] on synthetic T-shape images with texture.</figcaption></figure>

<table><caption>Table 3</caption><thead><tr><th>k</th><th>Dice(T<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>min(det(∇φ<sub>k</sub>))</th><th>max det(∇φ<sub>k</sub>))</th><th>Re-SSD(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>Re-SSD(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th></tr></thead><tbody><tr><td>0</td><td>0.4805</td><td>0.9848</td><td>0.8170</td><td>0.7674</td><td>1.5001</td><td>0.5357</td><td>0.0357</td></tr><tr><td>1</td><td>0.4856</td><td>0.9860</td><td>0.8142</td><td>0.7592</td><td>1.4674</td><td>0.5156</td><td>0.0391</td></tr><tr><td>2</td><td>0.4890</td><td>0.9870</td><td>0.8221</td><td>0.7600</td><td>1.4766</td><td>0.4321</td><td>0.0416</td></tr><tr><td>3</td><td>0.5390</td><td>0.9541</td><td>0.8931</td><td>0.6118</td><td>1.6482</td><td>0.2125</td><td>0.0652</td></tr><tr><td>4</td><td>0.5919</td><td>0.9478</td><td>0.9132</td><td>0.4290</td><td>1.7419</td><td>0.1463</td><td>0.0625</td></tr><tr><td>5</td><td>0.6107</td><td>0.9390</td><td>0.9177</td><td>0.4835</td><td>1.7345</td><td>0.1582</td><td>0.0941</td></tr><tr><td>6</td><td>0.6110</td><td>0.9360</td><td>0.9206</td><td>0.5260</td><td>1.7323</td><td>0.1582</td><td>0.1393</td></tr><tr><td>7</td><td>0.6252</td><td>0.9373</td><td>0.9230</td><td>0.5546</td><td>1.7342</td><td>0.1546</td><td>0.1595</td></tr></tbody></table>

Quantitative analysis of the multiscale registration model on synthetic T-shape images with texture (for comparison, Dice(T, R) = 0.7020).

<table><caption>Table 4</caption><thead><tr><th>k</th><th>Dice(T<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>Re-SSD(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>min(det(∇φ<sub>k</sub>))</th><th>max det(∇φ<sub>k</sub>))</th></tr></thead><tbody><tr><td>0</td><td>0.7826</td><td>0.8547</td><td>0.1911</td><td>0.7420</td><td>1.2366</td></tr><tr><td>1</td><td>0.7243</td><td>0.8777</td><td>0.1092</td><td>0.7010</td><td>1.2712</td></tr><tr><td>2</td><td>0.7063</td><td>0.9575</td><td>0.0529</td><td>0.6546</td><td>1.3111</td></tr></tbody></table>

Quantitative analysis of the multiscale registration model from [8] on synthetic T-shape images with texture (for comparison, Dice(T, R) = 0.7020).

948 grids do not exhibit overlaps and therefore confirm the theoretical topology preservation

949 property of our model, meaning that the produced deformations are physically relevant
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Multiscale registration results on liver MRI from full inhalation to full exhalation (size: 195×166, time: 78 minutes): each row represents a scale of the deformation; the first column displays the Template image at scale k, i.e. $\bar{T}_k$, the second column shows the Reference image at scale k, i.e. $\bar{R}_k$, the third one illustrates the deformed Template obtained at scale k, i.e. $\bar{T}_k \circ \phi_k$, the fourth one exhibits the absolute difference $|\bar{T}_k \circ \phi_k - \bar{R}_k|$ at scale k, the fifth column presents the inverse deformation at scale k, i.e. $\psi_k \approx \phi_k^{-1} \approx (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k)^{-1}$, the sixth column represents the composition of deformations at scale k, i.e. $\phi_k \approx \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k$, and finally the last column displays the deformation obtained at scale k, i.e. $\varphi_k$. A zoom on $\varphi_k$ is proposed to see its local nature as the scale grows.</figcaption></figure>
with a positive Jacobian determinant. One can also see, thanks to the last column, that the hierarchical decomposition of the deformations obtained with our model behaves as expected. That is, the deformations from the first scales are large and global, representing the movements of the main organs, i.e. the liver and the kidney, meanwhile as the scale grows, the deformation becomes more localised and refined to model the motion of small features or structures inside the organs, i.e. the blood vessels here imaged as white dots. The zoom-in view helps to see this refined transformation.

We observe that at each scale, the Jacobian determinant remains positive, which supports

<table><caption>**Table 5**</caption><thead><tr><th>k</th><th>Dice(T<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T∘φ<sub>k</sub>,R)</th><th>min(det(∇φ<sub>k</sub>))</th><th>max det(∇φ<sub>k</sub>))</th><th>Re-SSD(T∘φ<sub>k</sub>,R)</th><th>Re-SSD(T<sub>k</sub>∘φ<sub>k</sub>,R<sub>k</sub>)</th></tr></thead><tbody><tr><td>0</td><td>0.7009</td><td>0.9850</td><td>0.9144</td><td>0.8063</td><td>1.5046</td><td>0.6024</td><td>0.4756</td></tr><tr><td>1</td><td>0.7115</td><td>0.980</td><td>0.9165</td><td>0.7729</td><td>1.5282</td><td>0.5401</td><td>0.4626</td></tr><tr><td>2</td><td>0.7164</td><td>0.9838</td><td>0.9189</td><td>0.6617</td><td>1.5332</td><td>0.5022</td><td>0.4709</td></tr><tr><td>3</td><td>0.7209</td><td>0.9802</td><td>0.9211</td><td>0.5785</td><td>1.5374</td><td>0.4793</td><td>0.4768</td></tr><tr><td>4</td><td>0.7281</td><td>0.9749</td><td>0.9265</td><td>0.5729</td><td>1.5454</td><td>0.4498</td><td>0.4698</td></tr><tr><td>5</td><td>0.7323</td><td>0.9700</td><td>0.9301</td><td>0.5603</td><td>1.5605</td><td>0.4392</td><td>0.4595</td></tr><tr><td>6</td><td>0.7349</td><td>0.9623</td><td>0.9335</td><td>0.5572</td><td>1.5940</td><td>0.4320</td><td>0.4476</td></tr><tr><td>7</td><td>0.7303</td><td>0.9539</td><td>0.9347</td><td>0.5632</td><td>1.6538</td><td>0.4338</td><td>0.4434</td></tr><tr><td>8</td><td>0.7283</td><td>0.9504</td><td>0.9366</td><td>0.5752</td><td>1.6342</td><td>0.4357</td><td>0.4397</td></tr><tr><td>9</td><td>0.7271</td><td>0.9452</td><td>0.9375</td><td>0.5604</td><td>1.6368</td><td>0.4429</td><td>0.4451</td></tr></tbody></table>

Quantitative analysis of the multiscale registration model on free-breathing liver MRI images (for comparison, Dice(T,R) = 0.7511).

our previous claim that our model produces physically meaningful and reliable deformations. Also, at each scale, Dice($\bar{T}_k \circ \phi_k$, $\bar{R}_k$) is always greater than Dice($\bar{T}_k$, $\bar{R}_k$) and close to one, and Re-SSD($\bar{T}_k \circ \phi_k$, $\bar{R}_k$) is close to 0 which reinforces the quality and accuracy of the registration process at each scale. Finally, one can see that as the scale grows, Dice($T \circ \phi_k$, R) becomes closer to 1 and Re-SSD($T \circ \phi_k$, R) closer to 0 which can be interpreted as the alignment refinement through the scales. Indeed, in Figure 5, we see that as the scale grows, more localised and refined deformations are added to correct for small features displacements. This is further justified by the tendency of the determinant range to widen as the scale increases.

**4.6. Cardiac MRI.** We were supplied with a whole cardiac MRI examination of a patient (courtesy of Caroline Petitjean from the LITIS, University of Rouen Normandie, France). It is made of 280 images divided into 14 levels of slice and 20 images per cardiac cycle of size 150 × 150. A cardiac cycle is composed of a contraction phase (40% of the cycle duration), followed by a dilation phase (60% of the cycle duration). In order to assess the accuracy of the proposed algorithm in handling large and nonlinear deformations, we propose to register a pair of the type: Reference corresponding to end diastole (ED), that is when the heart is the most dilated, and Template corresponding to end systole (ES), that is when the heart is the most contracted. This corresponds to the results depicted in Figure 6. We can see that at each level the deformed truncated Template and the truncated Reference images are well-matched. This visual inspection is confirmed by the Dice coefficients Dice($\bar{T}_k \circ \phi_k$, $\bar{R}_k$) close to one and Re-SSD($\bar{T}_k \circ \phi_k$, $\bar{R}_k$) close to 0 in Table 6. Our algorithm also gives us as outputs the global deformation grids, the global inverse deformation grids and the refined deformation grids at each scale, plotted respectively in column 5, 6 and 7 of Figure 6. We see that none of them exhibit overlaps meaning that the
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Multiscale registration results on CINE cardiac MRI from the end of systole to the end of diastole (size: 150×150, time: 21 minutes): each row represents a scale of the deformation; the first column displays the Template image at scale k, i.e. $\bar{T}_k$, the second column shows the Reference image at scale k, i.e. $\bar{R}_k$, the third one illustrates the deformed Template obtained at scale k, i.e. $\bar{T}_k \circ \phi_k$, the fourth one exhibits the absolute difference $|\bar{T}_k \circ \phi_k - \bar{R}_k|$ at scale k, the fifth column presents the inverse deformation at scale k, i.e. $\psi_k \approx \phi_k^{-1} \approx (\varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k)^{-1}$, the sixth column represents the composition of deformations at scale k, i.e. $\phi_k \approx \varphi_0 \circ \varphi_1 \circ \dots \circ \varphi_k$, and finally the last column displays the deformation obtained at scale k, i.e. $\varphi_k$. A zoom on $\varphi_k$ is proposed to see its local nature as the scale grows. A few scales have been removed to improve the readability of the figure.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Multiscale registration results from [8] on CINE cardiac MRI images.</figcaption></figure>

produced deformations are invertible and preserve topology. This is corroborated by the fact that the Jacobian determinant of the global deformation remains positive at each scale in Table 6. Finally, the last column displays the deformation produced at each scale and we notice that in the first scales the deformation is global encompassing the movements of the main parts of the heart while as the scale grows the deformations become more localised and refined and correspond to the motion of small black structures as shown in the zoom-in views. This allows to correct the registration to capture smaller displacements and subsequently, to improve the matching accuracy of the initial Template with the initial Reference. This is highlighted by the fact that $\text{Dice}(T \circ \phi_k, \mathbb{R})$ grows and gets closer to 1 as $k$ increases and that the $\text{Re-SSD}(T \circ \phi_k, \mathbb{R})$ decreases and becomes closer to 0 in Table 6.
Here again, the determinant minimum tends to decrease, while the determinant maximum increases as the scale becomes bigger, which means that locally the contractions/expansions are bigger to correct the alignment between the original Template and Reference.

<table><caption>**Table 6**</caption><thead><tr><th>k</th><th>Dice(T<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T<sub>k</sub> ◦ φ<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T ◦ φ<sub>k</sub>,R)</th><th>min(det(∇φ<sub>k</sub>))</th><th>max(det(∇φ<sub>k</sub>))</th><th>Re-SSD(T ◦ φ<sub>k</sub>,R)</th><th>Re-SSD(T<sub>k</sub> ◦ φ<sub>k</sub>,R<sub>k</sub>)</th></tr></thead><tbody><tr><td>0</td><td>0.3202</td><td>0.8504</td><td>0.8528</td><td>0.6398</td><td>1.9041</td><td>0.3124</td><td>0.3419</td></tr><tr><td>1</td><td>0.3489</td><td>0.8837</td><td>0.8646</td><td>0.6248</td><td>1.9174</td><td>0.2808</td><td>0.2957</td></tr><tr><td>2</td><td>0.3589</td><td>0.8917</td><td>0.8693</td><td>0.6110</td><td>1.8978</td><td>0.2709</td><td>0.2729</td></tr><tr><td>3</td><td>0.3675</td><td>0.8947</td><td>0.8746</td><td>0.2611</td><td>1.9012</td><td>0.2611</td><td>0.2591</td></tr><tr><td>4</td><td>0.3774</td><td>0.9001</td><td>0.8824</td><td>0.5857</td><td>1.9129</td><td>0.2500</td><td>0.2444</td></tr><tr><td>5</td><td>0.3833</td><td>0.9007</td><td>0.8895</td><td>0.5852</td><td>1.9324</td><td>0.2365</td><td>0.2314</td></tr><tr><td>6</td><td>0.3893</td><td>0.9033</td><td>0.8946</td><td>0.5848</td><td>1.9551</td><td>0.2261</td><td>0.2219</td></tr><tr><td>7</td><td>0.3936</td><td>0.9014</td><td>0.8978</td><td>0.5761</td><td>1.9706</td><td>0.2203</td><td>0.2175</td></tr><tr><td>8</td><td>0.3952</td><td>0.9035</td><td>0.8992</td><td>0.5535</td><td>1.9567</td><td>0.2178</td><td>0.2164</td></tr><tr><td>9</td><td>0.3959</td><td>0.9049</td><td>0.9049</td><td>0.5137</td><td>2.0167</td><td>0.2086</td><td>0.2074</td></tr></tbody></table>

Quantitative analysis of the multiscale registration model on CINE cardiac MRI images (for comparison, Dice(T,R) = 0.4059).

<table><caption>**Table 7**</caption><thead><tr><th>k</th><th>Dice(T<sub>k</sub>,R<sub>k</sub>)</th><th>Dice(T<sub>k</sub> ◦ φ<sub>k</sub>,R<sub>k</sub>)</th><th>Re-SSD(T<sub>k</sub> ◦ φ<sub>k</sub>,R<sub>k</sub>)</th><th>min(det(∇φ<sub>k</sub>))</th><th>max(det(∇φ<sub>k</sub>))</th></tr></thead><tbody><tr><td>0</td><td>0.7391</td><td>0.8322</td><td>0.4066</td><td>0.5458</td><td>2.3574</td></tr><tr><td>1</td><td>0.7210</td><td>0.9183</td><td>0.1729</td><td>0.4848</td><td>2.5065</td></tr><tr><td>2</td><td>0.7217</td><td>0.9453</td><td>0.1321</td><td>0.4821</td><td>2.5957</td></tr><tr><td>3</td><td>0.7122</td><td>0.9431</td><td>0.1236</td><td>0.4776</td><td>2.7164</td></tr></tbody></table>

Quantitative analysis of the multiscale registration model from [8] on CINE cardiac MRI images (for comparison, Dice(T,R) = 0.4059).

**4.7. Comparative assessment with a well-established method.** In order to substantiate the relevancy of the proposed approach in comparison to well-established methods, a comparative assessment is carried out between our method and [8] which is based on hyperelasticity principles and multiresolution techniques. Indeed, a parallel can be drawn between our approach and multiresolution techniques as described in [21, Chapter 13, Section 6]. In this latter framework, starting with the coarsest resolution, the deformation pairing the two images is computed, generally requiring low computational costs. The coarse resolution deformation is then extended by interpolation on a finer grid, and serves as initial condition for the registration task at finer resolution. Apart from reducing the computational burden, this approach brings some regularisation as structural main changes are captured by all scales, while more subtle displacements are only encoded in finer scales, according to the same principle as the one at work in our model. However, these two methods differ in the meaning given to multiscale/multiresolution. While multiresolution techniques refer primarily to the size of the images manipulated (which are increasingly large at finer scales) and to the implemented interpolation technique to move from one scale to a finer one, our proposed work focuses more on the image content and on the scale of the features (the size of the images manipulated being the same through the process). It seems to us that it is more closely related to the information encoded by the images whereas multiresolution techniques remove arbitrarily some image data.

The FAIR code ([22]) (courtesy of Pr. Modersitzki (Institute of Mathematics and Image Computing, University of Lübeck) and Dr. Ruthotto (Emory University, Department of Mathematics)) is used, implementing the method developed in [8]. Our analysis focuses
on the two experiments *T-shapes with texture* and *Cardiac MRI* (see Table 4 and Table 7).
For the former, three scales are computed in the multiresolution setting of [8], whereas four
are considered in the latter.

Note that, in the code [8], no boundary conditions are prescribed while in our setting, Dirichlet boundary conditions are enforced. Visually, this may have a slight impact on the reconstructed deformations in particular near the boundaries. Finally, if we consider one step of our algorithm (going from scale $k-1$ to scale $k$) as an execution of the multiresolution algorithm [8], the computation times are comparable. At last, the deformation grids have been designed with the same spacing for both methods.

Three angles of inquiry are addressed:

(i) first of all, of a conceptual nature. As stressed in introduction, multiresolution techniques refer primarily to the size of the images manipulated (which are increasingly large at finer scales) and subsequently to the reduction of the computation cost as well as to the implemented interpolation technique to move from one scale to a finer one. On the contrary, our proposed work focuses more on the image content and on the scale of the features (the size of the images manipulated being the same through the process). It thus seems to us that our approach is more closely related to the information encoded by the images and to the interpretation thereof whereas multiresolution techniques remove arbitrarily some image data. From our point of view, the deformations that result from [8] rather encode the structural main changes and subtle localised displacements cannot be discriminated clearly from one step to another (see in particular the deformation grids of Fig. 4 and Fig. 7). A deformation obtained at a given scale appears rather as an upsampling of the deformation achieved at the previous scale).

(ii) Second, a qualitative/visual comparison of the results produced by both algorithms. In each case, the deformations generated by the algorithms are smooth and the deformed Templates are faithful to reality (quantified in point (iii)). Note that imposing the deformation to be equal to the identity mapping on the boundary in our approach is a strong constraint (but consistent with our theoretical model), which explains the differences in deformation behaviour that can be observed near the boundaries.

Without drawing generalised conclusions, we can nevertheless observe that a slight artefact appears in the deformed Template (scale $k = 3$, bright region of the left side of the right ventricular cavity) in Fig. 7, while our result is closer to reality.

(iii) At last, a quantitative analysis is provided. Again, the figures should be analysed with care as we are not comparing exactly the same thing. Quantitative measures, whether it be Dice coefficient or Re-SSD are slightly better with [8]. Several hypotheses can be put forward in addition to the different boundary conditions. The first is once again linked to the very nature of the multiresolution approach, which is not exactly in line with the philosophy of our approach. While in the last step of the algorithm [8], the exact data (i.e. the original images) are processed, we deal, in our case, with versions of these from which very small details have been removed. Additionally, in method [8], a regridding technique is at work in an un-
derlying way since at a given $k$ step, the initialisation is done with the deformed
Template resulting from the composition of the Template at scale $k$ with the up-
sampled (interpolated) deformation obtained at $k-1$. In our setting however, the
whole composition $\phi_k$ is computed at step $k$ from which we derive $\phi_k$: there is no
regridding involved.

**5. Conclusion.** To conclude, we have introduced a multiscale deformation representa-
tion consisting of the composition of intermediate deformations: the coarser one encodes
the movements of the main structural elements computed from the truncated Template and
Reference, reflecting only the essential features, while the finer one encompasses the local
and refined motion of small items. The proposed variational model relies on hyperelastic-
ity principles to ensure the produced deformations are bi-Lipschitz homeomorphisms and
therefore physically meaningful. Theoretical results including the existence of minimisers
for the model at each scale and an asymptotic result are provided to support the math-
ematical and computational soundness of our approach. Several numerical experiments
are conducted on both synthetic and medical images to show the ability of our model to
produce accurate hierarchical representations of deformations and to deal with noisy data.
A natural extension of our work in 3D is the object of future work, together with an as-
ymptotic analysis of the proposed numerical algorithm, which is a work in progress ([13]).
Replacing the hyperelastic-based regularisation term in our model by other state of the art
regularisation for the registration process could also be interesting to study.

**Acknowledgements.** The authors would like to thank Pr. Modersitzki (Institute of Mathematics and Image Computing, University of Lübeck), Dr. Ruthotto (Emory University, Department of Mathematics), Saskia Neuber (Institute of Mathematics and Image Computing, University of Lübeck) and Pia Franziska Schulz (Institute of Mathematics and Image Computing, University of Lübeck) for providing us with FAIR code associated to [8] and for their valuable advice on parameter estimation.

1086 REFERENCES

1087 [1] L. AMBROSIO, N. FUSCO, AND D. PALLARA, *Functions of Bounded Variation and Free Discontinuity Problems*, Oxford university press, 2000.

1089 [2] K. ASTALA, T. IWANIEC, AND G. MARTIN, *Elliptic Partial Differential Equations and Quasiconformal Mappings in the Plane*, Princeton University Press, 2009.

1091 [3] P. ATHAVALE, R. XU, P. RADAU, A. NACHMAN, AND G. A. WRIGHT, *Multiscale properties of weighted total variation flow with applications to denoising and registration*, Med. Image Anal., 23 (2015), pp. 28–42.

1094 [4] J. M. BALL, *Global invertibility of Sobolev functions and the interpenetration of matter*, P. Roy. Soc. Edin. A, 88 (1981), pp. 315–328.

1096 [5] M. BEG, M. MILLER, A. TROUVÉ, AND L. YOUNES, *Computing large deformation metric mappings via geodesic flows of diffeomorphisms*, Int. J. Comput. Vis., 61 (2005), pp. 139–157.

1099 [6] H. BREZIS, *Analyse fonctionnelle*, Dunod Paris, 2005.

1100 [7] H. BREZIS, *Functional Analysis, Sobolev Spaces and Partial Differential Equations*, Universitext, Springer New York, 2010.

1102 [8] M. BURGER, J. MODERSITZKI, AND L. RUTHOTTO, *A Hyperelastic Regularization Energy for Image Registration*, SIAM J Sci Comput, 35 (2013), p. B132–B148.
1104 [9] P. CIARLET, Three-Dimensional Elasticity, Mathematical Elasticity, Elsevier Science, 1994.

1105 [10] B. DACOROGNA, Direct Methods in the Calculus of Variations, Second Edition, Springer, 2008.

1106 [11] N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCHÖNLIEB, *A Variational Model Dedicated to Joint Segmentation, Registration, and Atlas Generation for Shape Analysis*, SIAM J. Imaging Sci., 13 (2020), pp. 351-380.

1109 [12] N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCHÖNLIEB, *A variational model dedicated to joint segmentation, registration, and atlas generation for shape analysis*, SIAM J. Imaging Sci., 13 (2020), pp. 351-380.

1112 [13] N. DEBROUX, C. LE GUYADER, AND L. VESE, *Asymptotic result for a decoupled nonlinear-elasticity-based registration model*, In preparation, (2022).

1114 [14] N. DEBROUX, C. LE GUYADER, AND L. A. VESE, *Multiscale registration*, in Scale Space and Variational Methods in Computer Vision, A. Elmoataz, J. Fadili, Y. Quéau, J. Rabin, and L. Simon, eds., Cham, 2021, Springer International Publishing, pp. 115-127.

1117 [15] F. DEMENGEL, G. DEMENGEL, AND R. ERNÉ, *Functional Spaces for the Theory of Elliptic Partial Differential Equations*, Universitext, Springer London, 2012.

1119 [16] L. R. DICE, *Measures of the Amount of Ecologic Association Between Species*, Ecology, (1945), pp. 297-302.

1121 [17] B. GRIS, S. DURRLEMAN, AND A. TROUVÉ, *A Sub-Riemannian Modular Framework for Diffeomorphism-Based Analysis of Shape Ensembles*, SIAM J. Imaging Sci., 11 (2018), pp. 802-833.

1124 [18] P. L. LIONS AND B. MERCIER, *Splitting algorithms for the sum of two nonlinear operators*, Siam Journal on Numerical Analysis - SIAM J NUMER ANAL, 16 (1979), pp. 964-979.

1126 [19] K. C. LAM, T. C. NG, AND L. M. LUI, *Multiscale Representation of Deformation via Beltrami Coefficients*, Multiscale Model. Simul., 15 (2017), pp. 864-891.

1128 [20] O. LEHTO, *Quasiconformal Mappings*, Springer New York, 1987.

1129 [21] J. MODERSITZKI, *Numerical Methods for Image Registration*, Oxford University Press, 2004.

1130 [22] J. MODERSITZKI, *FAIR: Flexible Algorithms for Image Registration*, SIAM, 2009.

1131 [23] K. MODIN, A. NACHMAN, AND L. RONDI, *A multiscale theory for image registration and nonlinear inverse problems*, Adv. Math., 346 (2019), pp. 1009-1066.

1133 [24] D. PAQUIN, D. LEVY, E. SCHREIBMANN, AND L. XING, *Multiscale Image Registration*, Math. Biosci. Eng., 3 (2006), pp. 389-418.

1135 [25] D. PAQUIN, D. LEVY, AND L. XING, *Hybrid multiscale landmark and deformable image registration*, Math. Biosci. Eng., 4 (2007), pp. 711-737.

1137 [26] D. PAQUIN, D. LEVY, AND L. XING, *Multiscale deformable registration of noisy medical images*, Math. Biosci. Eng., 5 (2008), pp. 125-144.

1139 [27] L. RISSER, F.-X. VIALARD, R. WOLZ, M. MURGASOVA, D. D. HOLM, AND D. RUECKERT, *Simultaneous Multi-scale Registration Using Large Deformation Diffeomorphic Metric Mapping*, IEEE T. Med. Imaging, 30 (2011), pp. 1746-1759.

1142 [28] S. SOMMER, F. LAUZE, M. NIELSEN, AND X. PENNEC, *Sparse Multi-Scale Diffeomorphic Registration: The Kernel Bundle Framework*, J. Math. Imaging Vis., 46 (2013), pp. 292-308.

1144 [29] A. SOTIRAS, C. DAVATZIKOS, AND N. PARAGIOS, *Deformable medical image registration: A survey*, IEEE Trans. Med. Imaging, 32 (2013), pp. 1153-1190.

1146 [30] D. M. STRONG, J.-F. AUJOL, AND T. F. CHAN, *Scale Recognition, Regularization Parameter Selection, and Meyer's G Norm in Total Variation Regularization*, Multiscale Model Simul., 5 (2006), pp. 273-303.

1149 [31] E. TADMOR, S. NEZZAR, AND L. VESE, *A Multiscale Image Representation Using Hierarchical (BV, L²) Decompositions*, Multiscale Model. Simul., 2 (2004), pp. 554-579.

1151 [32] M. VON SIEBENTHAL, G. SZÉKELY, U. GAMPER, P. BOESIGER, A. LOMAX, AND P. CATTIN, *4D MR imaging of respiratory organ motion and its variability*, Physics in Medicine & Biology, 52 (2007), p. 1547.