Noémie Debroux

John Aston
email: j.aston@statslab.cam.ac.uk

Fabien Bonardi
email: fabien.bonardi@univ-evry.fr

Alistair Forbes
email: alistair.forbes@npl.co.uk

Carole Le Guyader
email: carole.le-guyader@insa-rouen.fr

Marina Romanchikova
email: ma-rina.romanchikova@npl.co.uk

Carola-Bibiane Schönlieb

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

Keywords: 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

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 ∞ 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.

1. Introduction. In recent years, joint image processing models have experienced increasing attention, including combined segmentation/registration models [START_REF] Droske | Multiscale joint segmentation and registration of image morphology[END_REF][START_REF] Han | Urbach, Mumford-Shah Model for One-to-One Edge Matching[END_REF] (joint phase field approximation and registration), [START_REF] Lord | Simultaneous registration and parcellation of bilateral hippocampal surface pairs for local asymmetry quantification[END_REF] (model based on metric structure comparison), [START_REF] Debroux | A non-local topology-preserving segmentation guided registration model[END_REF][START_REF] Swierczynski | A level-set approach to joint image segmentation and registration with application to CT lung imaging[END_REF] (level set formulation that merges the piecewise constant Mumford-Shah model with registration principles), [START_REF] Gooya | GLISTR: Glioma Image Segmentation and Registration[END_REF] (grounded in the expectation maximisation algorithm), [START_REF] Debroux | A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors[END_REF] (based on a nonlocal characterisation of weighted-total variation and nonlocal shape descriptors), or [START_REF] An | A Variational PDE Based Level Set Method for a Simultaneous Segmentation and Non-Rigid Registration[END_REF][START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF][START_REF] Ozeré | Joint segmentation/registration model by shape alignment via weighted total variation minimization and nonlinear elasticity[END_REF][START_REF] Rumpf | A nonlinear elastic shape averaging approach[END_REF][START_REF] Vemuri | Image registration via level-set motion: applications to atlas-based segmentation[END_REF][START_REF] Wyatt | MAP MRF joint segmentation and registration of medical images[END_REF]; joint image reconstruction and motion estimation [START_REF] Blume | Joint Reconstruction of Image and Motion in Gated Positron Emission Tomography[END_REF][START_REF] Burger | A variational model for joint motion estimation and image reconstruction[END_REF][START_REF] Chun | Joint Image Reconstruction and Nonrigid Motion Estimation with a Simple Penalty that Encourages Local Invertibility[END_REF][START_REF] Odille | Joint Reconstruction of Multiple Images and Motion in MRI: Application to Free-Breathing Myocardial T2Quantification[END_REF][START_REF] Schumacher | Combined Reconstruction and Motion Correction in SPECT Imaging[END_REF][START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF][START_REF] Burger | A variational reconstruction method for undersampled dynamic X-ray tomography based on physical motion models[END_REF][START_REF] Lucka | Enhancing compressed sensing 4d photoacoustic tomography by simultaneous motion estimation[END_REF][START_REF] Aviles-Rivero | Compressed sensing plus motion (CS+M): A new perspective for improving undersampled MR image reconstruction[END_REF]; joint reconstruction and registration for post-acquisition motion correction [25] with the goal to reconstruct a single motion-free image and retrieve the physiological dynamics through the deformation maps, joint optical flow estimation with phase field segmentation of the flow field [START_REF] Brune | Detection of intensity and motion edges within optical flow via multidimensional control[END_REF], or joint segmentation/optimal transport models [START_REF] Boink | Combined modelling of optimal transport and segmentation revealing vascular properties[END_REF] (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 feature 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 - [START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Guyader | Self-Repelling Snakes for Topology-Preserving Segmentation Models[END_REF][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 accurate post-processing treatments based on mutual analysis: for instance, the representation 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- 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 variableness 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 [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF] (atlas generation problem phrased in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF]), [START_REF] Cazelles | Geodesic PCA versus log-PCA of histograms in the Wasserstein space[END_REF] (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), [START_REF] Zhang | Elastic registration and shape analysis of functional objects[END_REF] (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), [START_REF] Arsigny | A log-Euclidean framework for statistics on diffeomorphisms[END_REF] (statistics performed on the space of diffeomorphisms), [START_REF] Hong | Enforcing local context into shape statistics[END_REF] (the use of a kernel descriptor that characterises local shape properties ensures geometrically meaningful correspondence between shapes with statistical studies of the deformations), [START_REF] Rumpf | A nonlinear elastic shape averaging approach[END_REF][START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF] (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 [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF]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 [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF][START_REF] Modersitzki | FAIR: Flexible Algorithms for Image Registration[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] 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 - [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], [START_REF] Broit | Optimal Registration of Deformed Images[END_REF], [START_REF] Burger | A hyperelastic regularization energy for image registration[END_REF], [START_REF] Christensen | Deformable Shape Models for Anatomy[END_REF], [26], [START_REF] Derfoul | 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[END_REF], [START_REF] Droske | A variational approach to non-rigid morphological registration[END_REF], [START_REF] Fischer | Curvature based image registration[END_REF], [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF], [START_REF] Ozeré | Joint segmentation/registration model by shape alignment via weighted total variation minimization and nonlinear elasticity[END_REF], [START_REF] Rumpf | A nonlinear elastic shape averaging approach[END_REF] -, purely geometric ones - [START_REF] Ashburner | Nonlinear spatial normalization using basis functions[END_REF], [START_REF] Davis | A physics-based coordinate transformation for 3-D image matching[END_REF], [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF], [START_REF] Zagorchev | A comparative study of transformation functions for nonrigid image registration[END_REF] -, models including a priori knowledge ( [START_REF] Clatz | Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion els[END_REF]), 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 [START_REF] Rumpf | A nonlinear elastic shape averaging approach[END_REF], the shapes to be matched are not modeled by their closed contour but through a piecewise-constant partition (Potts model [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Potts | Some generalized order-disorder transformations[END_REF]), which constitutes the main difference with [START_REF] Rumpf | A nonlinear elastic shape averaging approach[END_REF]. 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 [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF], the second one using the Cauchy-stress tensors motivated by [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF], and the last one, more novel and on which the paper focuses, relying on tensor-based smoothing D m splines, influenced by [START_REF] Guyader | Spline approximation of gradient fields: Applications to wind velocity fields[END_REF]. 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 [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF], 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 nonlinear 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.

Depiction of the Model.

Let Ω be a connected bounded open subset of 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 : Ω → R the i-th template image with i = 1, • • • , 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 Ω to ensure that T i • ϕ 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,∞ (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 θ i : Ω → R -the variables {θ i } M i=1 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 θ R : Ω → 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 ϕ i : Ω → R 3 be the sought deformation between θ T i and the unknown mean segmentation θ R . Of course, in practice, ϕ i should be with values in Ω 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 Ω -based on Ball's results [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF]. A deformation is a smooth mapping that is orientation-preserving and injective, except possibly on ∂Ω, if self-contact is allowed. The deformation gradient is ∇ϕ i : Ω → M 3 (R), the set M 3 (R) being the set of real square matrices of order 3. The associated displacement field is denoted by u i such that ϕ i = Id + u i , and ∇ϕ i = I 3 + ∇u i , with Id, the identity mapping, and I 3 , the 3 × 3 identity matrix. We also need the following notations: A : B = trA T B, the matrix inner product and ||A|| = √ A : A, the related matrix norm (Frobenius norm).

Following the joint model philosophy in a variational framework, the sought deformations ϕ i , the segmentations θ T i , for all i = 1, • • • , M , as well as the segmented atlas θ R are seen as minimal arguments of a specifically designed cost function. It comprises a regularisation on ϕ i , for all i = 1, • • • , M , prescribing the nature of the deformations, a penalisation on θ T i , for all i = 1, • • • , M , and θ R , favouring piecewise constant segmentations, a segmentation fidelity term ensuring the closeness of the θ T i to the initial available image T i , and a data-driven term measuring the alignment between the deformed segmentations {θ T i • ϕ i } M i=1 and θ R , intertwining 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 [START_REF] Rumpf | A nonlinear elastic shape averaging approach[END_REF]: "the arithmetic mean x of observations {x i } M i=1 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 ( [START_REF] Ciarlet | Elasticité Tridimensionnelle[END_REF]). 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 ϕ 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 [START_REF] Ciarlet | Elasticité Tridimensionnelle[END_REF][START_REF] Dret | Notes de Cours de DEA. Méthodes mathématiques en élasticité[END_REF]) is given by

W O (F ) = K 1 i=1 a i F γ i + K 2 j=1 b j CofF β j + Γ(detF ), with a i > 0, b j > 0, γ i ≥ 0, β j ≥ 0, for all i = 1, • • • , K 1 and all j = 1, • • • , K 2 , and Γ : ]0, +∞[→ R being a convex function satisfying lim δ→0 + Γ(δ) = lim δ→+∞ Γ(δ) = +∞.
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 ) = a 1 F 4 + a 2 CofF 4 + a 3 (detF -1) 2 + a 4 (detF ) 10 -3(a 1 + a 2 ) -a 4 if detF > 0, +∞ otherwise,
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 [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF] 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 (Ω, M 3 (R)), we complement this stored energy function

W Op by the term 1 { . L ∞ (Ω,M 3 (R)) ≤α} (F ) + 1 { . L ∞ (Ω,M 3 (R)) ≤β} (F -1
), with α ≥ 1, and β ≥ 1, where 1 A denotes the convex characteristic function of a convex set A. Therefore, the regularisation can be written as

W (F ) = Ω W Op (F ) dx + 1 { . {L ∞ (Ω,M 3 (R)) ≤α} (F ) + 1 { . L ∞ (Ω,M 3 (R)) ≤β} (F -1 ). Remark 2.1. In terms of functional spaces, if ϕ ∈ W 1,∞ (Ω, R 3 ) (suitable space owing to the L ∞ hard constraints), Cof ∇ϕ and det ∇ϕ are automatically elements of L ∞ (Ω, M 3 (R))
and L ∞ (Ω) respectively, since L ∞ (Ω, R 3 ) has a structure of commutative Banach algebra. Penalising the L ∞ norm of ∇ϕ thus entails control over the Jacobian determinant. This additional term implicitly gives an upper and lower bound on the Jabobian 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 [START_REF] Lam | Landmark-and Intensity-Based Registration with Large Deformations via Quasiconformal Maps[END_REF][START_REF] Lam | Quasi-Conformal Hybrid Multi-modality Image Registration and its Application to Medical Image Fusion[END_REF][START_REF] Lee | Landmark-Matching Transformation with Large Deformation Via N-dimensional Quasi-conformal Maps[END_REF][START_REF] Chan | Topology-Preserving Image Segmentation by Beltrami Representation of Shapes[END_REF]. 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 ( [START_REF] Lehto | Quasiconformal Mappings in the Plane[END_REF]) (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) ≤ K < ∞, then f will be called K-quasiconformal.'
As from the one hand, K(G) ≥ 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 : Ω → Ω in W 1,2 loc (Ω) is K-quasiconformal. This characterization serves as the core of the models [START_REF] Lam | Landmark-and Intensity-Based Registration with Large Deformations via Quasiconformal Maps[END_REF][START_REF] Lam | Quasi-Conformal Hybrid Multi-modality Image Registration and its Application to Medical Image Fusion[END_REF][START_REF] Lee | Landmark-Matching Transformation with Large Deformation Via N-dimensional Quasi-conformal Maps[END_REF][START_REF] Chan | Topology-Preserving Image Segmentation by Beltrami Representation of Shapes[END_REF] which aim to control the Beltrami coefficient.

Theorem 2.3. [6, Theorem 2.5.4] Suppose f : Ω → Ω is a homeomorphic W 1,2 loc -mapping. Then f is K-quasiconformal if and only if ∂f ∂ z (z) = µ(z) ∂f ∂z (z) for almost everywhere z ∈ Ω,
where µ, called the Beltrami coefficient of f is a bounded measurable function satisfying

µ ∞ ≤ K -1 K + 1 < 1 With ∂f ∂z = 1 2 ( ∂f ∂x -i ∂f ∂y ) and ∂f ∂ z = 1 2 ( ∂f ∂x +i ∂f ∂y ), setting f (z) = f (x+iy) = u(x, y)+iv(x, y)
(so in the context of registration, the sought deformation is ϕ = (u, v) T ), the Jacobian is de-

fined by J f (z) = ∂f ∂z (z) 2 - ∂f ∂ z (z) 2 = ∂u ∂x (x, y) ∂v ∂y (x, y) - ∂u ∂y (x, y) ∂v ∂x (x, y) = det ∇ϕ(x, y). Consequently, if f , W 1,2 loc -homeomorphism, is K-quasiconformal, J f (z) = ∂f ∂z (z) 2 - ∂f ∂ z (z) 2 , = ∂f ∂z (z) 2 1 -|µ(z)| 2 , entailing that J f (z) = det ∇ϕ(x, y) > 0 a.e., since µ ∞ < 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 ( µ ∞ < 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 ( [START_REF] Potts | Some generalized order-disorder transformations[END_REF]). The latter, also known as piecewise-constant Mumford-Shah model [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] with N ∈ N phases/shapes (N is thus a prior), is written, for an observed image f , as

inf u∈U E P otts (u) = N l=1 α 2 T V (u l ) + Ω N l=1 u l (c l -f ) 2 dx, with U = {u = (u l ) l=1,••• ,N ∈ (BV (Ω, {0, 1})) N , N l=1
u l = 1, a.e. on Ω}, and

c l = Ω f u l dx Ω u l dx if Ω u l dx = 0, 0 otherwise,
, α > 0 being a weighting parameter balancing the fidelity term and the regularisation. The notation T V denotes the classical Total Variation, measuring the perimeter length of the set defined by {x ∈ Ω | u l (x) = 1} thanks to the coarea formula ( [START_REF] Demengel | Functional Spaces for the Theory of Elliptic Partial Differential Equations[END_REF][START_REF] Evans | Measure Theory and Fine Properties of Functions[END_REF]). The segmentation/partitioning is then retrieved by ũ = N l=1 u l c l , ũ 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 piecewisesmooth approximation instead of a piecewise-constant approximation (see [START_REF] Vese | Variational Methods in Image Processing[END_REF]) or in terms of histograms (see [START_REF] Papadakis | Convex histogram-based joint image segmentation with regularized optimal transport cost[END_REF]) 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 dist ((θ T i , ϕ i ) i=1,••• ,M , θ R ) = 1 2M M i=1 N l=1 T V (θ T i ,l • ϕ i -θ R,l ),
with notations consistent with the definition of

U, i.e. ∀i ∈ {1, • • • , M }, θ T i = (θ T i ,l ) l=1,••• ,N ∈ (BV (Ω, {0, 1})) N and θ R = (θ R,l ) l=1,••• ,N ∈ (BV (Ω, {0, 1})) N .
Remark 2.5. Consistently with Remark 2.4, we could also envision a model including both the deformations ϕ 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 [START_REF] Vercauteren | Diffeomorphic Demons: Efficient Nonparametric Image Registration[END_REF] 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

inf F 1 (θ R , {θ T i , ϕ i } M i=1 ) = 1 M M i=1 γ T 2 N l=1 T V (θ T i ,l ) + Ω N l=1 θ T i ,l (c T i ,l -T i ) 2 dx + γ R 2 N l=1 T V (θ R,l ) + Ω N l=1 θ R,l (c R,l -T i • ϕ i ) 2 dx + λ 2 N l=1 T V (θ T i ,l • ϕ i -θ R,l ) + W (∇ϕ i ) , (P) with c T i ,l = Ω θ T i ,l (x)T i (x) dx Ω θ T i ,l (x) dx if Ω θ T i ,l (x) dx = 0 0 otherwise , c R,l =    1 M M i=1 Ω θ R,l (x)T i •ϕ i (x) dx Ω θ R,l (x) dx if Ω θ R,l (x) dx = 0 0 otherwise , α ≥ 1 and β ≥ 1.
An illustration of the overall components of the model as well as the pipeline of the resulting analysis is given in Figure 2. 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:

• Ŵ = {ψ ∈ Id + W 1,∞ 0 (Ω, R 3 ), 1 det∇ψ ∈ L 10 (Ω), det ∇ψ > 0 a.e. in Ω, ∇ψ L ∞ (Ω,M 3 (R)) ≤ α, (∇ψ) -1 L ∞ (Ω,M 3 (R)) ≤ β}, The infimum is searched for θ R ∈ U, θ T i ∈ U, and ϕ i ∈ Ŵ for all i ∈ {1, • • • , 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 ∀i ∈ {1, • • • , M }, θ T i is an unknown of the problem, when addressing sub-problem 1 in the algorithm -optimisation over θ T i -, we need to rephrase explicitly the objective function in terms of ϕ -1 i in order to make the computations feasible. It implies that from a mathematical viewpoint, the functional space the ϕ i 's belong to should ensure that the ϕ -1 i 's exist and that the deformations are with values in Ω, which is guaranteed here with the devised stored energy function W Op (F ) and thanks to Ball's results. In particular, ϕ i is a bi-Lipschitz homeomorphism such that det ∇ϕ i > 0 a.e. (ii) According to [1, Theorem 2.6.], if F : U → V is a bi-Lipschitz homeomorphism between open subsets of R N such that, either det Df ≥ 0 H N -a.e. on U or det Df ≤ 0 H N -a.e. on U , ∀f ∈ BV (U ), the function f • F -1 belongs to BV (V ). Our framework fulfills these criteria, so θ T i ,l • ϕ i ∈ BV (Ω) and the vector space structure of BV gives us immediately that for all i ∈ {1,

• • • , M }, θ T i ,l • ϕ i -θ R,l ∈ BV (Ω).
Proof. The proof is based on the theory of the calculus of variations, and relies on Ball's results [START_REF] Ball | Global invertibility of Sobolev functions and the interpenetration of matter[END_REF] and arguments inspired by [START_REF] Wirth | On the Gamma-limit of joint image segmentation and registration functionals based on phase fields[END_REF]. See section SM1 for the detailed proof.

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

Numerical Method of Resolution.

3.1. Description and Analysis of the Numerical Method. Inspired by a prior work by Negrón Marrero [START_REF] Marrero | A numerical method for detecting singular minimizers of multidimensional problems in nonlinear elasticity[END_REF] 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 ϕ from its gradient ∇ϕ 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 ϕ i for each i, W i simulating the inverse Jacobian (∇ϕ i ) -1 for all i, and ∀i ∈ {1,

• • • , M }, ∀l ∈ {1, • • • , N }, θ Ti ,l = θ T i ,l •ϕ i -θ R,l
, to simplify numerical computations, and derive a functional minimisation problem phrased in terms of (θ T i , ϕ i , V i , W i , θ Ti ) i=1,••• ,M , θ R . However, we do not impose equality constraints as in [START_REF] Marrero | A numerical method for detecting singular minimizers of multidimensional problems in nonlinear elasticity[END_REF], 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 F 1,γ ({ϕ i , θ T i , V i , W i } M i=1 , (θ Ti ,l ) i=1,••• ,M l=1,••• ,N , θ R ) = 1 M M i=1 γ T 2 N l=1 T V (θ T i ,l ) + Ω N l=1 θ T i ,l (c T i ,l -T i ) 2 dx + γ R 2 N l=1 T V (θ R,l ) + Ω N l=1 θ R,l (c R,l -T i • ϕ i ) 2 dx + λ 2 N l=1 T V (θ Ti ,l ) + γ Ω N l=1 |θ Ti ,l -(θ T i ,l • ϕ i -θ R,l )| dx + Ω W Op (V i ) dx + γ 4 V i -∇ϕ i 4 L 4 (Ω,M 3 (R)) + 1 { . L ∞ (Ω,M 3 (R)) ≤α} (V i ) + 1 { . L ∞ (Ω,M 3 (R)) ≤β} (W i ) + γ 2 W i -V -1 i 2 L 2 (Ω,M 3 (R)) , (DP) with ∀l ∈ {1, • • • , N }, c T i ,l = Ω θ T i ,l (x)T i (x) dx Ω θ T i ,l (x) dx if Ω θ T i ,l (x) dx = 0 0 otherwise , c R,l =    1 M M i=1 Ω θ R,l (x)T i •ϕ i (x) dx Ω θ R,l (x) dx if Ω θ R,l (x) dx = 0 0 otherwise , α ≥ 1 and β ≥ 1. We address this problem for ϕ i ∈ Id + W 1,4 0 (Ω, R 3 ), V i ∈ {ξ ∈ L ∞ (Ω, M 3 (R)) | detξ > 0 a.e. on Ω, 1 detξ ∈ L 10 (Ω), ξ L ∞ (Ω,M 3 (R)) ≤ α}, W i ∈ {ξ ∈ L 2 (Ω, M 3 (R)) | ξ L ∞ (Ω,M 3 (R)) ≤ β}, θ T i ∈ U such that θ T i ,l • ϕ i ∈ L 1 (Ω) for all l ∈ {1, • • • , N } and for all i ∈ {1, • • • , M }, θ Ti ,l ∈ BV (Ω, {-1, 0, 1}) for all l ∈ {1, • • • , N } and for all i ∈ {1, • • • , M }, and θ R ∈ U.

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 (γ j ) j≥0 be an increasing sequence of positive real numbers such that lim

j→+∞ γ j = +∞. Let ({ϕ i,k j , θ T i ,k j , V i,k j , W i,k j } M i=1 , θ Ti ,l,k j i=1,••• ,M l=1,••• ,N , θ R,k j )
be a minimising sequence of the problem F 1,γ for γ = γ j . Then there exists a subsequence such

that ϕ i,k j W 1,4 (Ω,R 3 ) j→+∞ φi , θ T i ,k j (L 1 (Ω)) N -→ j→+∞ θT i , θ R,k j (L 1 (Ω)) N -→ j→+∞ θR , V i,k j * j→+∞ ∇ φi in L ∞ (Ω, M 3 (R)), W i,k j * j→+∞ (∇ φi ) -1 in L ∞ (Ω, M 3 (R)), θ Ti ,l,k j L 1 (Ω) -→ j→+∞ θT i ,l • φi -θR,l , for all l = 1, • • • , N , for all i = 1, • • • , M , and lim j→+∞ F 1,γ j ({ϕ i,k j , θ T i ,k j , V i,k j , W i,k j } M i=1 , θ Ti ,l,k j i=1,••• ,M l=1,••• ,N , θ R,k j ) = F 1 ( θR , { θT i , φi } M i=1 ) = inf F 1 , so that ( θR , { θT i , φi } M i=1 ) ∈ U M +1 × Ŵ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.

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, θ T i denotes the partition of T i into piecewise constant regions, i.e. θ T i = N l=1 c T i ,l θ T i ,l , the number of shapes, N being an unknown (see Remark 3.3) similarly for θ R -:

inf F 2,γ ({ϕ i , θ T i , V i , θ Ti , W i } M i=1 , θ R ) = 1 M M i=1 γ T ∇θ T i L 0 (Ω) + γ T ∇θ Ti L 0 (Ω) + λ T θ T i -T i 2 L 2 (Ω) + γ R ∇θ R L 0 (Ω) + λ R θ R -T i • ϕ i 2 L 2 (Ω) + γ 2 2 V i -∇ϕ i 2 L 2 (Ω,M 2 (R)) + γ 1 2 θ Ti -(θ T i • ϕ i -θ R ) 2 L 2 (Ω) + Ω W Op (V i , detV i ) dx + 1 { . L ∞ (Ω,M 2 (R)) ≤α} (V i ) + 1 { . L ∞ (Ω,M 2 (R)) ≤β} (W i ) + γ 3 2 W i -V -1 i 2 L 2 (Ω,M 2 (R)) , (DPb) with W Op (ψ, δ) = a 1 ψ 4 + a 2 (δ -1) 2 + a 3 δ 10 -2a 1 -a 3 if δ > 0 +∞ otherwise .
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 [START_REF] Debroux | A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors[END_REF]. Also BV (Ω) → L 2 (Ω) in 2 dimensions, so we can replace the L 1 -penalisation for the auxiliary variable θ Ti by an L 2 -penalisation term.

Remark 3.3. We have also opted for the discrete Potts model for the segmentation as in [START_REF] Storath | Fast partitioning of vector-valued images[END_REF] 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 [START_REF] Chambolle | A convex approach to minimal partitions[END_REF] 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 θ T i . For each i = 1, • • • , M , the problem in θ T i amounts to solve inf

θ T i γ T ∇θ T i L 0 (Ω) + λ T θ T i -T i 2 L 2 (Ω) + γ 1 2 θ Ti -θ T i • ϕ i + θ R 2 L 2 (Ω) , ⇔inf θ T i γ T ∇θ T i L 0 (Ω) + λ T θ T i -T i 2 L 2 (Ω) + γ 1 2 (θ Ti • ϕ -1 i -θ T i + θ R • ϕ -1 i )(det∇ϕ i ) -1 2 2 L 2 (Ω) , ⇔inf θ T i γ T ∇θ T i L 0 (Ω) + λ T + (det∇ϕ i ) -1 γ 1 2 θ T i - λ T T i + γ 1 2 (det∇ϕ i ) -1 (θ Ti • ϕ -1 i + θ R • ϕ -1 i ) λ T + (det∇ϕ i ) -1 γ 1 2 2 L 2 (Ω) .
This amounts to solve the Potts model with nonnegative weights and we use the algorithm in [START_REF] Storath | Fast partitioning of vector-valued images[END_REF] based on the Alternating Direction Method of Multipliers (ADMM) and linear programming.

• Sub-problem 2. Optimisation over θ Ti . For each i = 1, • • • , M , the sub-problem in θ Ti is the following one:

inf θ Ti γ T ∇θ Ti L 0 (Ω) + γ 1 2 θ Ti -(θ T i • ϕ i -θ R ) 2 L 2 (Ω) .
This is again the Potts model and we use the same algorithm [START_REF] Storath | Fast partitioning of vector-valued images[END_REF] to solve it in practice. • Sub-problem 3. Optimisation over θ R . By fixing all the other variables, the optimisation problem with respect to θ R becomes

inf θ R γ R ∇θ R L 0 (Ω) + 1 M M i=1 λ R θ R -T i • ϕ i L 2 (Ω) + γ 1 2 θ R -(θ T i • ϕ i -θ Ti ) 2 L 2 (Ω) , ⇔ inf θ R γ R ∇θ R L 0 (Ω) + (λ R + γ 1 2 )(θ R - ( 1 M M i=1 λ R T i • ϕ i + γ 1 2 (θ T i • ϕ i ) -θ Ti ) λ R + γ 1 2 ) 2 L 2 (Ω) .
This is again a Potts model that we solve with the algorithm [START_REF] Storath | Fast partitioning of vector-valued images[END_REF].

• Sub-problem 4. Optimisation over V i . For each i = 1, • • • , M , the sub-problem in V i reads inf V i F (V i ) + Reg(V i ) = Ω a 1 V i 4 + a 2 (detV i -1) 2 + a 3 (detV i ) 10 dx + γ 2 2 V i -∇ϕ i 2 L 2 (Ω,M 2 (R)) + γ 3 2 W i -V -1 i 2 L 2 (Ω) + 1 { . L ∞ (Ω,M 2 (R)) ≤α} (V i ).
This can be cast as a structured convex non-smooth optimisation problem of the sum of a proper closed convex function Reg(.) = 1 { . L ∞ (Ω,M 2 (R)) ≤α} (.) 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 [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]:

V k+1 i = prox γReg (V k i -γ∇F (V k i )),
with prox γReg (y) = min

x 1 2 x-y 2 2 +γReg(y) = min x 1 2 x-y 2 2 +γ1 { . L ∞ (Ω,M 2 (R)) ≤α} (y) = P { . L ∞ (Ω,M 2 (R)
) ≤α} (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 [START_REF] Liang | Improving FISTA: Faster, smarter and greedier[END_REF].

• Sub-problem 5. Optimisation over W i . For each i = 1, • • • , M , we solve the following minimisation problem inf

W i γ 3 2 W i -V -1 i 2 L 2 (Ω) + 1 { . L ∞ (Ω,M 2 (R)) ≤β} (W i ) = P { . L ∞ (Ω,M 2 (R)) ≤α} (V -1 i ).
• Sub-problem 6. Optimisation over ϕ i .

For each i = 1, • • • , M , the sub-problem in ϕ i reads inf ϕ i γ R θ R -T i • ϕ i 2 L 2 (Ω) + γ 1 2 θ Ti -θ T i • ϕ i + θ R 2 L 2 (Ω) + γ 2 2 V i -∇ϕ i 2 L 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 [START_REF] Laga | Elastic 3D shape analysis using square-root normal field representation[END_REF] 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/mechanicallyoriented and can be viewed as non-straightforward adaptations of [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF]. While in [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF] 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 ∀i ∈ {1, • • • , M }, ϕ i = Id on ∂Ω, no boundary stress is applied, and subsequently, inner volumetric stresses are considered. This constitutes a major difference with the work of [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF] 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 [START_REF] Rumpf | An elasticity-based covariance analysis of shapes[END_REF] : "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 ( [START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF]). 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 (Ω, 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 ϕ -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 -:

min v∈H 3 (Ω,R 2 ) |v| 2 3,Ω,R 2 + γ 2 ρ(∇v + ∇v T -∇u -∇u T -∇u T ∇u) 2 M 2 (R),N + ζ(v -u) 2 R 2 ,N , (4.1)
where

|.| 3,Ω,R 2 is the semi-norm on H 3 (Ω, R 2 ), ζ : H 3 (Ω, R 2 ) → R 2 N v → ζ(v) = (v(a 1 ), • • • , v(a N )) T , ρ : H 2 (Ω, M 2 (R)) → (M 2 (R)) N v → ρ(v) = (v(a 1 ), • • • , v(a N )) T . Also, a 1 , • • • , a N denote the image pixel
coordinates with N the total number of pixels, and

∀ξ ∈ (M 2 (R)) N , ∀η ∈ (M 2 (R)) N , ξ, η M 2 (R),N = N i=1 ξ i : η i , while ∀ξ ∈ (R 2 ) N , ∀η ∈ (R 2 ) N , ξ, η R 2 ,N = N i=1
ξ T i η i . However, stated as it is, the problem is not well defined since u ∈ W 1,∞ 0 (Ω, R 2 ) and does not belong to C 1 (Ω, R 2 ) preventing us from extracting isolated values of ∇u. Therefore, for the theoretical analysis of the model, we introduce

(f k ) ∈ C ∞ 0 (Ω, R 2 ) ∩ W 1,∞ (Ω, R 2
), the sequence from the density result such that,

f k -→ k→+∞ u in W 1,∞ (Ω, 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,••• ,N 0 and A 1 = {b i } i=1,••• ,N 1 be
two sets of N 0 and N 1 points of Ω respectively, containing both a P 1 -unisolvent subset. Let us denote by ρ 0 the operator defined by

ρ 0 : H 3 (Ω, R 2 ) → (R 2 ) N 0 v → ρ 0 (v) = (v(a i )) T i=1,••• ,N 0 ,
and by ρ 1 the operator defined by

ρ 1 : H 2 (Ω, M 2 (R)) → (M 2 (R)) N 1 v → ρ 1 (v) = (v(b i )) T i=1,••• ,N 1 .
We introduce the functionals

F ,k :      H 3 (Ω, R 2 ) → R v → ρ 0 (v -f k ) 2 R 2 ,N 0 + γ 2 ρ 1 (∇v + ∇v T -∇f k -∇f T k -∇f T k ∇f k 2 M 2 (R),N 1 + |v| 2 3,Ω,R 2
, and consider the problem

Search for u ∈ H 3 (Ω, R 2 ) such that: ∀v ∈ H 3 (Ω, R 2 ), F ,k (u ) ≤ F ,k (v) . (4.2)
We omit the explicit dependency of u 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:

   Search for u ∈ H 3 (Ω, R 2 ) such that ∀v ∈ H 3 (Ω, R 2 ), ρ 0 (u ), ρ 0 (v) R 2 ,N 0 + γ 2 ρ 1 (∇u + ∇u T ), ρ 1 (∇v + ∇v T ) M 2 (R),N 1 + (u , v) 3,Ω,R 2 = ρ 0 (v), ρ 0 (f k ) R 2 ,N 0 + γ 2 ρ 1 (∇v + ∇v T ), ρ 1 (∇f k + ∇f T k + ∇f T k ∇f k ) M 2 (R),N 1 . (4.3)
Proof. The detailed proof is available in section SM6.

We now define a new norm equivalent to the classical norm on H 3 (Ω, 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

. A 0 ,3,Ω,R 2 : H 3 (Ω, R 2 ) → R f → f A 0 ,3,Ω,R 2 = ( ρ 0 (f ) 2 R 2 ,N 0 + ρ 1 (∇f + ∇f T ) 2 M 2 (R),N 1 + |f | 2 3,Ω,R 2 ) 1 2
, is a Hilbert norm equivalent to the norm .

3,Ω,R 2 in H 3 (Ω, R 2 ).
Proof. The detailed proof is given in section SM7.

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.

We now focus on a convergence result. Let D be a subset of ]0, +∞[ for which 0 is an accumulation point. For any d ∈ D, let A d be a set of N = N (d) distinct points from Ω that contains a P 1 -unisolvent subset. We assume that sup

ρ d : H 3 (Ω, R 2 ) → (R 2 ) N v → ρ d (v) = ((v(a)) a∈A d ) T ,
and by . A d ,3,Ω,R 2 , the norm defined by

f A d ,3,Ω,R 2 = [ ρ d (f ) 2 N + ρ 1 (∇f + ∇f T ) 2 M 2 (R),N 1 + |f | 3,Ω,R 2 ] 1 2 .
As shown in the previous lemma, .

A d ,3,Ω,R 2 is equivalent to the norm . 3,Ω,R 2 in H 3 (Ω, R 2 ) -but not uniformly in d-. Lemma 4.5. Let B 1 = {b 01 , • • • , b 0,N } be a fixed P 1 -unisolvent subset of Ω. By hypothesis, 0 ∈ D, and lim d→0 sup x∈Ω δ(x, A d ) = 0 holds, so ∀j = 1, • • • , N , ∃(a d 0,j ) d∈D , (∀d ∈ D, a d 0,j ∈ A d ) and b 0j = lim d→0 a d 0j .
For any

d ∈ D, let A d 0 be the set {a d 01 , • • • , a d 0,N } and let . A d 0 ,3,Ω,R 2 be the norm defined by ∀f ∈ H 3 (Ω, R 2 ), f A d 0 ,3,Ω,R 2 = [ N j=1 f (a d 0j ) 2 R 2 + N 1 i=1 ∇f (b i ) + ∇f (b i ) T 2 + |f | 2 3,Ω,R 2 ] 1 2 .
Then there exists η > 0 such that for any d ≤ η, . A d 0 ,3,Ω,R 2 is a norm on H 3 (Ω, R 2 ) uniformly equivalent on D∩]0, η] to the norm . 3,Ω,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: Proof. The detailed proof is available in section SM10.

         Search for u d ∈ H 3 (Ω, R 2 ) such that ∀v ∈ H 3 (Ω, R 2 ), ρ d (u d -f k ) 2 R 2 ,N + γ 2 ρ 1 (∇u d + (∇u d ) T -∇f k -∇f T k -∇f T k ∇f k ) 2 M 2 (R),N 1 + |u d | 2 3,Ω,R 2 ≤ ρ d (v -f k ) 2 R 2 ,N + γ 2 ρ 1 (∇v + (∇v) T -∇f k -∇f T k -∇f T k ∇f k ) 2 M 2 (R),N 1 + |v| 2 3,Ω,R 2 , . (4.4)
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

) i=1,••• ,M j=1,••• ,M = Ω v 1,j v 1,i + v 2,j v 2,i dx,
where (v 1,i v 2,i ) T for i = 1, • • • , M are the displacement fields obtained in a linear space for each image, and the resulting displacement fields are denoted by (v 1,pca,j , v 2,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 θ Ti -θ

T i • ϕ i + θ R 2 L 2 (Ω)
). We do not question here the relevancy of the regularisation on the ϕ i 's since in practice, the ϕ i 's need to be invertible (as in the expression θ T i • ϕ i , θ 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 γ T and γ R weighting the L 0 norm control thus the balance between the L 0 -component and the L 2 -penalisation together with λ R and λ T . High values of γ T /γ R favor few large partitions, while small values yield an approximation exhibiting more jumps. λ R , γ T , and γ 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. σ pca 2 ,i represents the i th resulting Cauchy tensor from the second method and (v 1,pca 2 ,i,δ v 2,pca 2 ,i,δ ) T the associated displacement field obtained with weighting parameter δ 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 θ R • (Id + 5c.10

6 v 1,pca 1 ,i v 2,pca 1 ,i ), θ R • (Id + 5c v 1,pca 2 ,i,δ v 2,pca 2 ,i,δ
), and

θ R • (Id + 50c v 1,pca 3 ,i v 2,pca 3 ,i
), 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.

a 1 a 2 a 3 γ 1 γ 2 γ 3 α β γ R γ T λ T λ R dt nbIter T-shape 1

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-bottomright 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 [START_REF] Von Siebenthal | 4D MR imaging of respiratory organ motion and its variability[END_REF]. 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 ( [START_REF] Siebenthal | Analysis and Modelling of Respiratory Liver Motion using 4DMRI[END_REF]). It exemplifies the ability of the method to alleviate ghosting artefacts, in particular splitting into two effects, through parameter γ 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 the physics ( [START_REF] Siebenthal | Analysis and Modelling of Respiratory Liver Motion using 4DMRI[END_REF]): 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 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.

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. [START_REF] Ambrosio | DC calculus[END_REF][START_REF] Astala | Elliptic Partial Differential Equations and Quasiconformal Mappings in the Plane[END_REF][START_REF] Blume | Joint Reconstruction of Image and Motion in Gated Positron Emission Tomography[END_REF][START_REF] Brune | Detection of intensity and motion edges within optical flow via multidimensional control[END_REF] 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.

Figure 1 .

 1 Figure 1. Illustration of the positive impact of joint approaches against sequential ones.

Figure 2 .

 2 Figure 2. Overview of our framework

  x∈Ω δ(x, A d ) = d, where δ is the Euclidean distance in R 2 . Thus d is the radius of the biggest sphere included in Ω that contains no point from A d . Also d is bounded and lim d→0 sup x∈Ω δ(x, A d ) = 0. For any d ∈ D, let us denote by ρ d the mapping defined by

Theorem 4 . 6 (l = 0 such that u d l l→+∞ f k in H 3 (

 463 Convergence). For any d ∈ D, we denote by u d the unique solution to problem (4.4) for fixed. Then under the above assumptions, there exists a subsequence (u d l ) with lim l→+∞ d Ω, R 2 ), and lim k→+∞ lim l→+∞ u d l -u 1,∞ = 0.

Figure 4 .

 4 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.

Figure 5 .

 5 Figure 5. First four modes of variation obtained with our method based on approximation modelling via contour representations.

Figure 6 .

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

Figure 7 .

 7 Figure 7. Some of the initial slices[START_REF] Ambrosio | DC calculus[END_REF][START_REF] Astala | Elliptic Partial Differential Equations and Quasiconformal Mappings in the Plane[END_REF][START_REF] Blume | Joint Reconstruction of Image and Motion in Gated Positron Emission Tomography[END_REF][START_REF] Brune | Detection of intensity and motion edges within optical flow via multidimensional control[END_REF] 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.

Figure 8 .

 8 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.

Figure 9 .

 9 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.

Figure 11 .

 11 Figure 11. Some of the initial images[START_REF] Ambrosio | DC calculus[END_REF][START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF][START_REF] Astala | Elliptic Partial Differential Equations and Quasiconformal Mappings in the Plane[END_REF][START_REF] Aviles-Rivero | Compressed sensing plus motion (CS+M): A new perspective for improving undersampled MR image reconstruction[END_REF] with the associated segmentations given by our joint algorithm underneath. The last line 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.

Figure 12 .

 12 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.

Figure 13 .

 13 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.

Figure 14 .

 14 Figure 14. Comparison of the first mode of variation for 2 different methods.

Table 1

 1 Parameters.

			5.10 3	0.01	1	8.10 4	1	10	100	3	0.5	1	1	0.001	100
	Heart ED(108)-ES(101)	5	1.10 3	4	1	8.10 4	1	100	100	0.02	0.03	1.5	1.5	0.01	500
	Liver -slice 12	5	1.10 3	4	1.5	8.10 4	1	10	100	0.05	0.05	1.5	1.5	0.01	100
	Execution time		Atlas generation PCA 1 st method	PCA 2 nd method PCA 3 rd method
	T-shape			7 min			2 sec			4 min			44 sec
	Heart ED(108)-ES(101)		49 min			3 sec			10 min		3 min 27 sec
	Liver slice 12			26 min			6 sec			7 min 28 sec		13 min

Table 2

 2 Execution times Figure Input images.

N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCH ÖNLIEB

N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCH ÖNLIEB

N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCH ÖNLIEB

N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCH ÖNLIEB

N. DEBROUX, J. ASTON, F. BONARDI, A. FORBES, C. LE GUYADER, M. ROMANCHIKOVA, AND C.-B. SCH ÖNLIEB

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.