# A Variational Model Dedicated to Joint Segmentation, Registration, and Atlas Generation for Shape Analysis

Noémie Debroux, John Aston, Fabien Bonardi, Alistair Forbes, Carole Le Guyader, Marina Romanchikova, Carola-Bibiane Schönlieb

► To cite this version:

Noémie Debroux, John Aston, Fabien Bonardi, Alistair Forbes, Carole Le Guyader, et al. A Variational Model Dedicated to Joint Segmentation, Registration, and Atlas Generation for Shape Analysis. SIAM Journal on Imaging Sciences, 2020, 13 (1), pp.351-380. 10.1137/19M1271907. hal-03034316

HAL Id: hal-03034316

https://uca.hal.science/hal-03034316

Submitted on 1 Dec 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A Variational Model Dedicated to Joint Segmentation, Registration and Atlas Generation for Shape Analysis*

Noémie Debroux†, John Aston‡, Fabien Bonardi§, Alistair Forbes¶, Carole Le Guyader||, Marina Romanchikova||, and Carola-Bibiane Schönlieb†

**Abstract.** In medical image analysis, constructing an atlas, i.e. a mean representative of an ensemble of images, is a critical task for practitioners to estimate variability of shapes inside a population, and to characterise and understand how structural shape changes have an impact on health. This involves identifying significant shape constituents of a set of images, a process called segmentation, and mapping this group of images to an unknown mean image, a task called registration, making a statistical analysis of the image population possible. To achieve this goal, we propose treating these operations jointly to leverage their positive mutual influence, in a hyperelasticity setting, by viewing the shapes to be matched as Ogden materials. The approach is complemented by novel hard constraints on the $L^\infty$ norm of both the Jacobian and its inverse, ensuring that the deformation is a bi-Lipschitz homeomorphism. Segmentation is based on the Potts model, which allows for a partition into more than two regions, i.e. more than one shape. The connection to the registration problem is ensured by the dissimilarity measure that aims to align the segmented shapes. A representation of the deformation field in a linear space equipped with a scalar product is then computed in order to perform a geometry-driven Principal Component Analysis (PCA) and to extract the main modes of variations inside the image population. Theoretical results emphasizing the mathematical soundness of the model are provided, among which existence of minimisers, analysis of a numerical method, asymptotic results and a PCA analysis, as well as numerical simulations demonstrating the ability of the model to produce an atlas exhibiting sharp edges, high contrast and a consistent shape.

**Key words.** Segmentation, registration, nonlinear elasticity, Ogden materials, Potts model, atlas generation, asymptotic results, $D^m$-splines, geometric PCA

**AMS subject classifications.** 68U10, 49, 65D18, 41A15, 62H25.

**1. Introduction.** In recent years, joint image processing models have experienced increasing attention, including combined segmentation/registration models [33, 37] (joint phase field approximation and registration), [52] (model based on metric structure comparison), [29, 68] (level set formulation that merges the piecewise constant Mumford-Shah model with registration principles), [36] (grounded in the expectation maximisation algorithm), [28] (based

*Submitted to the editors DATE.

**Funding:** This project is co-financed by the European Union with the European regional development fund (ERDF, HN0002137) and by the Normandie Regional Council via the M2NUM project.

†Department of Applied Mathematics and Theoretical Physics (DAMPT), Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK (nd448@cam.ac.uk, cbs31@cam.ac.uk).

‡Statslab, Department of Pure Mathematics and Mathematical Statistics, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK (j.aston@statslab.cam.ac.uk).

§IBISC, Université d'Évry, 36, Rue du Pelvoux, CE1455 Courcouronnes, 91020 Evry Cédex, France (fabien.bonardi@univ-evry.fr).

¶National Physical Laboratory, Hampton Road, Teddington, Middlesex, TW11 0LW UK (alistair.forbes, marina.romanchikova@npl.co.uk).

||Normandie Univ, Institut National des Sciences Appliquées de Rouen, Laboratory of Mathematics, 76000 Rouen, France (carole.le-guyader@insa-rouen.fr).
on a nonlocal characterisation of weighted-total variation and nonlocal shape descriptors), or
[2, 48, 59, 62, 70, 76]; joint image reconstruction and motion estimation [11, 16, 22, 58, 64,
69, 15, 53, 8]; joint reconstruction and registration for post-acquisition motion correction [25]
with the goal to reconstruct a single motion-free image and retrieve the physiological dynam-
ics through the deformation maps, joint optical flow estimation with phase field segmentation
of the flow field [14], or joint segmentation/optimal transport models [12] (to determine the
velocity of blood flow in vascular structures). This can be attributed to several factors:

(i) the will to limit error propagation. Indeed, addressing the considered tasks in a unified joint framework (or multitasking) and thus exploiting the strong correlation between them reduces the propagation of uncertainty, contrary to a sequential treatment that may amplify errors from step to step;

(ii) Second, —and this is a corollary of the previous point —, performing simultaneously these tasks yields positive mutual influence and benefit on the obtained results as shown in Figure 1. To exemplify this observation, we can think first of joint models for image reconstruction and registration: not only does the approach correct the misalignment problem, but it also allows for alleviating ghosting artefacts.

In the case of joint segmentation/registration models —the case that will be addressed
more thoroughly afterwards —, as salient component pairing, shape/geometrical fea-
ture matching and intensity distribution comparison drive registration, processing
these tasks simultaneously in a single framework may in particular reduce the influence
of noise since the mapping can be done through the pairing of significant structures,
e.g., by transferring the edges, and not only through intensity correlation.

Besides, the registration operation can be viewed as the inclusion of priors to guide the segmentation process, in particular, for the questions of topology-preservation (the unknown deformation is substituted for the classical evolving curve of the segmentation process —[40, 47][72, Chapter 9] for instance —and the related Jacobian determinant is subject to positivity constraints) and geometric priors (since the registration allows to overcome the issue raised by weak boundary definition due to noise sources in the acquisition device, to degradation of the image contents during reconstruction, etc., by restoring them). In return, relevant segmented structures contribute to fostering accurate registration, providing then a reliable estimation of the deformation between the encoded structures, not only based on intensity matching (taking the form of a local criterion), but also on geometrical/shape pairing (having a nonlocal character).

(iii) Lastly, the pooling of the various results produced by the joint model allows for accu-
rate post-processing treatments based on mutual analysis: for instance, the represen-
tation of the true underlying anatomy of an organ from a set of multiple acquisitions
corrupted by motion, when tackling simultaneously reconstruction and registration, or
the generation of an atlas in the context of joint segmentation and registration. Atlas
refers to a specific model for an ensemble of images and serves as a benchmark, i.e. a
meaningful statistical image, to account for the variability (e.g., different shapes and
sizes for organs in medical imaging) that might be observed in a population.

The proposed work adopts this joint model philosophy. It aims at addressing the issue of
designing a unified variational model for joint segmentation, registration and atlas generation
by exploiting the strong correlation between the two former tasks thus reducing error propaga-
Example of joint segmentation and registration model and the effect on registration (red circles).

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

Template

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

Reference

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

Deformed template
Joint segmentation
and registration [25].

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

Deformed template
Registration alone [28].

Example of our joint model and the impact on the atlas quality (yellow arrows).

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

Segmented atlas via
our joint approach.

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

Segmented atlas via a
sequential approach.

Our joint model creates an atlas with less blurring artefacts and sharper edges than a sequential approach.

Assessement of the segmentation and registration quality :

Joint model : Dice(Reference, deformed Template) = 94.31%,

Registration alone model : Dice(Reference, deformed Template) = 93.71%.

<caption>**Figure 1.** Illustration of the positive impact of joint approaches against sequential ones.</caption>

tion, in the medical imaging setting. The latter one requires the mapping of a group of images to a *mean* representative, which is an additional unknown of the problem, the subsequent goal being to extract a relevant hidden structure from this ensemble of images. As in medical images the variability between individuals is significant, constructing a meaningful statistical image of the global underlying anatomy of an organ from a set of images to measure this variability is of great interest. It allows for the derivation of image statistics, the retrieval of the inherent dynamics of a single individual's organ, the estimation of the probability that a particular spatial location takes on a certain label, the detection and quantisation of abnormalities, that is, more generally, it allows to characterise and understand how geometrical and structural changes influence health. A large body of papers feeds the field of atlas generation and shape statistics among which [39] (atlas generation problem phrased in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [10]), [18] (the shapes to be analysed are modeled as random histograms and in order to learn principal modes of variation from such data, the Wasserstein distance between probability measures is introduced), [78] (dedicated to elastic shape analysis; a unified registration/parameterised object statistical analysis framework is tackled, based on square-root transformations and able to process data as diverse as curves, functions, surfaces and images), [4] (statistics performed on the space of diffeomorphisms), [38] (the use of a kernel descriptor that characterises local shape properties ensures geometrically meaningful correspondence between shapes with statistical studies of the deformations), [62, 63] (the shapes are viewed as closed contours approximated by phase fields, and shape averaging and covariance analysis are carried out in a nonlinear elasticity setting), to name a few.

The difficulty in designing the model arises from the complexity of the formulation that
is generally underconstrained, involves nonlinearity and non-convexity, and is dictated by the given application. While segmentation attempts to reproduce the ability of human beings to track down significant patterns and automatically gather them into significant constituents (see [7, Chapter 4] or [72, Part II] for a relevant analysis of this problem), it remains a challenging and ill-posed task (as emphasised by Zhu et al. ([79])) since the definition of an object encompasses various acceptations: it can be something material —a thing—or a periodic pattern, this heterogeneity entailing the design of suitable methodologies for each specific application. Similarly, for the registration assignment (see [54, 55, 66] for the registration counterpart with Matlab software), the sought deformation is usually viewed as a minimal argument (uniqueness defaults in general) of a specifically tailored cost function that has a polymorphous character in nature. For images acquired on different devices and depicting various physical phenomena, the quality of registration is not measured by intensity distribution alignment, but by the degree of shape/geometrical feature pairing. Also, several stances can be adopted to describe the setting the objects to be matched fall within (physical models —[10], [13], [17], [21], [26], [31], [32], [35], [48], [59], [62] —, purely geometric ones —[5], [27], [65], [77] —, models including a priori knowledge ([24]), depending on the assumption regarding the properties of the deformation to be recovered) and to devise the measure of alignment (that is, how the available data are exploited to drive the registration process), increasing thus the complexity of the problem. To meet these criteria, we devise, in a variational framework, a theoretically well-motivated and physically relevant combined model, capable of handling large deformations, reliable in terms of pairing of the shapes encoded in the images, and efficient in extracting a relevant underlying structure decomposed into shapes from the considered set of images. Statistical shape analysis is then performed by means of a Principal Component Analysis (PCA) on the obtained deformations to retrieve the main modes of variations inside the dataset, after finding a suitable representative of the deformation in a linear space (i.e. in order that the recovered deformation lives in a vector space).

The results are obtained through the use of the hyperelasticity setting and the design of an original geometric dissimilarity measure ensuring alignment of the (possibly nested) shapes for the combined model —thus favouring the matching of shapes rather than the coupling of grey levels with the underlying goal to potentially process images with different modalities—, and the introduction of a tensor-based approximation problem for the statistical analysis. Unlike [62], the shapes to be matched are not modeled by their closed contour but through a piecewise-constant partition (Potts model [56, 61]), which constitutes the main difference with [62]. Not only does the shape pairing rely on the object outer envelope matching, but also on the internal structure alignment. This way of looking at shapes entails substantial modifications in the design of the functional to be minimised and in the search for an appropriate representative of the deformation in a vector space.

More precisely, the novelty of the paper rests upon: (i) an original modeling involving the stored energy function of an Ogden material complemented by two new hard constraints on the Jacobian and its inverse (in addition to the theoretical utility of these constraints, it also allows to control changes of length), the Potts model for segmentation, and an original discrepancy measure ensuring edge mapping; (ii) the derivation of theoretical results encompassing
non straightforward mathematical tools; (iii) the analysis and comparison of three different methods to perform statistical analysis on the obtained deformation: the first one, based on linearised elasticity principles largely inspired by [63], the second one using the Cauchy-stress tensors motivated by [63], and the last one, more novel and on which the paper focuses, relying on tensor-based smoothing $D^m$ splines, influenced by [46]. Section 2 is devoted to the analysis of the mathematical model including a theorem of existence of minimisers, while Section 3 is dedicated to the theoretical analysis of a numerical method of resolution based on a splitting approach and implying the Alternating Direction Method of Multipliers (ADMM) techniques and proximal gradient descent algorithms. Section 4 deals with the resulting geometry-driven statistical analysis, which requires finding a fitting representation of the obtained deformation in a linear space before performing a PCA. As already mentioned, our motivation is to investigate how the linear elasticity based approach and the Cauchy-stress based method compare to the $D^m$ splines approximation based procedure. The first two are non-straightforward adaptations of the techniques envisioned in [63], but the view we take to model the shapes —phase field rather than closed contours —implies substantial changes in the physics of the problem, while the emphasis is put on the last one for its novelty. Section 5 focuses on numerical simulations with a thorough comparison between sequential treatments and the proposed joint model, demonstrating the ability of our model to handle large deformations and to produce in the end, an atlas with sharp edges, high contrast and reflecting a realistic shape.

Let us emphasise that the focus of the paper is on the mathematical presentation of a non-linear elasticity-based unified segmentation, registration, and atlas generation model. Hence, the computational results are currently still restricted to two dimensions due in practice to the applications that were presented to us by clinicians. However, as will be seen next, the proposed algorithm can be easily adapted to the three-dimensional case.

## 2. Mathematical Modeling.

2.1. **Depiction of the Model.** Let $\Omega$ be a connected bounded open subset of $\mathbb{R}^3$ with boundary sufficiently smooth (convenient way of saying that in a given definition, the smoothness of the boundary is such that all arguments make sense and enabling us to use compact Sobolev embeddings among others). Let us denote by $T_i : \bar{\Omega} \to \mathbb{R}$ the $i$-th template image with $i = 1, \dots, M$ —available data in our problem —, $M$ being the total number of initial images. For theoretical and numerical purposes, we assume that each $T_i$ is compactly supported on $\Omega$ to ensure that $T_i \circ \varphi_i$ is always defined and we assume that $T_i$ is Lipschitz continuous. It can thus be considered as an element of the Sobolev space $W^{1,\infty}(\mathbb{R}^3)$, and the chain rule applies. The partitioning of each template $T_i$ into regions with homogeneous intensities, defining shapes, is encoded in the variable $\theta_i : \bar{\Omega} \to \mathbb{R}$ —the variables $\{\theta_i\}_{i=1}^M$ belonging to the set of unknowns of the problem and being read as simplified versions of the images $T_i$ that encompass the geometrical shapes —, and $\theta_R : \bar{\Omega} \to \mathbb{R}$ is the unknown segmented atlas generated by our model. As will be seen later, these variables allow making the connection between segmentation and registration. Also, using these schematic versions of the images tends to favour shape pairing rather than grey level mapping. Let $\varphi_i : \bar{\Omega} \to \mathbb{R}^3$ be the sought deformation between $\theta_{T_i}$ and the unknown mean segmentation $\theta_R$. Of course, in practice, $\varphi_i$ should be with values in $\bar{\Omega}$ but from a mathematical point of view, if we work with such spaces of functions, we lose the structure of vector space. Nonetheless, we can show that our
model retrieves deformations with values in $\bar{\Omega}$ — based on Ball's results [9]. A deformation is a smooth mapping that is orientation-preserving and injective, except possibly on $\partial\Omega$, if self-contact is allowed. The deformation gradient is $\nabla\varphi_i : \bar{\Omega} \to M_3(\mathbb{R})$, the set $M_3(\mathbb{R})$ being the set of real square matrices of order 3. The associated displacement field is denoted by $u_i$ such that $\varphi_i = \text{Id} + u_i$, and $\nabla\varphi_i = I_3 + \nabla u_i$, with $\text{Id}$, the identity mapping, and $I_3$, the $3 \times 3$ identity matrix. We also need the following notations: $A : B = \text{tr}A^T B$, the matrix inner product and $||A|| = \sqrt{A : A}$, the related matrix norm (Frobenius norm).

Following the joint model philosophy in a variational framework, the sought deformations $\varphi_i$, the segmentations $\theta_{T_i}$, for all $i = 1, \dots, M$, as well as the segmented atlas $\theta_R$ are seen as minimal arguments of a specifically designed cost function. It comprises a regularisation on $\varphi_i$, for all $i = 1, \dots, M$, prescribing the nature of the deformations, a penalisation on $\theta_{T_i}$, for all $i = 1, \dots, M$, and $\theta_R$, favouring piecewise constant segmentations, a segmentation fidelity term ensuring the closeness of the $\theta_{T_i}$ to the initial available image $T_i$, and a data-driven term measuring the alignment between the deformed segmentations $\{\theta_{T_i} \circ \varphi_i\}_{i=1}^M$ and $\theta_R$, intervening then segmentation and registration.

In this work, we view all the template images and their respective partitioning as deformed versions of a single image/segmentation. Inspired by the observation in [62]: "the arithmetic mean $x$ of observations $\{x_i\}_{i=1}^M$ can be interpreted as the minimiser of the total elastic deformation energy in a system where the average $x$ is connected to each $x_i$ by an elastic spring, under the Hooke's law", a natural choice for the definition of the mean segmentation is given by the particular deformed configuration that minimises the total nonlinear hyperelastic deformation energy required to align each segmentation to this mean configuration. To allow large deformations, the shapes to be matched are viewed as isotropic (exhibiting the same mechanical properties in every direction), homogeneous (showing the same behaviour everywhere inside the material), and hyperelastic (with a stress-strain relation derived from a strain energy density) materials, and more precisely as Ogden ones ([23]). Note that rubber, filled elastomers, and biological tissues are often modeled within the hyperelastic framework, which motivates our modeling. This perspective drives the design of the regularisation on the deformations $\varphi_i$ which is thus based on the stored energy function of an Ogden material, prescribing then a physically-meaningful nature.

We recall that the general expression for the stored energy function of an Ogden material (see [23, 45]) is given by

$$W_O(F) = \sum_{i=1}^{K_1} a_i ||F||^{\gamma_i} + \sum_{j=1}^{K_2} b_j ||\text{Cof}F||^{\beta_j} + \Gamma(\det F),$$

with $a_i > 0$, $b_j > 0$, $\gamma_i \ge 0$, $\beta_j \ge 0$, for all $i = 1, \dots, K_1$ and all $j = 1, \dots, K_2$, and $\Gamma: ]0, +\infty[ \to \mathbb{R}$ being a convex function satisfying $\lim_{\delta \to 0^+} \Gamma(\delta) = \lim_{\delta \to +\infty} \Gamma(\delta) = +\infty$. The first term penalises changes in length, the second one controls the changes in area while the third one restricts the changes in volume. The latter also ensures preservation of topology by imposing positivity of the Jacobian determinant almost everywhere. In this work, we focus
on the following particular energy:

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

fulfilling the previous assumptions. The third and fourth terms govern the distribution of the Jacobian determinant: the latter prevents singularities and large contractions by penalising small values of the determinant, while the former promotes values of the determinant close to 1 avoiding thus expansions and contractions that are too large. The choice of the remaining terms is motivated by the theoretical results in [9] to ensure that the deformations are homeomorphisms. The constants are added to fulfill the energy property $W_{Op}(I_3) = 0$. In order to avoid singularity as much as possible, to get deformations that are bi-Lipschitz homeomorphisms, and to obtain Cauchy-stress tensors (whose formal definition will be given in Section SM4) in the linear space $L^2(\Omega, M_3(\mathbb{R}))$, we complement this stored energy function $W_{Op}$ by the term $1_{\{\|\cdot\|_{L^\infty(\Omega, M_3(\mathbb{R}))} \le \alpha\}}(F) + 1_{\{\|\cdot\|_{L^\infty(\Omega, M_3(\mathbb{R}))} \le \beta\}}(F^{-1})$, with $\alpha \ge 1$, and $\beta \ge 1$, where $1_A$ denotes the convex characteristic function of a convex set $A$. Therefore, the regularisation can be written as

$$W(F) = \int_{\Omega} W_{Op}(F) dx + 1_{\{\|\cdot\|_{L^{\infty}(\Omega, M_3(\mathbb{R}))} \le \alpha\}}(F) + 1_{\{\|\cdot\|_{L^{\infty}(\Omega, M_3(\mathbb{R}))} \le \beta\}}(F^{-1}).$$

**Remark 2.1.** In terms of functional spaces, if $\varphi \in W^{1,\infty}(\Omega, \mathbb{R}^3)$ (suitable space owing to the $L^\infty$ hard constraints), Cof $\nabla\varphi$ and $\det\nabla\varphi$ are automatically elements of $L^\infty(\Omega, M_3(\mathbb{R}))$ and $L^\infty(\Omega)$ respectively, since $L^\infty(\Omega, \mathbb{R}^3)$ has a structure of commutative Banach algebra. Penalising the $L^\infty$ norm of $\nabla\varphi$ thus entails control over the Jacobian determinant. This additional term implicitly gives an upper and lower bound on the Jacobian determinant ensuring thus topology preservation.

**Remark 2.2.** In an effort to ensure orientation preservation, a new framework based on quasiconformal mappings has been introduced, as developed in [43, 44, 49, 20]. This approach takes root in mathematical considerations (the deformation is ruled by mathematical criteria —it satisfies the Beltrami equation —rather than physical ones) and relies on the fine properties of quasiconformal mappings. Quasiconformal mappings can be defined as follows ([50]) (we restrict ourselves to quasiconformal mappings that are homeomorphisms between plane domains): ‘A sense-preserving homeomorphism $f$ of the domain $G$ is called quasiconformal if its maximal dilation $K(G)$ is finite. If $K(G) \le K < \infty$, then $f$ will be called $K$-quasiconformal.’ As from the one hand, $K(G) \ge 1$ for non-conformal sense-preserving homeomorphisms, and from the other hand, if $w$ is conformal, $K(G) = 1$, $K(G)$ can be viewed as a measure of deviation from conformality. The following result gives a necessary and sufficient condition ensuring that a given homeomorphism $f: \Omega \to \Omega'$ in $W_{loc}^{1,2}(\Omega)$ is $K$-quasiconformal. This characterization serves as the core of the models [43, 44, 49, 20] which aim to control the Beltrami coefficient.
**Theorem 2.3.** [6, Theorem 2.5.4] Suppose $f: \Omega \to \Omega'$ is a homeomorphic $W_{loc}^{1,2}$-mapping. Then $f$ is $K$-quasiconformal if and only if

$$ \frac{\partial f}{\partial \bar{z}}(z) = \mu(z) \frac{\partial f}{\partial z}(z) \text{ for almost everywhere } z \in \Omega, $$

where $\mu$, called the Beltrami coefficient of $f$ is a bounded measurable function satisfying

$$ \|\mu\|_{\infty} \leq \frac{K-1}{K+1} < 1 $$

With $\frac{\partial f}{\partial z} = \frac{1}{2}(\frac{\partial f}{\partial x} - i \frac{\partial f}{\partial y})$ and $\frac{\partial f}{\partial \bar{z}} = \frac{1}{2}(\frac{\partial f}{\partial x} + i \frac{\partial f}{\partial y})$, setting $f(z) = f(x+iy) = u(x,y)+iv(x,y)$ (so in the context of registration, the sought deformation is $\varphi = (u,v)^T$), the Jacobian is defined by $J_f(z) = \left|\frac{\partial f}{\partial z}(z)\right|^2 - \left|\frac{\partial f}{\partial \bar{z}}(z)\right|^2 = \frac{\partial u}{\partial x}(x,y) \frac{\partial v}{\partial y}(x,y) - \frac{\partial u}{\partial y}(x,y) \frac{\partial v}{\partial x}(x,y) = \det \nabla \varphi(x,y)$. Consequently, if $f, W_{loc}^{1,2}$-homeomorphism, is $K$-quasiconformal,

$$ \begin{aligned} J_f(z) &= \left|\frac{\partial f}{\partial z}(z)\right|^2 - \left|\frac{\partial f}{\partial \bar{z}}(z)\right|^2, \\ &= \left|\frac{\partial f}{\partial z}(z)\right|^2 (1 - |\mu(z)|^2), \end{aligned} $$

entailing that $J_f(z) = \det \nabla \varphi(x,y) > 0$ a.e., since $\|\mu\|_{\infty} < 1$. This observation allows to make a connection between our model and the quasiconformal setting in the sense that ensuring boundedness of the Beltrami coefficient ($\|\mu\|_{\infty} < 1$) implies positivity of the related deformation Jacobian determinant, which is the constraint we aim to fulfill.

The aforementioned regulariser is then applied along with a discrepancy measure, which allows intertwining the segmentation and registration tasks, and a segmentation part comprising a fidelity term and a sparsity measure on the paired edges based on the Potts model ([61]). The latter, also known as piecewise-constant Mumford-Shah model [56] with $N \in \mathbb{N}$ phases/shapes ($N$ is thus a prior), is written, for an observed image $f$, as

$$ \inf_{u \in \mathcal{U}} E_{Potts}(u) = \sum_{l=1}^{N} \frac{\alpha}{2} TV(u_l) + \int_{\Omega} \sum_{l=1}^{N} u_l (c_l - f)^2 dx, $$

with $\mathcal{U} = \{u = (u_l)_{l=1,...,N} \in (BV(\Omega, \{0,1\}))^N, \sum_{l=1}^{N} u_l = 1, \text{ a.e. on } \Omega\}$, and

$c_l = \begin{cases} \frac{\int_{\Omega} f u_l \, dx}{\int_{\Omega} u_l \, dx} & \text{if } \int_{\Omega} u_l \, dx \neq 0, \\ 0 & \text{otherwise,} \end{cases}$ , $\alpha > 0$ being a weighting parameter balancing the fidelity term and the regularisation. The notation $TV$ denotes the classical Total Variation, measuring the perimeter length of the set defined by $\{x \in \Omega | u_l(x) = 1\}$ thanks to the coarea formula ([30, 34]). The segmentation/partitioning is then retrieved by $\tilde{u} = \sum_{l=1}^{N} u_l c_l$, $\tilde{u}$ being a decomposition
of the initial image $f$ into $N$ shapes defined by the characteristic functions $u_l$ with constant intensity values $c_l$, each one corresponding to an object of interest under the assumption that it is defined by a homogeneous region with close intensity values.

**Remark 2.4.** Extensions to homogeneous regions in terms of texture with a piecewise-smooth approximation instead of a piecewise-constant approximation (see [72]) or in terms of histograms (see [60]) are possible, depending on the nature of the considered images, but this is not the scope of this paper.

The characteristic functions $u_l$ give a good representation of the geometric features inside the images, and can be seen as nonlocal shape descriptors that will help the registration process. In that prospect, we introduce this novel geometric dissimilarity measure whose aim is to align the salient structures based on the previous decomposition without taking into account the intensity values —thus favouring shape pairing —:

$$E_{\text{dist}}((\theta_{T_i}, \varphi_i)_{i=1, \dots, M}, \theta_R) = \frac{1}{2M} \sum_{i=1}^{M} \sum_{l=1}^{N} TV(\theta_{T_i, l} \circ \varphi_i - \theta_{R, l}),$$

with notations consistent with the definition of $\mathcal{U}$, i.e. $\forall i \in \{1, \dots, M\}, \theta_{T_i} = (\theta_{T_i, l})_{l=1, \dots, N} \in (BV(\Omega, \{0, 1\}))^N$ and $\theta_R = (\theta_{R, l})_{l=1, \dots, N} \in (BV(\Omega, \{0, 1\}))^N$.

**Remark 2.5.** Consistently with Remark 2.4, we could also envision a model including both the deformations $\varphi_i$ pairing the structures (i.e. viewed as global deformations) and additional components reflecting better the more local deformations. This results mathematically in a composition of deformations as in [71] in which they design a topology-preserving registration model based on diffusion principles and demons' algorithm consisting of composing small deformations to obtain the final global one. Again, this is not the scope of the proposed work.

It thus allows for the registration of images acquired through different mechanisms and is more robust to small changes of intensities that can happen even for images of the same modality, especially in medical images. It measures the perimeter length of the misaligned region for each structure of interest and thus drives the registration process by mapping the shapes.

In the end, the overall problem denoted by (P) is stated by

$$\begin{array}{l}
\inf \mathcal{F}_1(\theta_R, \{\theta_{T_i}, \varphi_i\}_{i=1}^M) = \frac{1}{M} \sum_{i=1}^M \left( \frac{\gamma_T}{2} \sum_{l=1}^N TV(\theta_{T_i, l}) + \int_{\Omega} \sum_{l=1}^N \theta_{T_i, l} (c_{T_i, l} - T_i)^2 dx \right. \\
\qquad \left. + \frac{\gamma_R}{2} \sum_{l=1}^N TV(\theta_{R, l}) + \int_{\Omega} \sum_{l=1}^N \theta_{R, l} (c_{R, l} - T_i \circ \varphi_i)^2 dx + \frac{\lambda}{2} \sum_{l=1}^N TV(\theta_{T_i, l} \circ \varphi_i - \theta_{R, l}) + W(\nabla \varphi_i) \right), \\
\text{with } c_{T_i, l} = \begin{cases} \frac{\int_{\Omega} \theta_{T_i, l}(x) T_i(x) dx}{\int_{\Omega} \theta_{T_i, l}(x) dx} & \text{if } \int_{\Omega} \theta_{T_i, l}(x) dx \neq 0, \\ 0 & \text{otherwise} \end{cases}
\end{array}$$
$$c_{R,l} = \begin{cases} \frac{1}{M} \sum_{i=1}^{M} \frac{\int_{\Omega} \theta_{R,l}(x) T_i \circ \varphi_i(x) dx}{\int_{\Omega} \theta_{R,l}(x) dx} & \text{if } \int_{\Omega} \theta_{R,l}(x) dx \neq 0, \alpha \ge 1 \text{ and } \beta \ge 1. \\ 0 & \text{otherwise} \end{cases}$$

An illustration of the overall components of the model as well as the pipeline of the resulting analysis is given in Figure 2.

<figure><img src="image_8.png" /><figcaption>Figure 2. Overview of our framework</figcaption></figure>

## 2.2. Theoretical Results.
In this subsection, we theoretically analyse problem (P) by showing its well-definedness. In that purpose, we prove the existence of minimisers in the following theorem.

**Theorem 2.6 (Existence of minimisers).**

We introduce the functional space:

$$\begin{gathered} \hat{W} = \{ \psi \in \mathrm{Id} + W_0^{1,\infty}(\Omega, \mathbb{R}^3), \frac{1}{\det \nabla \psi} \in L^{10}(\Omega), \det \nabla \psi > 0 \text{ a.e. in } \Omega, \\ \|\nabla \psi\|_{L^{\infty}(\Omega, M_3(\mathbb{R}))} \le \alpha, \|(\nabla \psi)^{-1}\|_{L^{\infty}(\Omega, M_3(\mathbb{R}))} \le \beta \}, \end{gathered}$$

The infimum is searched for $\theta_R \in \mathcal{U}$, $\theta_{T_i} \in \mathcal{U}$, and $\varphi_i \in \hat{W}$ for all $i \in \{1, \dots, M\}$. There exists at least one minimiser to this problem.

**Remark 2.7.** We would like to highlight that the design of the considered functional spaces are dictated by mathematical and practical reasons, and that we cannot fundamentally relax the constraints in their formulation.

(i) As $\forall i \in \{1, \dots, M\}$, $\theta_{T_i}$ is an unknown of the problem, when addressing sub-problem 1 in the algorithm —optimisation over $\theta_{T_i}$ —, we need to rephrase explicitly the objective function in terms of $\varphi_i^{-1}$ in order to make the computations feasible. It implies that from a mathematical viewpoint, the functional space the $\varphi_i$'s belong to should ensure
that the $\varphi_i^{-1}$'s exist and that the deformations are with values in $\bar{\Omega}$, which is guaranteed here with the devised stored energy function $W_{Op}(F)$ and thanks to Ball's results. In particular, $\varphi_i$ is a bi-Lipschitz homeomorphism such that $\det \nabla \varphi_i > 0$ a.e.

(ii) According to [1, Theorem 2.6.], if $F: U \to V$ is a bi-Lipschitz homeomorphism between open subsets of $\mathbb{R}^N$ such that, either $\det Df \ge 0$ $\mathcal{H}^N$-a.e. on $U$ or $\det Df \le 0$ $\mathcal{H}^N$-a.e. on $U$, $\forall f \in BV(U)$, the function $f \circ F^{-1}$ belongs to $BV(V)$. Our framework fulfills these criteria, so $\theta_{T_i,l} \circ \varphi_i \in BV(\Omega)$ and the vector space structure of $BV$ gives us immediately that for all $i \in \{1, \dots, M\}$, $\theta_{T_i,l} \circ \varphi_i - \theta_{R,l} \in BV(\Omega)$.

*Proof.* The proof is based on the theory of the calculus of variations, and relies on Ball's results [9] and arguments inspired by [75]. See section SM1 for the detailed proof. ■

We now investigate an original numerical method for the resolution of (P).

### 3. Numerical Method of Resolution.

**3.1. Description and Analysis of the Numerical Method.** Inspired by a prior work by Negrón Marrero [57] in which the author describes and analyses a numerical method detecting singular minimisers and avoiding the Lavrentiev phenomenon for 3D problems in nonlinear elasticity, we introduce auxiliary variables and split the original problem into sub-problems that are computationally more tractable. The idea of Marrero's work is to decouple the deformation $\varphi$ from its gradient $\nabla\varphi$ and to formulate a related decoupled problem under equality constraints, moving thus the nonlinearity in the Jacobian to this new variable. With this in mind, we introduce the following auxiliary variables: $V_i$ simulating the Jacobian of $\varphi_i$ for each $i$, $W_i$ simulating the inverse Jacobian $(\nabla\varphi_i)^{-1}$ for all $i$, and $\forall i \in \{1, \dots, M\}$, $\forall l \in \{1, \dots, N\}$, $\theta_{\tilde{T}_i,l} = \theta_{T_i,l} \circ \varphi_i - \theta_{R,l}$, to simplify numerical computations, and derive a functional minimisation problem phrased in terms of $(\theta_{T_i}, \varphi_i, V_i, W_i, \theta_{\tilde{T}_i})_{i=1,\dots,M}, \theta_{R}$. However, we do not impose equality constraints as in [57], but integrate instead, $L^p$-type penalisations ($p = 1$ or $p = 2$; the choice for the $L^1$-penalisation will be discussed later) into the functional, partially relaxing a constrained problem under both equality and inequality constraints by a problem under inequality constraints only. The decoupled problem is thus defined as follows:

$$ \inf \left\{ \begin{aligned} & \mathcal{F}_{1,\gamma}(\{\varphi_i, \theta_{T_i}, V_i, W_i\}_{i=1}^M, (\theta_{\tilde{T}_i,l})_{\substack{i=1,\dots,M \\ l=1,\dots,N}}, \theta_R) \\ & \quad + \frac{1}{M} \sum_{i=1}^M \left( \frac{\gamma_T}{2} \sum_{l=1}^N TV(\theta_{T_i,l}) \right. \\ & \quad \qquad \left. + \int_{\Omega} \sum_{l=1}^N \theta_{T_i,l} (c_{T_i,l} - T_i)^2 dx + \frac{\gamma_R}{2} \sum_{l=1}^N TV(\theta_{R,l}) + \int_{\Omega} \sum_{l=1}^N \theta_{R,l} (c_{R,l} - T_i \circ \varphi_i)^2 dx \right. \\ & \quad \qquad \left. + \frac{\lambda}{2} \sum_{l=1}^N TV(\theta_{\tilde{T}_i,l}) + \gamma \int_{\Omega} \sum_{l=1}^N |\theta_{\tilde{T}_i,l} - (\theta_{T_i,l} \circ \varphi_i - \theta_{R,l})| dx + \int_{\Omega} W_{Op}(V_i) dx \right. \\ & \quad \qquad \left. + \frac{\gamma}{4} \|V_i - \nabla \varphi_i\|_{L^4(\Omega, M_3(\mathbb{R}))}^4 + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_3(\mathbb{R}))} \le \alpha\}}(V_i) + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_3(\mathbb{R}))} \le \beta\}}(W_i) \right\} \\ & \quad + \frac{\gamma}{2} \|W_i - V_i^{-1}\|_{L^2(\Omega, M_3(\mathbb{R}))}^2 \end{aligned} \right\} \quad (\text{DP}) $$
with $\forall l \in \{1, \dots, N\}, c_{T_i,l} = \begin{cases} \frac{\int_{\Omega} \theta_{T_i,l}(x) T_i(x) dx}{\int_{\Omega} \theta_{T_i,l}(x) dx} & \text{if } \int_{\Omega} \theta_{T_i,l}(x) dx \neq 0, \\ 0 & \text{otherwise} \end{cases}$,

$c_{R,l} = \begin{cases} \frac{1}{M} \sum_{i=1}^{M} \frac{\int_{\Omega} \theta_{R,l}(x) T_i(x) \varphi_i(x) dx}{\int_{\Omega} \theta_{R,l}(x) dx} & \text{if } \int_{\Omega} \theta_{R,l}(x) dx \neq 0, \alpha \ge 1 \text{ and } \beta \ge 1. \text{ We address this} \\ 0 & \text{otherwise} \end{cases}$

problem for $\varphi_i \in \text{Id} + W_0^{1,4}(\Omega, \mathbb{R}^3)$, $V_i \in \{\xi \in L^\infty(\Omega, M_3(\mathbb{R})) \mid \det\xi > 0 \text{ a.e. on } \Omega, \frac{1}{\det\xi} \in L^{10}(\Omega), \|\xi\|_{L^\infty(\Omega, M_3(\mathbb{R}))} \le \alpha\}$, $W_i \in \{\xi \in L^2(\Omega, M_3(\mathbb{R})) \mid \|\xi\|_{L^\infty(\Omega, M_3(\mathbb{R}))} \le \beta\}$, $\theta_{T_i} \in \mathcal{U}$ such that $\theta_{T_i,l} \circ \varphi_i \in L^1(\Omega)$ for all $l \in \{1, \dots, N\}$ and for all $i \in \{1, \dots, M\}$, $\theta_{\bar{T}_i,l} \in BV(\Omega, \{-1, 0, 1\})$ for all $l \in \{1, \dots, N\}$ and for all $i \in \{1, \dots, M\}$, and $\theta_R \in \mathcal{U}$.

**3.2. Theoretical Results.** In this subsection, we theoretically analyse problem (DP) and show an asymptotic result relating the decoupled problem (DP) to the initial problem (P).

**Theorem 3.1 (Asymptotic result).** Let $(\gamma_j)_{j \ge 0}$ be an increasing sequence of positive real numbers such that $\lim_{j \to +\infty} \gamma_j = +\infty$. Let $(\{\varphi_{i,k_j}, \theta_{T_i,k_j}, V_{i,k_j}, W_{i,k_j}\}_{i=1}^M, (\theta_{\bar{T}_i,l,k_j})_{\substack{i=1,\dots,M \\ l=1,\dots,N}})$ be a minimising sequence of the problem $\mathcal{F}_{1,\gamma}$ for $\gamma = \gamma_j$. Then there exists a subsequence such that $\varphi_{i,k_j} \stackrel{W^{1,4}(\Omega,\mathbb{R}^3)}{\underset{j \to +\infty}{\longrightarrow}} \bar{\varphi}_i, \theta_{T_i,k_j} \stackrel{(L^1(\Omega))^N}{\underset{j \to +\infty}{\longrightarrow}} \bar{\theta}_{T_i}, \theta_{R,k_j} \stackrel{(L^1(\Omega))^N}{\underset{j \to +\infty}{\longrightarrow}} \bar{\theta}_R, V_{i,k_j} \stackrel{*}{\underset{j \to +\infty}{\longrightarrow}} \nabla \bar{\varphi}_i \text{ in } L^\infty(\Omega, M_3(\mathbb{R}))$, $W_{i,k_j} \stackrel{*}{\underset{j \to +\infty}{\longrightarrow}} (\nabla \bar{\varphi}_i)^{-1} \text{ in } L^\infty(\Omega, M_3(\mathbb{R}))$, $\theta_{\bar{T}_i,l,k_j} \stackrel{L^1(\Omega)}{\underset{j \to +\infty}{\longrightarrow}} \bar{\theta}_{T_i,l} \circ \bar{\varphi}_i - \bar{\theta}_{R,l}$, for all $l = 1, \dots, N$, for all $i = 1, \dots, M$, and $\lim_{j \to +\infty} \mathcal{F}_{1,\gamma_j}(\{\varphi_{i,k_j}, \theta_{T_i,k_j}, V_{i,k_j}, W_{i,k_j}\}_{i=1}^M, (\theta_{\bar{T}_i,l,k_j})_{\substack{i=1,\dots,M \\ l=1,\dots,N}}) = \mathcal{F}_1(\bar{\theta}_R, \{\bar{\theta}_{T_i}, \bar{\varphi}_i\}_{i=1}^M) = \inf \mathcal{F}_1$, so that $(\bar{\theta}_R, \{\bar{\theta}_{T_i}, \bar{\varphi}_i\}_{i=1}^M) \in \mathcal{U}^{M+1} \times \hat{\mathcal{W}}^M$ is a minimiser of the initial problem (P).

*Proof.* This proof is divided into three parts. The first one consists of deriving a coercivity inequality. The second one shows the convergence of a minimising sequence and the last one is dedicated to the lower semi-continuity of the functional. See section *section SM2* for a detailed proof. ■

Equipped with this material and argument, we propose the following discretised numerical scheme.

**3.3. Numerical Scheme.** In this subsection, we restrict ourselves to the two-dimensional case and make some minor changes to the model for the purpose of simplicity but the extension to the three-dimensional case shouldn't induce additional challenges.

We now consider the following discrete two-dimensional decoupled problem —note that from now on, $\theta_{T_i}$ denotes the partition of $T_i$ into piecewise constant regions, i.e. $\theta_{T_i} =$
$\sum_{l=1}^{N} c_{T_i,l} \theta_{T_i,l}$, the number of shapes, $N$ being an unknown (see Remark 3.3) similarly for $\theta_R$—

$$ \inf \left\{ \begin{array}{l} \mathcal{F}_{2,\gamma}(\{\varphi_i, \theta_{T_i}, V_i, \theta_{\tilde{T}_i}, W_i\}_{i=1}^M, \theta_R) = \frac{1}{M} \sum_{i=1}^M \gamma_T \|\nabla \theta_{T_i}\|_{L^0(\Omega)} + \gamma_{\tilde{T}} \|\nabla \theta_{\tilde{T}_i}\|_{L^0(\Omega)} \\ \qquad + \lambda_T \|\theta_{T_i} - T_i\|_{L^2(\Omega)}^2 + \gamma_R \|\nabla \theta_R\|_{L^0(\Omega)} + \lambda_R \|\theta_R - T_i \circ \varphi_i\|_{L^2(\Omega)}^2 \\ \qquad + \frac{\gamma_2}{2} \|V_i - \nabla \varphi_i\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 + \frac{\gamma_1}{2} \|\theta_{\tilde{T}_i} - (\theta_{T_i} \circ \varphi_i - \theta_R)\|_{L^2(\Omega)}^2 + \int_{\Omega} W'_{Op}(V_i, \det V_i) dx \\ \\ (\text{DPb}) \quad + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(V_i) + \mathbb{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}}(W_i) + \frac{\gamma_3}{2} \|W_i - V_i^{-1}\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 \end{array} \right\} $$

$$ \text{with } W'_{Op}(\psi, \delta) = \begin{cases} a_1 \|\psi\|^4 + a_2(\delta - 1)^2 + \frac{a_3}{\delta^{10}} - 2a_1 - a_3 & \text{if } \delta > 0 \\ +\infty & \text{otherwise} \end{cases}. $$

**Remark 3.2.** In the two-dimensional case, the cofactor matrix vanishes and we only need an $L^2$-penalisation to get the asymptotic result as in [28]. Also $BV(\Omega) \hookrightarrow L^2(\Omega)$ in 2 dimensions, so we can replace the $L^1$-penalisation for the auxiliary variable $\theta_{\tilde{T}_i}$ by an $L^2$-penalisation term.

**Remark 3.3.** We have also opted for the discrete Potts model for the segmentation as in [67] since it does not require any prior knowledge on the number of shapes in the image. If the number of shapes is known a priori, another approach based on convexification as in [19] can be applied.

We address this optimisation problem by an alternating scheme in which we fix all the variables except one and solve the subproblem related to the remaining unknown iteratively.

* **Sub-problem 1. Optimisation over $\theta_{T_i}$**. For each $i = 1, \dots, M$, the problem in $\theta_{T_i}$ amounts to solve

$$ \begin{align*} & \inf_{\theta_{T_i}} \gamma_T \|\nabla \theta_{T_i}\|_{L^0(\Omega)} + \lambda_T \|\theta_{T_i} - T_i\|_{L^2(\Omega)}^2 + \frac{\gamma_1}{2} \|\theta_{\tilde{T}_i} - \theta_{T_i} \circ \varphi_i + \theta_R\|_{L^2(\Omega)}^2, \\ \Leftrightarrow & \inf_{\theta_{T_i}} \gamma_T \|\nabla \theta_{T_i}\|_{L^0(\Omega)} + \lambda_T \|\theta_{T_i} - T_i\|_{L^2(\Omega)}^2 \\ & \qquad + \frac{\gamma_1}{2} \|(\theta_{\tilde{T}_i} \circ \varphi_i^{-1} - \theta_{T_i} + \theta_R \circ \varphi_i^{-1})(\det \nabla \varphi_i)^{-\frac{1}{2}}\|_{L^2(\Omega)}^2, \\ \Leftrightarrow & \inf_{\theta_{T_i}} \gamma_T \|\nabla \theta_{T_i}\|_{L^0(\Omega)} \\ & \qquad + \left\| \sqrt{\lambda_T + (\det \nabla \varphi_i)^{-1} \frac{\gamma_1}{2}} \theta_{T_i} - \frac{\lambda_T T_i + \frac{\gamma_1}{2} (\det \nabla \varphi_i)^{-1} (\theta_{\tilde{T}_i} \circ \varphi_i^{-1} + \theta_R \circ \varphi_i^{-1})}{\sqrt{\lambda_T + (\det \nabla \varphi_i)^{-1} \frac{\gamma_1}{2}}} \right\|_{L^2(\Omega)}^2. \end{align*} $$

This amounts to solve the Potts model with nonnegative weights and we use the algorithm in [67] based on the Alternating Direction Method of Multipliers (ADMM) and linear programming.
* **Sub-problem 2. Optimisation over $\theta_{\bar{T}_i}$.** For each $i = 1, \dots, M$, the sub-problem in $\theta_{\bar{T}_i}$ is the following one:

$$ \inf_{\theta_{\bar{T}_i}} \gamma_{\bar{T}} \|\nabla \theta_{\bar{T}_i}\|_{L^0(\Omega)} + \frac{\gamma_1}{2} \|\theta_{\bar{T}_i} - (\theta_{T_i} \circ \varphi_i - \theta_R)\|_{L^2(\Omega)}^2. $$

This is again the Potts model and we use the same algorithm [67] to solve it in practice.

* **Sub-problem 3. Optimisation over $\theta_R$.** By fixing all the other variables, the optimisation problem with respect to $\theta_R$ becomes

$$ \inf_{\theta_R} \gamma_R \|\nabla \theta_R\|_{L^0(\Omega)} + \frac{1}{M} \sum_{i=1}^{M} \lambda_R \|\theta_R - T_i \circ \varphi_i\|_{L^2(\Omega)} + \frac{\gamma_1}{2} \|\theta_R - (\theta_{T_i} \circ \varphi_i - \theta_{\bar{T}_i})\|_{L^2(\Omega)}^2, $$

$$ \Leftrightarrow \inf_{\theta_R} \gamma_R \|\nabla \theta_R\|_{L^0(\Omega)} + \|(\lambda_R + \frac{\gamma_1}{2})(\theta_R - \frac{(\frac{1}{M}\sum_{i=1}^M \lambda_R T_i \circ \varphi_i + \frac{\gamma_1}{2}(\theta_{T_i} \circ \varphi_i) - \theta_{\bar{T}_i})}{\lambda_R + \frac{\gamma_1}{2}})\|_{L^2(\Omega)}^2. $$

This is again a Potts model that we solve with the algorithm [67].

* **Sub-problem 4. Optimisation over $V_i$.** For each $i = 1, \dots, M$, the sub-problem in $V_i$ reads

$$ \inf_{V_i} F(V_i) + Reg(V_i) = \int_{\Omega} a_1 \|V_i\|^4 + a_2(\det V_i - 1)^2 + \frac{a_3}{(\det V_i)^{10}} dx \\ + \frac{\gamma_2}{2} \|V_i - \nabla \varphi_i\|_{L^2(\Omega, M_2(\mathbb{R}))}^2 + \frac{\gamma_3}{2} \|W_i - V_i^{-1}\|_{L^2(\Omega)}^2 + \mathbf{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(V_i). $$

This can be cast as a structured convex non-smooth optimisation problem of the sum of a proper closed convex function $Reg(\cdot) = \mathbf{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(\cdot)$ and a smooth function $F$ corresponding to the remaining of the functional. This is a classical optimisation problem and several schemes have been developed to solve it. In practice, we use the simple iterative forward-backward splitting algorithm [41]:

$$ V_i^{k+1} = \text{prox}_{\gamma Reg}(V_i^k - \gamma \nabla F(V_i^k)), $$

with $\text{prox}_{\gamma Reg}(y) = \min_x \frac{1}{2}\|x-y\|_2^2 + \gamma Reg(y) = \min_x \frac{1}{2}\|x-y\|_2^2 + \gamma \mathbf{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(y) = P_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(y), P_C$ being the projection operator onto the convex set $C$. This could be improved in future work by using for instance the algorithm proposed in [51].

* **Sub-problem 5. Optimisation over $W_i$.** For each $i = 1, \dots, M$, we solve the following minimisation problem

$$ \inf_{W_i} \frac{\gamma_3}{2} \|W_i - V_i^{-1}\|_{L^2(\Omega)}^2 + \mathbf{1}_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \beta\}}(W_i) = P_{\{\|\cdot\|_{L^\infty(\Omega, M_2(\mathbb{R}))} \le \alpha\}}(V_i^{-1}). $$

* **Sub-problem 6. Optimisation over $\varphi_i$.** For each $i = 1, \dots, M$, the sub-problem in $\varphi_i$ reads

$$ \inf_{\varphi_i} \gamma_R \|\theta_R - T_i \circ \varphi_i\|_{L^2(\Omega)}^2 + \frac{\gamma_1}{2} \|\theta_{\bar{T}_i} - \theta_{T_i} \circ \varphi_i + \theta_R\|_{L^2(\Omega)}^2 + \frac{\gamma_2}{2} \|V_i - \nabla \varphi_i\|_{L^2(\Omega)}^2. $$

We propose to solve the associated Euler-Lagrange equation using an $L^2$-gradient flow scheme with an implicit Euler time stepping.
The overall algorithm is summarised in section SM3.

*Remark 3.4.* Under mild assumptions —by replacing the $L^0$-penalization by an $L^1$ one —, we can prove the convergence of the algorithm as well as an asymptotic result.

We now turn to the geometry-driven statistical analysis.

**4. Representation of the Deformations in a Linear Space and Geometry-driven PCA.**

In this section, we focus on the performance of 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. The main hindrance is that our deformation maps live in a nonlinear space whereas classical statistical tools require the objects to be in a linear space. Therefore, we first need to find a good representation of our deformations in a linear space equipped with a scalar product (in order to compute the covariance operator), enabling us to perform a Principal Component Analysis (PCA) on these representatives afterwards. The fundamental axiom of elasticity stating that the energy required to deform an object from a state of reference to another equilibrium state is the same whatever the chosen path is, prevents a straightforward definition of geodesics. Therefore, the use of Riemannian geometry principles as in [42] cannot be envisioned.

In the following, we propose, study and compare three different strategies to get a relevant depiction of our deformations in a linear space. The first two ones are physically/mechanically-oriented and can be viewed as non-straightforward adaptations of [63]. While in [63] the shapes are modelled through their boundaries and subject to boundary stresses, our framework involves the whole image as the object to be deformed. Note that with the prescribed boundary conditions $\forall i \in \{1, \cdots, M\}, \varphi_i = \text{Id}$ on $\partial\Omega$, no boundary stress is applied, and subsequently, inner volumetric stresses are considered. This constitutes a major difference with the work of [63] and as demonstrated in section SM4, it entails substantial adaptations in the mathematical developments. The first two methods rely on fundamental notions of elastic behaviour and the following observation made in [63] : "*the classical covariance tensor can be identified with the covariance tensor of the displacements obtained by adding a small fraction of the i-th spring force under the Hooke's law*". Whilst the first method is based on the linearisation of the stored energy function around the identity, which might result in the loss of the initial nonlinear nature of the deformations but has the advantage of being fast, the second approach is more intricate. It retrieves the whole nature of the deformations by performing the PCA on the Cauchy stress tensors, relying on the locally underlying one-to-one relation between this tensor and the deformation, but requires the resolution of a highly nonlinear and non convex problem similar to the one studied previously to come back to the deformation space.

Our goal was to design an alternative method that would be a good compromise between rendering the nonlinear nature of the deformation and in terms of numerical complexity. The first objective is achieved by handling both the deformation field and the deformation tensors that encode the local deformation state resulting from stresses. This constitutes another novelty of the proposed work. We have moved toward a completely different point of view since the problem is no longer explored as a physical one but is now identified as an approximation one in the $D^m$-spline setting ([3]). The first two methods thus serve as benchmark to assess the interest of this new vision.
Due to page number limitation and as the third method proves to be a proper trade-off between ability to reproduce the nonlinear nature of the deformations and intermediate computation time, the mathematical details of the first two methods are postponed in the document *Supplementary Materials* section SM4 and we only focus on the third one.

**4.1. Third Approach: Approximation Modelling.** This section is devoted to the analysis of a novel method in which the linear representation problem is seen as an approximation one in the $D^m$-spline setting. Since the deformation tensor suitably characterises the local deformation (amplitude, direction, etc.), we aim at finding an appropriate approximation of our deformations in a linear space $H^3(\Omega, \mathbb{R}^2)$ that also approximates well the deformation tensors. For the sake of clarity, we omit the indices $i$ in the following. Henceforth, $u$ denotes the displacement field related to $\varphi^{-1}$, inverse deformation field obtained at the outcome of the first algorithmic stage. We consider the following problem in the two-dimensional case —case of interest in the numerical part, but straightforwardly extendable to 3D —:

$$ (4.1) \qquad \min_{v \in H^3(\Omega, \mathbb{R}^2)} \epsilon|v|_{3,\Omega,\mathbb{R}^2}^2 + \frac{\gamma}{2} \langle \rho(\nabla v + \nabla v^T - \nabla u - \nabla u^T - \nabla u^T \nabla u) \rangle_{M_2(\mathbb{R}),N}^2 \\ + \langle \zeta(v-u) \rangle_{\mathbb{R}^2,N}^2, $$

where $|\cdot|_{3,\Omega,\mathbb{R}^2}$ is the semi-norm on $H^3(\Omega, \mathbb{R}^2)$, $\zeta: \begin{array}{ccc} H^3(\Omega, \mathbb{R}^2) & \to & (\mathbb{R}^2)^N \\ v & \mapsto & \zeta(v) = (v(a_1), \dots, v(a_N))^T \end{array}$, $\rho: \begin{array}{ccc} H^2(\Omega, M_2(\mathbb{R})) & \to & (M_2(\mathbb{R}))^N \\ v & \mapsto & \rho(v) = (v(a_1), \dots, v(a_N))^T \end{array}$. Also, $a_1, \dots, a_N$ denote the image pixel coordinates with $N$ the total number of pixels, and $\forall \xi \in (M_2(\mathbb{R}))^N, \forall \eta \in (M_2(\mathbb{R}))^N$, $\langle \xi, \eta \rangle_{M_2(\mathbb{R}),N} = \sum_{i=1}^N \xi_i : \eta_i$, while $\forall \xi \in (\mathbb{R}^2)^N, \forall \eta \in (\mathbb{R}^2)^N$, $\langle \xi, \eta \rangle_{\mathbb{R}^2,N} = \sum_{i=1}^N \xi_i^T \eta_i$. However, stated as it is, the problem is not well defined since $u \in W_0^{1,\infty}(\Omega, \mathbb{R}^2)$ and does not belong to $C^1(\Omega, \mathbb{R}^2)$ preventing us from extracting isolated values of $\nabla u$. Therefore, for the theoretical analysis of the model, we introduce $(f_k) \in C_0^\infty(\Omega, \mathbb{R}^2) \cap W^{1,\infty}(\Omega, \mathbb{R}^2)$, the sequence from the density result such that,

$$ f_k \underset{k \to +\infty}{\longrightarrow} u \text{ in } W^{1,\infty}(\Omega, \mathbb{R}^2). $$

In practice however, we solve (4.1) and we give details on the implementation in section SM12.

**Remark 4.1.** An alternative approach would consist in using Lebesgue-Besicovitch differentiation theorem that states that for almost every point, the value of an integrable function is the limit of infinitesimal averages taken about the point.

Let $A_0 = \{a_i\}_{i=1,\dots,N_0}$ and $A_1 = \{b_i\}_{i=1,\dots,N_1}$ be two sets of $N_0$ and $N_1$ points of $\bar{\Omega}$ respectively, containing both a $P^1$-unisolvent subset. Let us denote by $\rho_0$ the operator defined by

$$ \rho_0: \begin{array}{ccc} H^3(\Omega, \mathbb{R}^2) & \to & (\mathbb{R}^2)^{N_0} \\ v & \mapsto & \rho_0(v) = (v(a_i))_{i=1,\dots,N_0}^T, \end{array} $$
and by $\rho_1$ the operator defined by

$$ \rho_1 : \left| \begin{array}{l} H^2(\Omega, M_2(\mathbb{R})) \to (M_2(\mathbb{R}))^{N_1} \\ v \mapsto \rho_1(v) = (v(b_i))_{i=1, \dots, N_1}^T \end{array} \right. . $$

We introduce the functionals

$$ \mathcal{F}_{\epsilon,k} : \left\{ \begin{array}{l} H^3(\Omega, \mathbb{R}^2) \to \mathbb{R} \\ v \mapsto \langle \rho_0(v - f_k) \rangle_{\mathbb{R}^2, N_0}^2 + \frac{\gamma}{2} \langle \rho_1(\nabla v + \nabla v^T - \nabla f_k - \nabla f_k^T - \nabla f_k^T \nabla f_k) \rangle_{M_2(\mathbb{R}), N_1}^2 , \\ + \epsilon |v|_{3, \Omega, \mathbb{R}^2}^2 \end{array} \right. $$

and consider the problem

$$ (4.2) \qquad \left\{ \begin{array}{l} \text{Search for } u_{\epsilon} \in H^3(\Omega, \mathbb{R}^2) \text{ such that:} \\ \forall v \in H^3(\Omega, \mathbb{R}^2), \mathcal{F}_{\epsilon,k}(u_{\epsilon}) \le \mathcal{F}_{\epsilon,k}(v) \end{array} \right. . $$

We omit the explicit dependency of $u_\epsilon$ on $k$. In the sequel, we theoretically study the model (4.2), and start by proving its equivalence with a variational formulation.

**Theorem 4.2 (Equivalence of problems).** The problem (4.2) is equivalent to the following variational problem:

$$ (4.3) \quad \left\{ \begin{array}{l} \text{Search for } u_{\epsilon} \in H^3(\Omega, \mathbb{R}^2) \text{ such that } \forall v \in H^3(\Omega, \mathbb{R}^2), \\ \langle \rho_0(u_{\epsilon}), \rho_0(v) \rangle_{\mathbb{R}^2, N_0} + \frac{\gamma}{2} \langle \rho_1(\nabla u_{\epsilon} + \nabla u_{\epsilon}^T), \rho_1(\nabla v + \nabla v^T) \rangle_{M_2(\mathbb{R}), N_1} + \epsilon(u_{\epsilon}, v)_{3, \Omega, \mathbb{R}^2} \\ = \langle \rho_0(v), \rho_0(f_k) \rangle_{\mathbb{R}^2, N_0} + \frac{\gamma}{2} \langle \rho_1(\nabla v + \nabla v^T), \rho_1(\nabla f_k + \nabla f_k^T + \nabla f_k^T \nabla f_k) \rangle_{M_2(\mathbb{R}), N_1}. \end{array} \right. $$

*Proof.* The detailed proof is available in section SM6.

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

We now define a new norm equivalent to the classical norm on $H^3(\Omega, \mathbb{R}^2)$, which will be useful in the following. We make the dependency on the set $A_0$ explicit, while the set $A_1$ is supposed to be fixed once and for all.

**Lemma 4.3 (Equivalence of norms).** The mapping defined by

$$ \| \cdot \|_{A_0, 3, \Omega, \mathbb{R}^2} : \left\{ \begin{array}{l} H^3(\Omega, \mathbb{R}^2) \to \mathbb{R} \\ f \mapsto \|f\|_{A_0, 3, \Omega, \mathbb{R}^2} = (\langle \rho_0(f) \rangle_{\mathbb{R}^2, N_0}^2 + \langle \rho_1(\nabla f + \nabla f^T) \rangle_{M_2(\mathbb{R}), N_1}^2 + |f|_{3, \Omega, \mathbb{R}^2}^2)^{\frac{1}{2}} , \end{array} \right. $$

is a Hilbert norm equivalent to the norm $\| \cdot \|_{3, \Omega, \mathbb{R}^2}$ in $H^3(\Omega, \mathbb{R}^2)$.

*Proof.* The detailed proof is given in section SM7.

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

We are now able to prove the existence and uniqueness of the minimizer.

**Theorem 4.4 (Existence and uniqueness of a minimizer).** The variational problem (4.3) admits a unique solution.

*Proof.* The detailed proof is given in section SM8.

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

We now focus on a convergence result. Let $D$ be a subset of $]0, +\infty[$ for which $0$ is an accumulation point. For any $d \in D$, let $A^d$ be a set of $N = N(d)$ distinct points from $\bar{\Omega}$ that
contains a $P^1$-unisolvent subset. We assume that $\sup_{x \in \Omega} \delta(x, A^d) = d$, where $\delta$ is the Euclidean distance in $\mathbb{R}^2$. Thus $d$ is the radius of the biggest sphere included in $\Omega$ that contains no point from $A^d$. Also $d$ is bounded and $\limsup_{d \to 0} \delta(x, A^d) = 0$. For any $d \in D$, let us denote by $\rho^d$ the mapping defined by

$$ \rho^d : \begin{cases} H^3(\Omega, \mathbb{R}^2) \to (\mathbb{R}^2)^N \\ v \mapsto \rho^d(v) = ((v(a))_{a \in A^d})^T \end{cases}, $$

and by $\|\cdot\|_{A^d, 3, \Omega, \mathbb{R}^2}$, the norm defined by

$$ \|f\|_{A^d, 3, \Omega, \mathbb{R}^2} = [\langle \rho^d(f) \rangle_N^2 + \langle \rho_1(\nabla f + \nabla f^T) \rangle_{M_2(\mathbb{R}), N_1}^2 + |f|_{3, \Omega, \mathbb{R}^2}]^{\frac{1}{2}}. $$

As shown in the previous lemma, $\|\cdot\|_{A^d, 3, \Omega, \mathbb{R}^2}$ is equivalent to the norm $\|\cdot\|_{3, \Omega, \mathbb{R}^2}$ in $H^3(\Omega, \mathbb{R}^2)$ —but not uniformly in $d$—.

**Lemma 4.5.** Let $B_1 = \{b_{01}, \dots, b_{0,N}\}$ be a fixed $P^1$-unisolvent subset of $\bar{\Omega}$. By hypothesis, $0 \in \bar{D}$, and $\limsup_{d \to 0} \delta(x, A^d) = 0$ holds, so

$$ \forall j = 1, \dots, N, \exists (a_{0,j}^d)_{d \in D}, (\forall d \in D, a_{0,j}^d \in A^d) \text{ and } b_{0j} = \lim_{d \to 0} a_{0j}^d. $$

For any $d \in D$, let $A_0^d$ be the set $\{a_{01}^d, \dots, a_{0,N}^d\}$ and let $\|\cdot\|_{A_0^d, 3, \Omega, \mathbb{R}^2}$ be the norm defined by $\|f\|_{H^3(\Omega, \mathbb{R}^2)}$,

$$ \|f\|_{A_0^d, 3, \Omega, \mathbb{R}^2} = \left[ \sum_{j=1}^{N} \langle f(a_{0j}^d) \rangle_{\mathbb{R}^2}^2 + \sum_{i=1}^{N_1} \| \nabla f(b_i) + \nabla f(b_i)^T \|^2 + |f|_{3, \Omega, \mathbb{R}^2}^2 \right]^{\frac{1}{2}}. $$

Then there exists $\eta > 0$ such that for any $d \le \eta$, $\|\cdot\|_{A_0^d, 3, \Omega, \mathbb{R}^2}$ is a norm on $H^3(\Omega, \mathbb{R}^2)$ uniformly equivalent on $D \cap [0, \eta]$ to the norm $\|\cdot\|_{3, \Omega, \mathbb{R}^2}$.

*Proof.* The detailed proof is available section SM9. ■

Equipped with this result, we are able to prove a convergence result on the following problem:

$$ (4.4) \quad \left\{ \begin{array}{ll} \text{Search for } u_{\epsilon}^d \in H^3(\Omega, \mathbb{R}^2) \text{ such that } \forall v \in H^3(\Omega, \mathbb{R}^2), & \\ \langle \rho^d(u_{\epsilon}^d - f_k) \rangle_{\mathbb{R}^2, N}^2 + \frac{\gamma}{2} \langle \rho_1(\nabla u_{\epsilon}^d + (\nabla u_{\epsilon}^d)^T - \nabla f_k - \nabla f_k^T - \nabla f_k^T \nabla f_k) \rangle_{M_2(\mathbb{R}), N_1}^2 & \\ + \epsilon |u_{\epsilon}^d|_{3, \Omega, \mathbb{R}^2}^2 \le \langle \rho^d(v - f_k) \rangle_{\mathbb{R}^2, N}^2 + \frac{\gamma}{2} \langle \rho_1(\nabla v + (\nabla v)^T - \nabla f_k - \nabla f_k^T) \rangle_{M_2(\mathbb{R}), N_1}^2 & \\ - \nabla f_k^T \nabla f_k \rangle_{M_2(\mathbb{R}), N_1}^2 + \epsilon |v|_{3, \Omega, \mathbb{R}^2}^2, & \end{array} \right. $$

**Theorem 4.6 (Convergence).** For any $d \in D$, we denote by $u_{\epsilon}^d$ the unique solution to problem (4.4) for $\epsilon$ fixed. Then under the above assumptions, there exists a subsequence $(u_{\epsilon}^{d_l})$ with $\lim_{l \to +\infty} d_l = 0$ such that

$$ u_{\epsilon}^{d_l} \underset{l \to +\infty}{\to} f_k $$

in $H^3(\Omega, \mathbb{R}^2)$, and $\lim_{k \to +\infty} \lim_{l \to +\infty} \|u_{\epsilon}^{d_k} - u\|_{1, \infty} = 0$.
*Proof.* The detailed proof is available in section SM10.

An alternative convergence study is given in the *Supplementary material* section SM11.
A description of the numerical resolution is provided in section SM12.

**Remark 4.7.** A classical PCA is then performed on the obtained displacement fields using the $L^2$ scalar product for the covariance operator, i.e. $(C_{i,j})_{\substack{i=1,\dots,M \\ j=1,\dots,M}} = \int_{\Omega} v_{1,j}v_{1,i} + v_{2,j}v_{2,i} dx$, where $(v_{1,i} v_{2,i})^T$ for $i = 1, \dots, M$ are the displacement fields obtained in a linear space for each image, and the resulting displacement fields are denoted by $(v_{1,\mathrm{pca},j}, v_{2,\mathrm{pca},j})^T$, where $j$ stands for the mode indexation.

# 5. Numerical Simulations.

**5.1. General framework.** This section is devoted to the analysis of numerical experiments. First, on a set of 19 binary images 'device8-1' from the MPEG7 shape database (http://www.dabi.temple.edu/~shape/MPEG7/dataset.html) and then on medical images: cardiac MRI made of 8 frames per slice of size 150x150, the first image of the sequence reflecting the case where the heart is most dilated (end diastole - ED), while the last one illustrating the case where the heart is most contracted (end systole - ES), and liver dynamic MRI made of 14 frames per slice of size 195x166.

The computations have been made on an Intel Core i7 computer with 2.60 GHz and 8 GB memory, using MUMPS packages, in a C implementation.

The question of assessing the proposed model encompasses several angles of inquiry:

(i) the qualitative evaluation of the obtained atlas in comparison to a sequential treatment of the segmentation and registration tasks (note that the model requires the segmentation step be processed first since it involves the penalisation $\|\theta_{\hat{T}_i} - \theta_{T_i} \circ \varphi_i + \theta_R\|_{L^2(\Omega)}^2$). We do not question here the relevancy of the regularisation on the $\varphi_i$'s since in practice, the $\varphi_i$'s need to be invertible (as in the expression $\theta_{T_i} \circ \varphi_i$, $\theta_{T_i}$ is also an unknown), which is guaranteed with the proposed regulariser. As the average shape is an unknown of the problem and ground truth is not provided, the evaluation of the method itself primarily relies on visual inspection and empirical arguments consistent with the biological phenomena involved;

(ii) the evaluation of the PCA in capturing strongly nonlinear geometric variations.

These two main levels of discussion dictate the structure of the section. Each subsection focuses on a specific dataset and provides both the atlas generated by our joint model, together with the first principal modes of variation. The gain of the combined approach in terms of sharp edges and in removing ghosting artefacts (blurring / splitting into two effects) is emphasised compared to a sequential treatment, as well as the accuracy of our proposed model in reflecting the high nonlinear geometric variations in comparison to the two more physically/mechanically oriented methods. This stage is delicate since as stated below, it requires comparing the three methodologies in the most efficient and impartial manner, and subsequently, setting the involved tuning parameters (e.g., the weight balancing the loading forces for the Cauchy stress tensor based PCA) adequately. That said, results demonstrate nevertheless and as expected, that linearised elasticity produces very small displacements, which tends to favor comparisons with the Cauchy stress tensor based PCA.
For the sake of reproducibility, we provide in Table 1 the values of the tuning parameters. The coefficients $a_1$, $a_2$ and $a_3$ involved in the Ogden stored energy function affect respectively the averaged local change of length, and the averaged local change of area, impacting subsequently on the rigidity of the deformation. The higher the $a_i$'s are, the more rigid the deformation is. The ranges of these parameters are rather stable for the medical experiments. Parameters $\gamma_T$ and $\gamma_R$ weighting the $L^0$ norm control thus the balance between the $L^0$-component and the $L^2$-penalisation together with $\lambda_R$ and $\lambda_T$. High values of $\gamma_T/\gamma_R$ favor few large partitions, while small values yield an approximation exhibiting more jumps. $\lambda_R$, $\gamma_T$, and $\gamma_1$ weight the fidelity term in the registration task, and thus the higher they are, the closer the deformed templates and deformed segmentations are to the mean segmentation. The visualisation of the main modes of variation is done as follows. We denote by $(v_{1,pca_{1,i}}, v_{2,pca_{1,i}})^T$ the $i$th resulting displacement field (related to the $i$th mode of variation) from the first method based on the linearisation of the stored energy function around the identity, providing in practice and as expected extremely small displacements. $\sigma_{pca_{2,i}}$ represents the $i$th resulting Cauchy tensor from the second method and $(v_{1,pca_{2,i,\delta}}, v_{2,pca_{2,i,\delta}})^T$ the associated displacement field obtained with weighting parameter $\delta$ balancing the inner forces chosen equal to 1, 2, 0.3 for the T-shape, the liver and the heart respectively, that prove to be suitable parameters for an unbiased analysis. At last, $(v_{1,pca_{3,i}}, v_{2,pca_{3,i}})^T$ stands for the $i$th resulting displacement field from the third method based on approximation theory. We propose visualising the $i$th mode of variation by showing $\theta_R \circ (\text{Id} + 5c.10^6 \begin{pmatrix} v_{1,pca_{1,i}} \\ v_{2,pca_{1,i}} \end{pmatrix})$, $\theta_R \circ (\text{Id} + 5c \begin{pmatrix} v_{1,pca_{2,i,\delta}} \\ v_{2,pca_{2,i,\delta}} \end{pmatrix})$, and $\theta_R \circ (\text{Id} + 50c \begin{pmatrix} v_{1,pca_{3,i}} \\ v_{2,pca_{3,i}} \end{pmatrix})$, for each method respectively, with $c$ varying from -5 to 5. The parameters are chosen in order to make the comparison as fair as possible. The computation times for each method and each example are provided in Table 2.

<table><caption>**Table 1**

Parameters.</caption><thead><tr><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>α</th><th>β</th><th>γ<sub>R</sub></th><th>γ<sub>T</sub></th><th>λ<sub>T</sub></th><th>λ<sub>R</sub></th><th>dt</th><th>nbIter</th></tr></thead><tbody><tr><td>T-shape</td><td>1</td><td>5.10<sup>3</sup></td><td>0.01</td><td>1</td><td>8.10<sup>4</sup></td><td>1</td><td>10</td><td>100</td><td>3</td><td>0.5</td><td>1</td><td>1</td><td>0.001</td><td>100</td></tr><tr><td>Heart ED(108)-ES(101)</td><td>5</td><td>1.10<sup>3</sup></td><td>4</td><td>1</td><td>8.10<sup>4</sup></td><td>1</td><td>100</td><td>100</td><td>0.02</td><td>0.03</td><td>1.5</td><td>1.5</td><td>0.01</td><td>500</td></tr><tr><td>Liver - slice 12</td><td>5</td><td>1.10<sup>3</sup></td><td>4</td><td>1.5</td><td>8.10<sup>4</sup></td><td>1</td><td>10</td><td>100</td><td>0.05</td><td>0.05</td><td>1.5</td><td>1.5</td><td>0.01</td><td>100</td></tr></tbody></table>

<table><caption>**Table 2**

Execution times</caption><thead><tr><th>Execution time</th><th>Atlas generation</th><th>PCA 1<sup>st</sup> method</th><th>PCA 2<sup>nd</sup> method</th><th>PCA 3<sup>rd</sup> method</th></tr></thead><tbody><tr><td>T-shape</td><td>7 min</td><td>2 sec</td><td>4 min</td><td>44 sec</td></tr><tr><td>Heart ED(108)-ES(101)</td><td>49 min</td><td>3 sec</td><td>10 min</td><td>3 min 27 sec</td></tr><tr><td>Liver slice 12</td><td>26 min</td><td>6 sec</td><td>7 min 28 sec</td><td>13 min</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Figure 3.** *Input images.*</figcaption></figure>

**5.2. T-shape example: 19 images.** The proposed method is first evaluated on a synthetic example (Figure 3) to emphasise the ability of the model to generate large deformations and to produce a physically sound average shape. As depicted in Figure 4, contrary to a sequential treatment of the tasks, the joint model creates a mean object with sharp edges perfectly matched by the deformed templates, see particularly the bottom of the T shape. The joint approach also tends to better preserve the original contrast of the images than the sequential approach. The principal component analysis (Figure 5) shows that the first two modes of variation have an effect on the undulation of both the vertical and horizontal bars (the first mode tends to represent the undulation in the direction *bottom-left corner-top-right corner* whereas the second one tends to capture the undulation in the direction *top-left corner-bottom-right corner*), while the third mode affects the thickness of the vertical bar particularly at the junction. The fourth mode, for its part, acts more locally on the curvature of the envelope of the shape, especially on the lower part of the horizontal bar that exhibits cavities. A first observation is that the proposed model allows uncorrelating the main tendencies, which is what is expected from such an analysis. It seems (Figure 6) that this decoupling property is not as well exemplified when applying the linearisation around identity which is not able to recover the nonlinear variations coming from the undulation or the Cauchy stress tensor based PCA mixing the first, the third and the fourth modes of variation from our third approach based on approximation modelling.

**5.3. Liver example: Slice 12 [74].** The second example is dedicated to right lobe liver dynamic MRI http://www.vision.ee.ethz.ch/~organmot/chapter_download.shtml, the shape of this organ being influenced by the surrounding structures such as the diaphragm ([73]). It exemplifies the ability of the method to alleviate ghosting artefacts, in particular splitting into two effects, through parameter $\gamma_R$ that influences the number of phases (Figure 7). The obtained atlas exhibits sharp edges with fewer artefacts (Figure 8) than with a sequential approach. The statistical analysis is then performed (Figure 9). The first mode of variation encodes the motion in superior/inferior direction of the liver, which is consistent with
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Comparison of the atlases generated by our joint model, and by a sequential approach with our discrepancy measure. Visual assessment in terms of blurring artefacts, and contrast are pointed out with yellow arrows.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 5.** First four modes of variation obtained with our method based on approximation modelling via contour representations.</figcaption></figure>

the physics ([73]): respiration is largely governed by the diaphragm, and the liver, located beneath the diaphragm, is thus strongly influenced by breathing (pushed downwards when the diaphragm is contracted and upwards when expanded). The same applies to the kidney located in the right lower part. As for the second mode of variation, it reflects the transversal motion, particularly visible on the right excrescences exhibited by the liver and the left wall. A comparison of the first mode of variation obtained with both the linearisation around the identity and the Cauchy stress tensor based PCA shows that our method reflects the geometric variability slightly better in the superior/inferior direction (Figure 10) particularly in the right lower part of the kidney.

**5.4. Heart example: frames ED(108)-ES(101).** Lastly, the joint algorithm is applied to the set of MRI representing a cardiac cycle. The obtained segmentations (Figure 11) are piecewise constant approximations of the initial images and reflect well the geometrical shapes of the template images. The deformed segmentations (involving large deformations) are close to each other and well aligned to the mean reference. The obtained atlas exhibits sharp edges, contrary to the result produced by a sequential treatment (Figure 12) that shows blurry artefacts. It is also anatomically consistent: the mean reference corresponds to a compromise between the full expansion and the full contraction, which is reasonable from a biological
<figure><img src="image_15.png" /></figure>

First method : Linearisation around the identity

Second method : PCA on the Cauchy tensor

Third method : Approximation modelling

**Figure 6.** *Comparison of the first mode of variation from three different methods using contour representations.*

standpoint. The second step of the algorithm is then applied. The obtained first modes of variation are again consistent with the anatomical dynamic of the heart (Figure 13): while the first mode encodes the dilation/contraction of the right ventricular chamber in the transverse direction, the second mode of variation conveys the vertical stretching of the left ventricular chamber. Our approximation-based PCA thus allows to disconnect these two movements, this property being less visible in the case of the Cauchy stress tensor based PCA (Figure 14) for which the first mode of variation encodes both tendencies.

**6. Conclusion.** This paper addressed the twofold question of finding an average representative of a dataset of different subjects and deriving then some statistics by identifying the main modes of variation. To achieve this goal, the problem is envisioned as a joint registration/segmentation one, based on nonlinear elasticity concepts. Once a linearisation is performed in order for the displacement fields to live in a vector space, a PCA is investigated to capture the meaningful geometric variations. The computational feasibility is exemplified through several applications, demonstrating that the generated atlas encodes the fine geometrical structures and exhibits sharp edges with fewer ghosting artefacts, while the proposed approximation based PCA uncorrelates properly the more significant tendencies. This work paves the way to several applications among which: (i) motion-correction problem from a set of multiple acquisitions corrupted by motion in this joint reconstruction and registration framework or (ii) multiscale segmentation/registration problem to extract the deformation pairing the structures (i.e. viewed as global deformations) and an additional deformation reflecting the more local variability.

REFERENCES

[1] L. AMBROSIO AND J. BERTRAND, *DC calculus*, Mathematische Zeitschrift, 288 (2018), pp. 1037–1080.

[2] J.-H. AN, Y. CHEN, F. HUANG, D. WILSON, AND E. GEISER, *MICCAI*, 2005, ch. A Variational PDE Based Level Set Method for a Simultaneous Segmentation and Non-Rigid Registration, pp. 286–293.

[3] R. ARCANGÉLI, M. LÓPEZ DE SILANES, AND J. TORRENS, *Multidimensional minimizing splines. Theory and applications*, Grenoble Sciences, Kluwer Academic Publishers, Boston, MA, 2004.

[4] V. ARSIGNY, O. COMMOWICK, X. PENNEC, AND N. AYACHE, *A log-Euclidean framework for statistics on diffeomorphisms*, in MICCAI, Springer, 2006, pp. 924–931.

[5] J. ASHBURNER AND K. J. FRISTON, *Nonlinear spatial normalization using basis functions*, Hum. Brain Mapp., 7 (1999), pp. 254–266.
<figure><img src="image_16.png" /><figcaption>Sum of all the difference maps between : (a) the initial segmentations and the segmented atlas;
(b) the deformed segmentations and the segmented atlas.</figcaption></figure>

<caption>**Figure 7.** Some of the initial slices (1,6,11,14) with the associated segmentations given by our joint algorithm underneath. The last row displays the sum of difference maps between the initial segmentation and the segmented atlas, the deformed/registered segmentations and the segmented atlas, and the segmented atlas generated by our joint model.</caption>

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

[7] G. AUBERT AND P. KORNPROBST, *Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations*, Springer-Verlag, 2001.

[8] A. I. AVILES-RIVERO, G. WILLIAMS, M. J. GRAVES, AND C.-B. SCHÖNLIEB, *Compressed sensing plus motion (CS+M): A new perspective for improving undersampled MR image reconstruction*, 2018, https://arxiv.org/abs/1810.10828.

[9] J. M. BALL, *Global invertibility of Sobolev functions and the interpenetration of matter*, P. Roy. Soc. Edin. A, 88 (1981), p. 315328.

[10] M. BEG, M. MILLER, A. TROUVÉ, AND L. YOUNES, *Computing large deformation metric mappings via*
<figure><img src="image_17.png" /><figcaption>**Figure 8.** Comparison of the atlases generated by our joint model, and by a sequential approach with our discrepancy measure. Visual assessment in terms of blurring and ghosting artefacts are indicated by yellow arrows.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 9.** First two modes of variation obtained with our method based on approximation modelling. The horizontal lines help visualise the vertical movement of the liver and the kidney.</figcaption></figure>

geodesic flows of diffeomorphisms, Int. J. Comput. Vis., 61 (2005), pp. 139–157.

[11] M. BLUME, A. MARTINEZ-MOLLER, A. KEIL, N. NAVAB, AND M. RAFECAS, *Joint Reconstruction of Image and Motion in Gated Positron Emission Tomography*, IEEE Trans. Med. Imag., 29 (2010), pp. 1892–1906.

[12] Y. BOINK, *Combined modelling of optimal transport and segmentation revealing vascular properties*, 2016.

[13] C. BROIT, *Optimal Registration of Deformed Images*, PhD thesis, Computer and Information Science, University of Pennsylvania, 1981.

[14] C. BRUNE, H. MAURER, AND M. WAGNER, *Detection of intensity and motion edges within optical flow via multidimensional control*, SIAM J. Imaging Sci., 2 (2009), pp. 1190–1210.

[15] M. BURGER, H. DIRKS, L. FRERKING, A. HAUPTMANN, T. HELIN, AND S. SILTANEN, *A variational reconstruction method for undersampled dynamic X-ray tomography based on physical motion models*, Inverse Problems, 33 (2017), p. 124008.
<figure><img src="image_19.png" /><figcaption>**Figure 10.** *Comparison of the first mode of variation for 3 different methods.*</figcaption></figure>

[16] M. BURGER, H. DIRKS, AND C. SCHÖNLIEB, *A variational model for joint motion estimation and image reconstruction*, SIAM J. Imaging Sci., 11 (2018), pp. 94–128.

[17] M. BURGER, J. MODERSITZKI, AND L. RUTHOTTO, *A hyperelastic regularization energy for image registration*, SIAM J. Sci. Comput., 35 (2013), pp. B132–B148.

[18] E. CAZELLES, V. SEGUY, J. BIGOT, M. CUTURI, AND N. PAPADAKIS, *Geodesic PCA versus log-PCA of histograms in the Wasserstein space*, SIAM J. Sci. Comput., 40 (2018), pp. B429–B456.

[19] A. CHAMBOLLE, D. CREMERS, AND T. POCK, *A convex approach to minimal partitions*, SIAM J. Imaging Sci., 5 (2012), pp. 1113–1158.

[20] H.-L. CHAN, S. YAN, L.-M. LUI, AND X.-C. TAI, *Topology-Preserving Image Segmentation by Beltrami Representation of Shapes*, J. Math. Imaging Vis., 60 (2018), pp. 401–421.

[21] G. E. CHRISTENSEN, *Deformable Shape Models for Anatomy*, PhD thesis, Washington University, Sever Institute of technology, USA, 1994.

[22] S. Y. CHUN AND J. A FESSLER, *Joint Image Reconstruction and Nonrigid Motion Estimation with a Simple Penalty that Encourages Local Invertibility*, vol. 7258, 02 2009.

[23] P. CIARLET, *Elasticité Tridimensionnelle*, Masson, 1985.

[24] O. CLATZ, M. SERMESANT, P.-Y. BONDIAU, H. DELINGETTE, S. K. WARFIELD, G. MALANDAIN, AND N. AYACHE, *Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion
<figure><img src="image_20.png" /><figcaption>Sum of all the difference maps between : (a) the initial segmentations and the segmented atlas; (b) the deformed segmentations and the segmented atlas.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 12.** Comparison of the atlases generated by our joint model, and by a sequential approach with our discrepancy measure. Visual assessment in terms of blurring artefacts are pointed out with yellow arrows.</figcaption></figure>
First Mode of Variation, eigenvalue : 85 422

<figure><img src="image_22.png" /><figcaption>**Figure 13.** First three modes of variation obtained with our method based on approximation modelling.
The red circles shows the moving part: dilation/contraction of the right ventricular chamber for the first mode,
and vertical stretching of the left ventricular chamber for the second mode.</figcaption></figure>

First Mode of Variation

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

Third method : Approximation modelling

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

Second method : PCA on the Cauchy tensor

<caption>**Figure 14.** Comparison of the first mode of variation for 2 different methods.</caption>

with *biomechanical deformation*, IEEE Trans. Med. Imaging, 24 (2005), pp. 1334–1346.

[25] V. CORONA, A. AVILES-RIVERO, N. DEBROUX, M. GRAVES, C. LE GUYADER, C.-B. SCHÖNLIEB, AND G. WILLIAMS, *Multi-tasking to Correct: Motion-Compensated MRI via Joint Reconstruction and Registration*, in Scale Space and Variational Methods in Computer Vision, 2019.

[26] C. DAVATZIKOS, *Spatial transformation and registration of brain images using elastically deformable mod-
els, Comput. Vis. Image Underst., 66 (1997), pp. 207-222.

[27] M. H. DAVIS, A. KHOTANZAD, D. P. FLAMIG, AND S. E. HARMS, *A physics-based coordinate transformation for 3-D image matching*, IEEE Trans. Med. Imaging, 16 (1997), pp. 317-328.

[28] N. DEBROUX AND C. LE GUYADER, *A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors*, SIAM J. Imaging Sci., 11 (2018), pp. 957-990.

[29] N. DEBROUX, S. OZERÉ, AND C. LE GUYADER, *A non-local topology-preserving segmentation guided registration model*, J. Math. Imaging Vision, (2017), pp. 1-24.

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

[31] R. DERFOUL AND C. LE GUYADER, *A relaxed problem of registration based on the Saint Venant-Kirchhoff material stored energy for the mapping of mouse brain gene expression data to a neuroanatomical mouse atlas*, SIAM J. Imaging Sci., 7 (2014), pp. 2175-2195.

[32] M. DROSKE AND M. RUMPF, *A variational approach to non-rigid morphological registration*, SIAM J. Appl. Math., 64 (2004), pp. 668-687.

[33] M. DROSKE AND M. RUMPF, *Multiscale joint segmentation and registration of image morphology*, IEEE Trans. Pattern Anal. Mach. Intell., 29 (2007), pp. 2181-2194.

[34] L. EVANS AND R. GARIEPY, *Measure Theory and Fine Properties of Functions*, CRC Press, 1992.

[35] B. FISCHER AND J. MODERSITZKI, *Curvature based image registration*, J. Math. Imaging Vis., 18 (2003), pp. 81-85.

[36] A. GOOYA, K. POHL, M. BILELLO, L. CIRILLO, G. BIROS, E. MELHEM, AND C. DAVATZIKOS, *GLISTR: Glioma Image Segmentation and Registration*, IEEE Trans. Med. Imaging, 31 (2012), pp. 1941-1954.

[37] J. HAN, B. BERKELS, M. DROSKE, J. HORNEGGER, M. RUMPF, C. SCHALLER, J. SCORZIN, AND H. URBACH, *Mumford-Shah Model for One-to-One Edge Matching*, IEEE Trans. Image Process., 16 (2007), pp. 2720-2732.

[38] B.-W. HONG, S. SOATTO, AND L. A. VESE, *Enforcing local context into shape statistics*, Adv. Comput. Math., 31 (2009), pp. 185-213.

[39] S. JOSHI, B. DAVIS, M. JOMIER, AND G. GERIG, *Unbiased diffeomorphic atlas construction for computational anatomy*, NeuroImage, 23 (2004), pp. S151-S160.

[40] M. KASS, A. WITKIN, AND D. TERZOPOULOS, *Snakes: Active contour models*, Int. J. COMPUT. VISION, 1 (1988), pp. 321-331.

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

[42] H. LAGA, I. H. JERMYN, S. KURTEK, AND A. SRIVASTAVA, *Elastic 3D shape analysis using square-root normal field representation*, in CDC, 2017 IEEE 56th Annual Conference, 2017, pp. 2711-2717.

[43] K. LAM AND L. LUI, *Landmark- and Intensity-Based Registration with Large Deformations via Quasi-conformal Maps*, SIAM J. Imaging Sci., 7 (2014), pp. 2364-2392.

[44] K. C. LAM AND L. M. LUI, *Quasi-Conformal Hybrid Multi-modality Image Registration and its Application to Medical Image Fusion*, in Advances in Visual Computing, 2015, pp. 809-818.

[45] H. LE DRET, *Notes de Cours de DEA. Méthodes mathématiques en élasticité*, 2003-2004.

[46] C. LE GUYADER, D. APPRATO, AND C. GOUT, *Spline approximation of gradient fields: Applications to wind velocity fields*, Math. Comput. Simul., 97 (2014), pp. 260-279.

[47] C. LE GUYADER AND L. VESE, *Self-Repelling Snakes for Topology-Preserving Segmentation Models*, IEEE Trans. Image Process., 17 (2008), pp. 767-779.

[48] C. LE GUYADER AND L. VESE, *A combined segmentation and registration framework with a nonlinear elasticity smoother*, Comput. Vis. Image Underst., 115 (2011), pp. 1689-1709.

[49] Y. T. LEE, K. C. LAM, AND L. M. LUI, *Landmark-Matching Transformation with Large Deformation Via N-dimensional Quasi-conformal Maps*, J. Sci. Comput., 67 (2016), pp. 926-954.

[50] O. LEHTO AND K. VIRTANEN, *Quasiconformal Mappings in the Plane*, Springer-Verlag, 1973.

[51] J. LIANG AND C.-B. SCHÖNLIEB, *Improving FISTA: Faster, smarter and greedier*, 2018.

[52] N. LORD, J. HO, B. VEMURI, AND S. EISENSCHENK, *Simultaneous registration and parcellation of bilateral hippocampal surface pairs for local asymmetry quantification*, IEEE Trans. Med. Imaging, 26 (2007), pp. 471-478.

[53] F. LUCKA, N. HUYNH, M. BETCKE, E. ZHANG, P. BEARD, B. COX, AND S. ARRIDGE, *Enhancing*
compressed sensing 4d photoacoustic tomography by simultaneous motion estimation, SIAM J. Imaging Sci., 11 (2018), pp. 2224-2253.

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

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

[56] D. MUMFORD AND J. SHAH, *Optimal approximations by piecewise smooth functions and associated variational problems*, Commun. Pure Appl. Anal., 42 (1989), pp. 577-685.

[57] P. NEGRÓN MARRERO, *A numerical method for detecting singular minimizers of multidimensional problems in nonlinear elasticity*, Numer. Math., 58 (1990), pp. 135-144.

[58] F. ODILLE, A. MENINI, J. ESCANYÉ, P. VUISSOZ, P. MARIE, M. BEAUMONT, AND J. FELBLINGER, *Joint Reconstruction of Multiple Images and Motion in MRI: Application to Free-Breathing Myocardial T₂ Quantification*, IEEE Trans. Med. Imag., 35 (2016), pp. 197-207.

[59] S. OZERÉ, C. GOUT, AND C. LE GUYADER, *Joint segmentation/registration model by shape alignment via weighted total variation minimization and nonlinear elasticity*, SIAM J. Imaging Sci., 8 (2015), pp. 1981-2020.

[60] N. PAPADAKIS AND J. RABIN, *Convex histogram-based joint image segmentation with regularized optimal transport cost*, J. Math. Imaging Vis., 59 (2017), pp. 161-186.

[61] R. B. POTTS, *Some generalized order-disorder transformations*, Math. Proc. Camb. Philos. Soc., 48 (1952), p. 106109.

[62] M. RUMPF AND B. WIRTH, *A nonlinear elastic shape averaging approach*, SIAM J. Imaging Sci., 2 (2009), pp. 800-833.

[63] M. RUMPF AND B. WIRTH, *An elasticity-based covariance analysis of shapes*, IJCV, 92 (2011), pp. 281-295.

[64] H. SCHUMACHER, J. MODERSITZKI, AND B. FISCHER, *Combined Reconstruction and Motion Correction in SPECT Imaging*, IEEE Trans. Nucl. Sci., 56 (2009), pp. 73-80.

[65] T. SEDERBERG AND S. PARRY, *Free-form deformation of solid geometric models*, SIGGRAPH Comput. Graph., 20 (1986), pp. 151-160.

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

[67] M. STORATH AND A. WEINMANN, *Fast partitioning of vector-valued images*, SIAM J. Imaging Sci., 7 (2014), pp. 1826-1852.

[68] P. SWIERCZYNSKI, B. W. PAPIEŻ, J. A. SCHNABEL, AND C. MACDONALD, *A level-set approach to joint image segmentation and registration with application to CT lung imaging*, Comput. Med. Imaging Graph., (2017).

[69] C. TOMASI AND T. KANADE, *Shape and motion from image streams under orthography: a factorization method*, Int. J. Comput. Vis., 9 (1992), pp. 137-154.

[70] B. VEMURI, J. YE, Y. CHEN, AND C. LEONARD, *Image registration via level-set motion: applications to atlas-based segmentation*, Med. Image Anal., 7 (2003), pp. 1-20.

[71] T. VERCAUTEREN, X. PENNEC, A. PERCHANT, AND N. AYACHE, *Diffeomorphic Demons: Efficient Non-parametric Image Registration*, NeuroImage, 45 (2008), pp. S61-72.

[72] L. VESE AND C. LE GUYADER, *Variational Methods in Image Processing*, Taylor & Francis, 2015.

[73] M. VON SIEBENTHAL, *Analysis and Modelling of Respiratory Liver Motion using 4DMRI*, PhD thesis, ETH Zurich, 2008.

[74] 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*, Phys. Med. Biol., 52 (2007), pp. 1547-1564.

[75] B. WIRTH, *On the Gamma-limit of joint image segmentation and registration functionals based on phase fields*, Interfaces Free Bound., 18 (2016), pp. 441-477.

[76] P. P. WYATT AND J. NOBLE, *MAP MRF joint segmentation and registration of medical images*, Med. Image Anal., 7 (2003), pp. 539-552.

[77] L. ZAGORCHEV AND A. GOSHTASBY, *A comparative study of transformation functions for nonrigid image registration*, IEEE Trans. Image Process., 15 (2006), pp. 529-538.

[78] Z. ZHANG, Q. XIE, AND A. SRIVASTAVA, *Elastic registration and shape analysis of functional objects*, John Wiley & Sons, Ltd, 2015, ch. 11, pp. 218-238.

[79] H. ZHU, F. MENG, J. CAI, AND S. LU, *Beyond pixels: A comprehensive survey from bottom-up to semantic image segmentation and cosegmentation*, J. Vis. Commun. Image R., 34 (2016), pp. 12-27.