A Time-Dependent Joint Segmentation and Registration Model: Application to Longitudinal Registration of Hepatic DCE-MRI Sequences

While segmentation consists in partitioning a given image into meaningful constituents in order to identify relevant structures such as homogeneous regions or edges, registration, given two images, aims at finding an optimal orientationpreserving one-to-one deformation aligning the structures visible in an image into the corresponding ones in the other. Recently, intertwining both tasks into a single framework has proven to yield better results in terms of accuracy —in particular when the images exhibit weak boundary definition —and increase of reliability of the encoded structure matching —since now, not only based on intensity distribution comparison but also on geometrical and topological features —. In line with this idea, we propose going a step further by adding explicitly some dynamics in the modelling, i.e., by making the minimization problem both space and time-dependent so that the correlation between both tasks is achieved through the process, connecting thus the problem to an interpolation one. The shapes to be matched are viewed as Saint Venant-Kirchhoff materials, a special instance of hyperelastic ones, and are implicitly modelled by level-set functions. These are evolved in order to minimize a functional including both a nonlinear-elasticity-based regularizer prescribing the physical nature of the deformation and a term penalizing the shape misalignment, thus promoting structure matching rather than intensity pairing. Theoretical results emphasizing the mathematical soundness of the model are provided, among which the existence of minimizers and the existence of a weak viscosity solution to the related evolution problem. The model is then applied to the longitudinal registration of hepatic dynamic contrast-enhanced MRI sequences and shows good performance. This application has an important impact on the computer-aided follow-up of patients suffering from liver cancers.


I. INTRODUCTION
Multitask frameworks and especially variational ones have demonstrated significant improvements over sequential approaches. The former cover a large spectrum of image processing problems including combined segmentation and registration models (see references herein below); joint image reconstruction and motion estimation [2], [8], [10], [26]; joint reconstruction and registration for post-acquisition motion correction [11] (with the goal to reconstruct a single motionfree corrected image and retrieve the physiological dynamics through the deformation maps), joint optical flow estimation with phase field segmentation of the flow field [7], or joint segmentation/optimal transport models [5] (to determine the velocity of blood flow in vascular structures). Sharing representation between tasks and carefully intertwining them allows to reduce error propagation, to create synergies, to compensate for some possible flaws such as image quality impairment, while increasing the accuracy of the outcomes and bridging the gap towards generalization. Joint segmentation and registration models such as [16], [18] (joint phase field approximation and registration), [21] (model based on metric structure comparison), [15], [25] (level set formulation that merges the piecewise constant Mumford-Shah model with registration principles), [17] (grounded in the expectation maximization algorithm), [14] (based on a nonlocal characterization of weighted-total variation and nonlocal shape descriptors), or [1], [20], [23], [24], [27], fall within this framework.
Registration can be seen as the incorporation of prior information -such as topological constraints for instance, since the one-to-one property of the deformation makes the moving shape homeomorphic to the original one -in the segmentation process. By the same token, accurate segmented structures enable one to drive the registration process correctly by transferring edges for instance. The primary scope of the paper is thus to go one step further in this bi-task formalism by explicitly introducing some dynamics in the modelling, i.e., by making the optimization problem both space and timedependent phrased then on a Sobolev space of Banach-spacevalued functions. The model promotes large deformations in a nonlinear-elasticity-based setting, and includes a fidelity measure fostering shape overlaps and relying on nonlocal shape descriptors from the piecewise constant Mumford-Shah model in a spatio-temporal setting.
This model is then applied to the longitudinal computeraided diagnosis and follow-up of hepatic diseases. In clinical routine, Dynamic Contrast Enhancement MRI (DCE-MRI) is a standard for liver cancers such as HCC (Hepato-Cellular Carcinoma) for instance [3]. In DCE-MRI four main phases may be visualized: at non-injected phase (phase I), the lesion is difficult to differentiate from healthy tissues; at contrasted arterial phase (phase II), the HCC is enhanced when compared with the surrounding liver tissues; and at the portal and delayed phases (phases III and IV), the signal of the HCC appears slightly lower than the surrounding liver, defined as wash out. In this context, our joint model serves as a registration tool of such sequences and shows good behavior with high Dice coefficient values. As a long-term vision, this will surely improve computerized analysis of DCE-MRI volumes in the follow-up of patients suffering from liver cancers as HCC.

A. Mathematical background
Let Ω be a connected bounded open subset of R 2 of class C 1 . Let us denote by R :Ω → R the Reference image assumed to be sufficiently smooth and by T :Ω → R the Template image. The shape contained in the Template image is assumed to be modelled by a Lipschitz continuous function Φ 0 (input of the problem) whose zero level line is the shape boundary. Denoting by C the zero level set of Φ 0 and by w ⊂ Ω the open set it delineates, Φ 0 is chosen such that For theoretical and numerical purposes, we may consider a linear extension operator (see [6, p. 158 , with C depending only on Ω. By this extension process, we consider then that Φ 0 ∈ W 1,∞ (R 2 ) to ensure that Φ 0 • ϕ -with ϕ introduced later -is always defined. We recall that, in the general case, if U is an open subset of R N , for 1 ≤ p ≤ +∞, the Sobolev space denoted by W 1,p (U ) consists of the functions in L p (U ) whose partial derivatives up to order 1, in the sense of distributions, can be identified with functions in L p (U ). Let ϕ :Ω → R 2 be the sought deformation. A deformation is a smooth mapping that is orientation-preserving and injective, except possibly on ∂Ω. The deformation gradient is ∇ϕ:Ω → M 2 (R), the set M 2 (R) being the set of all real square matrices of order 2 identified to R 4 . This deformation to be recovered is seen as the optimal solution of a specifically designed cost function, comprising a regularization on ϕ prescribing the nature of the deformation, and a term measuring alignment or how the available data are exploited to drive the registration process. These are depicted hereafter.

B. Deformation regularization
Nonlinear elasticity principles dictate the design of the smoothness on ϕ. The shapes to be matched are viewed as hyperelastic materials (capable of undergoing large deformations), and more precisely as Saint-Venant Kirchhoff ones. The considered regularizer is given, setting ξ = ∇ϕ, by : 2 +μ tr E 2 , the stored energy function of a Saint Venant-Kirchhoff material, λ and μ the Lamé coefficients, E = F T F − I /2 the Green-Saint Venant stress tensor measuring the deviation between ϕ and a rigid deformation, and with the following notation A : B = trA T B, the matrix inner product, and A = √ A : A the related matrix norm (Frobenius norm). Additionally, α = 2 λ+μ λ+2μ and Ψ is the convex mapping defined by Ψ : . Several arguments motivate this choice by comparison to the stored energy function W SV K alone: (i) first, although the stored energy function W SV K is the simplest one that agrees with the generic expression of the stored energy function of an isotropic, homogeneous, hyperelastic material, it lacks a term penalizing the determinant. It thus does not preclude deformations with negative Jacobian determinant; (ii) second, it is not rank-1 convex and consequently not quasiconvex (see [13,Chapter 9]), which raises some theoretical issues in terms of existence of minimizers, contrary to QW which exhibits this fine latter property; (iii) thirdly, we see in the expression of QW that when ξ 2 < α, a penalization on the determinant still remains, showing good behavior under compression. By introducing explicitly the time variable t ∈ [0,T ] in the model and the unknown ϕ = ϕ(x, t) being now space and time-dependent, it yields the following regularization on ϕ · 2 denoting the Euclidean norm in R 2 and with V = Ω × [0,T ].

C. Dissimilarity measure
Complying with the idea of promoting structure/shape matching rather than intensity-based rules, we introduce a metric fostering the overlap of shapes as follows. Recall that the shape contained in the Template image is assumed to be modelled by a Lipschitz continuous function Φ 0 whose zero level line is the shape boundary. The idea is to deform the function Φ 0 so that the zero level line of the deformed function, Φ 0 • ϕ, is aligned with the boundary of the shape belonging to the Reference. The trade-off between segmentation and registration -which makes our model a bi-task one -is ensured by processing at the same time the segmentation of the Reference via the level set functionΦ (see 4) and the deformation of Φ 0 , and these tasks are connected in a single framework to maximize the overlap of the two resulting shapes. More precisely, the dissimilarity measure is defined by: H ε denoting a C ∞ -regularization of the one-dimensional Heaviside function and the evolution ofΦ being dictated by the following evolution equation, which constitutes a revised version of [20] to unify local and global (region-based) features to segment the Reference : with (we dropped the dependency onΦ to lighten the expressions). Φ 0 is naturally taken to be the initial condition of this segmentation process. Functiong is an edge-detector function satisfying g(0) = 1,g strictly decreasing and lim r→+∞g (r) = 0. The first two components of the right hand side of equation (4) constitutes a balance between the classical geodesic active contour model ( [9]) and the piecewise constant Mumford-Shah model ( [22]) which allows to partition R into two phases w and Ω \ w with respective values c 1 and c 2 . The latter one prescribes some topological constraints ( [20]) in order for the evolving contour to be homeomorphic to the original one. The overall registration problem thus reads :

D. Theoretical results
We state the main theoretical result related to the existence of minimizers.
The well-definedness ofΦ is then investigated to ensure that the fidelity term exhibits sufficient regularity and makes sense. Equation (4) falls within the framework of viscosity solution theory ( [4], [12]) for equations with a measurable dependence in time (called L 1 -viscosity solution). Under mild assumptions, the following theorem holds. Theorem 2 (Existence of weak solutions of the considered evolution problem): Assuming that g :=g(|∇R|), g 1 2 and ∇g are bounded and Lipschitz continuous on R 2 , for fixed c 1 and c 2 , problem (4) admits at least one weak solution.

III. NUMERICAL RESOLUTION
We now aim to apply this theoretical model to the longitudinal study of an image sequence. Let (I n ) N n=0 be the temporal image sequence. The segmentation evolution problem now reads : For the registration problem, rather than considering a continuum in time which is not realistic, we assume that the problem is sampled in time and drop the regularization in time introduced only for theoretical purposes. We thus solve sequentially the subproblems : for i ∈ {1, . . . , Nζ}, ζ being the number of steps saved in the segmentation process between t n−1 and t n . In the end, the overall deformation is given by ϕ 1 • · · · • ϕ Nζ , and the deformation between the initial frame and the n-th frame by Remark 1: From a theoretical standpoint, the existence of minimizers for each subproblem is guaranteed: Rellich-Kondrachov's embedding theorem states that weak convergence in W 1,4 (Ω, R 2 ) leads to uniform convergence inΩ, an extension process as before can be applied on all ϕ k , k = 1, · · · , i−1 to ensure the well-definedness of the composition, and the continuous injection W 1,4 (R 2 , R 2 ) C 0 (R 2 , R 2 ) holds, these three elements combined allowing to handle the fidelity term. In order to deal with the nonlinearity in ∇ϕ, we propose introducing an auxiliary variable V i simulating the Jacobian deformation with a quadratic penalty method as in [15]. The decoupled problem becomes : for i ∈ {1, . . . , Nζ}. We then use an alternating minimization scheme to solve the problem. We refer the reader to [15,Section 4.3.] for an exhaustive description of the algorithm relying on the derivation of Euler-Lagrange equations solved by an L 2 gradient flow algorithm and an implicit/semi-implicit Euler time stepping scheme.

IV. NUMERICAL SIMULATIONS
The proposed method has been evaluated both qualitatively and quantitatively as a registration process to the longitudinal analysis of 2 DCE-MRI sequences (composed of 4 MRI phases each denoted I, II, III and IV) for a patient suffering from cirrhosis and HCC. DCE-MRI exams were performed on a 1.5T SIGNA™ Artist (General Electric, Milwaukee, WI) with a phased array coil. Manual segmentations of the liver have been elaborated thanks to a Slicer 3D plug-in we are developing [19]. The segmentation step of our method serves here as an interpolation process between the manual segmentations to guide the registration and improve both the matching quality and the topology-preservation property.
To assess the performance of our algorithm, we display the visual outputs in Fig. 1 and use two metrics: the Dice coefficient measuring set agreement with the highest score 1; the mutual information to measure image alignment with larger value meaning better matching. These quantitative measures are reported in Table I along with the max and min value of the deformation Jacobian determinant to appraise the topology preservation. We can see in Fig. 1 that the deformed contour is well aligned with the edge of the liver at all times and the deformation grid does not exhibit overlaps, highlighting the topology preservation property of our algorithm. This is particularly visible in the zoom-in views of the complex topology exhibited by the liver around the vena cava: the segmentation delineates well the concavity without violating topology preservation. It thus shows the capability of our model to handle large deformations on complex shapes and intensity variations between frames. This is further confirmed by the analysis of the Dice coefficients between the deformed initial frame (T • ϕ) and the corresponding temporal frame (R) and their mutual information. Indeed, the latter always significantly increases between the non-deformed frame T and the deformed one T • ϕ, while the Dice coefficient is greater than 98.7% at all times. Favoring shape matching rather than intensity pairing has therefore a positive impact on the registration quality when intensity changes are involved between the temporal frames. Furthermore, we can see in Table I that the determinant remains positive at all times which translates into the ability of our model to generate physically meaningful deformations along the longitudinal analysis. This experiment is thus a proof of concept that our method can be used for the longitudinal analysis of complex organs such as the liver including intensity variations.

V. DISCUSSION AND FUTURE WORKS
In this article, we have proposed a novel joint registrationsegmentation model that minimizes a functional by integrating both space and time-dependent terms. Our contribution enables one to generate large deformations thanks to a nonlinear elasticity regularization process and promotes structure matching rather than intensity-based similarity measuring.
We have successfully exploited our model in a computerized diagnosis and follow-up process by computing the registration of liver shapes within DCE-MRI sequences. By construction, these MRI volumes have variable intensity ranges (due to contrast agent absorption), which would have penalized standard approaches that consider more pixel values than geometrical or topological structures. In this context, we can handle concavities that are generally considered as strong topological constraints for registration. Our experiments show that our model can offer accurate performance by means of excellent Dice coefficient and mutual information measures.
As future works, we first plan to increase the capability of our model by developing a fully 3D approach rather than a slice-based one, and by studying multi-modal joint registration-segmentation through the context of personalized MRI/CT alignment. Moreover, we would like to tackle the problem of internal tissue registration by studying tumors throughout dynamic sequences. In this case, we can extend our framework by considering two level sets: one for the liver volume, and another one for the tumor. We can also produce an approximate registration of the tumors by applying the deformations obtained from the whole liver to internal tissues. Finally, we can exploit further our model by calculating a more precise segmentation of the liver and internal lesions within medical 3D volumes.  1. Visual assessment of the longitudinal analysis performed on the sequence 19-02-2013 Phase I to IV followed by 08-08-2013 Phase I to IV. For each time step, the image with the initial contour along with the image with the deformed contour and the deformation grid are shown. Zoom-in views are provided below to highlight the preservation of topology in complex shapes such as concavities, and the quality of the alignment between the deformed contour and the image edges.