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

Noémie Debroux, Guillaume Lienemann, Benoît Magnin, Carole Le Guyader, Antoine Vacavant

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

## ► To cite this version:

Noémie Debroux, Guillaume Lienemann, Benoît Magnin, Carole Le Guyader, Antoine Vacavant. A Time-Dependent Joint Segmentation and Registration Model: Application to Longitudinal Registration of Hepatic DCE-MRI Sequences. International Conference on Image Processing Theory, Tools and Applications, Nov 2020, Paris, France. hal-03034189

HAL Id: hal-03034189

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

Submitted on 1 Dec 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A Time-Dependent Joint Segmentation and Registration Model: Application to Longitudinal Registration of Hepatic DCE-MRI Sequences

Noémie Debroux

Université Clermont Auvergne
CNRS, SIGMA Clermont, Institut Pascal
F-63000, Clermont-Ferrand, France
noemie.debroux@uca.fr

Guillaume Lienemann

Université Clermont Auvergne
CHU, CNRS, SIGMA Clermont, Institut Pascal
F-63000, Clermont-Ferrand, France
glienemann@chu-clermontferrand.fr

Benoît Magnin

Université Clermont Auvergne
CHU, CNRS, SIGMA Clermont, Institut Pascal
F-63000, Clermont-Ferrand, France
bmagnin@chu-clermontferrand.fr

Carole Le Guyader

Normandie Univ
INSA Rouen, LMI
76000 Rouen, France
carole.le-guyader@insa-rouen.fr

Antoine Vacavant

Université Clermont Auvergne
CNRS, SIGMA Clermont, Institut Pascal
F-63000, Clermont-Ferrand, France
antoine.vacavant@uca.fr

**Abstract**—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 orientation-preserving 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.

**Index Terms**—registration, Saint Venant-Kirchhoff material, topology-preserving segmentation, time-dependent problems, viscosity solutions, hepatic DCE-MRI sequences.

## 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 motion-free 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.

This work was co-funded by the French Agence Nationale de la Recherche (grant ANR-18-CE45-0018, project R-Vessel-X, http://tgi.ip.uca.fr/r-vessel-x), by the European Union with the European Regional Development Fund (ERDF, 18P03390/18E01750/18P02733) and by the Haute-Normandie Regional Council via the M2SINUM project.
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 time-dependent phrased then on a Sobolev space of Banach-space-valued 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 computer-aided 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.

## II. MATHEMATICAL MODELLING

### A. Mathematical background

Let $\Omega$ be a connected bounded open subset of $\mathbb{R}^2$ of class $C^1$. Let us denote by $R: \bar{\Omega} \to \mathbb{R}$ the Reference image assumed to be sufficiently smooth and by $T: \bar{\Omega} \to \mathbb{R}$ the Template image. The shape contained in the Template image is assumed to be modelled by a Lipschitz continuous function $\Phi_0$ (input of the problem) whose zero level line is the shape boundary. Denoting by $C$ the zero level set of $\Phi_0$ and by $w \subset \Omega$ the open set it delineates, $\Phi_0$ is chosen such that $C = \{x \in \Omega | \Phi_0(x) = 0\}$, $w = \{x \in \Omega | \Phi_0(x) > 0\}$ and $\Omega \setminus \bar{w} = \{x \in \Omega | \Phi_0(x) < 0\}$. For theoretical and numerical purposes, we may consider a linear extension operator (see [6, p. 158]) $P: W^{1,\infty}(\Omega) \to W^{1,\infty}(\mathbb{R}^2)$ such that for all $\Phi \in W^{1,\infty}(\Omega)$, (i) $P\Phi|_{\Omega} = \Phi$, (ii) $\|P\Phi\|_{L^{\infty}(\mathbb{R}^2)} \le C \|\Phi\|_{L^{\infty}(\Omega)}$ and (iii) $\|P\Phi\|_{W^{1,\infty}(\mathbb{R}^2)} \le C \|\Phi\|_{W^{1,\infty}(\Omega)}$, with $C$ depending only on $\Omega$. By this extension process, we consider then that $\Phi_0 \in W^{1,\infty}(\mathbb{R}^2)$ to ensure that $\Phi_0 \circ \varphi$ — with $\varphi$ introduced later — is always defined.

We recall that, in the general case, if $U$ is an open subset of $\mathbb{R}^N$, for $1 \le p \le +\infty$, 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 $\varphi: \bar{\Omega} \to \mathbb{R}^2$ be the sought deformation. A deformation is a smooth mapping that is orientation-preserving and injective, except possibly on $\partial\Omega$. The deformation gradient is $\nabla\varphi: \bar{\Omega} \to M_2(\mathbb{R})$, the set $M_2(\mathbb{R})$ being the set of all real square matrices of order 2 identified to $\mathbb{R}^4$. This deformation to be recovered is seen as the optimal solution of a specifically designed cost function, comprising a regularization on $\varphi$ 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 $\varphi$. 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 $\xi = \nabla\varphi$, by:

$$ QW(\xi) = \begin{cases} W_{SVK}(\xi) + \mu (\det \xi - 1)^2 & \text{if } \|\xi\|^2 \ge \alpha, \\ \Psi(\det \xi) & \text{if } \|\xi\|^2 < \alpha, \end{cases} \quad (1) $$

with $W_{SVK}(F) = \frac{\lambda}{2} (\mathrm{tr} E)^2 + \mu \mathrm{tr} E^2$, the stored energy function of a Saint Venant-Kirchhoff material, $\lambda$ and $\mu$ the Lamé coefficients, $E = (F^T F - I)/2$ the Green-Saint Venant stress tensor measuring the deviation between $\varphi$ and a rigid deformation, and with the following notation $A: B = \mathrm{tr} A^T B$, the matrix inner product, and $\|A\| = \sqrt{A}: A$ the related matrix norm (Frobenius norm). Additionally, $\alpha = 2^{\frac{\lambda+\mu}{\lambda+2\mu}}$ and $\Psi$ is the convex mapping defined by $\Psi: s \mapsto -\frac{\mu}{2}s^2 + \mu(s-1)^2 + \frac{\mu(\lambda+\mu)}{2(\lambda+2\mu)}$. Several arguments motivate this choice by $:=\gamma$

comparison to the stored energy function $W_{SVK}$ alone: (i) first, although the stored energy function $W_{SVK}$ 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 $\|\xi\|^2 < \alpha$, a penalization on the determinant still remains, showing good behavior under compression.

By introducing explicitly the time variable $t \in [0, \bar{T}]$ in the model and the unknown $\varphi = \varphi(x, t)$ being now space and time-dependent, it yields the following regularization on $\varphi$

$$ \mathrm{Reg}(\varphi) = \iint_V \|\frac{\partial \varphi}{\partial t}\|_2^2 dx dt + \iint_V QW(\nabla_x \varphi(x, t)) dx dt, \quad (2) $$

$\| \cdot \|_2$ denoting the Euclidean norm in $\mathbb{R}^2$ and with $V = \Omega \times [0, \bar{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 $\Phi_0$ whose zero level line is the shape boundary. The idea is to deform the function $\Phi_0$ so that the zero level line of the deformed function, $\Phi_0 \circ \varphi$, 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 $\tilde{\Phi}$ (see 4) and the deformation of $\Phi_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:

$$ \mathrm{Fid}(\varphi) = \iint_V |H_{\varepsilon}(\Phi_0 \circ \varphi(x, t)) - H_{\varepsilon}(\tilde{\Phi}(x, t))|^2 dx dt, \quad (3) $$

$H_{\varepsilon}$ denoting a $C^{\infty}$-regularization of the one-dimensional Heaviside function and the evolution of $\tilde{\Phi}$ 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 :

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \tilde{\Phi}}{\partial t} = |\nabla \tilde{\Phi}| \left[ \mathrm{div} \left( \tilde{g}(|\nabla R|) \frac{\nabla \tilde{\Phi}}{|\nabla \tilde{\Phi}|} \right) + k \left( (R - c_2)^2 - (R - c_1)^2 \right) \right] + 4 \frac{\mu'}{d^2} \bar{H}(\tilde{\Phi}(x) + l) \bar{H}(l - \tilde{\Phi}(x)) \\ \displaystyle \int_{\Omega} \left[ (x-y, \nabla \tilde{\Phi}(y)) e^{-\|x-y\|^2/d^2} \bar{H}(\tilde{\Phi}(y)+l) \bar{H}(l-\tilde{\Phi}(y)) \right] dy, \\ \tilde{\Phi}(x, 0) = \Phi_0(x), \\ \frac{\partial \tilde{\Phi}}{\partial \vec{n}} = 0, \quad \text{on } \partial \Omega, \end{array} \right. \quad (4) $$

with $c_1 = \int_{\Omega} \bar{H}(\tilde{\Phi}) R(x) dx / \int_{\Omega} \bar{H}(\tilde{\Phi}) dx$ and $c_2 = \int_{\Omega} \bar{H}(-\tilde{\Phi}) R(x) dx / \int_{\Omega} \bar{H}(-\tilde{\Phi}) dx$ (we dropped the dependency on $\tilde{\Phi}$ to lighten the expressions). $\Phi_0$ is naturally taken to be the initial condition of this segmentation process. Function $\tilde{g}$ is an edge-detector function satisfying $\tilde{g}(0) = 1$, $\tilde{g}$ strictly decreasing and $\lim_{r \to +\infty} \tilde{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 $\Omega \setminus 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 :

$$ J(\varphi) = \iint_V \left\| \frac{\partial \varphi}{\partial t}(x, t) \right\|_2^2 dx dt + \iint_V [QW(\nabla_x \varphi(x, t)) + \frac{\nu}{2} (H_{\varepsilon}(\Phi_0 \circ \varphi(x, t)) - H_{\varepsilon}(\tilde{\Phi}(x, t)))^2] dx dt, \quad (5) $$

with $V = \Omega \times [0, \bar{T}]$.

### D. Theoretical results

We state the main theoretical result related to the existence of minimizers.

**Theorem 1 (Existence of minimizers):** There exists at least one minimizer of problem (5) in the functional space $W = \{\varphi \in L^4(0, \bar{T}; W^{1,4}(\Omega, \mathbb{R}^2)) \mid \frac{\partial \varphi}{\partial t} \in L^2(0, \bar{T}; L^2(\Omega, \mathbb{R}^2))\}$ endowed with the norm $\|\varphi\|_W = \|\varphi\|_{L^4(0, \bar{T}; W^{1,4}(\Omega, \mathbb{R}^2))} + \|\frac{\partial \varphi}{\partial t}\|_{L^2(0, \bar{T}; L^2(\Omega, \mathbb{R}^2))}$.

*Proof:* The proof is divided into the following three steps and relies on Aubin-Lions lemma:

* **Coercivity inequality:** By taking $\varphi(t) = \mathrm{Id} \quad \forall t \in [0; \bar{T}]$, we get $J(\varphi) < +\infty$ and the functional is proper. We first derive a coercivity inequality:

$$ \begin{aligned} J(\varphi) \ge & \|\frac{\partial \varphi}{\partial t}\|_{L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2])}^2 \\ & + \frac{\beta}{2(c^4 + 1)} \|\varphi\|_{L^4(0, \bar{T}; W^{1,4}(\Omega, \mathbb{R}^2))}^4 - \beta \alpha^2 \mathrm{meas}(\Omega) \bar{T} \\ & - \frac{\mu(\lambda + 3\mu)}{2(\lambda + 2\mu)} \mathrm{meas}(\Omega) \bar{T} + \gamma \mathrm{meas}(\Omega) \bar{T} + k' \bar{T}. \end{aligned} $$

So, the infimum is finite.

* **Convergence of a minimizing sequence:** Let $\{\varphi^k\}_{k \in \mathbb{N}}$ be a minimizing sequence. As the functional is proper, there exists $\hat{\varphi}$ such that for all $k$ large enough, $J(\varphi^k) \le J(\hat{\varphi}) + 1 < +\infty$. Then from the coercivity inequality, we can deduce that $\{\varphi^k\}_{k \in \mathbb{N}}$ is uniformly bounded in $W$ so we can extract a subsequence still denoted $\{\varphi^k\}$ such that $\varphi^k \to \bar{\varphi}$ in $W$. As $W^{1,4}(\Omega, \mathbb{R}^2) \subset L^2(\Omega, \mathbb{R}^2)$ with compact injection, Aubin-Lions lemma states that the embedding of $W$ in $L^4(0, \bar{T}; L^2(\Omega, \mathbb{R}^2))$ is compact. As $L^4(0, \bar{T}; L^2(\Omega, \mathbb{R}^2)) \subset L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2)$ is continuous then the embedding of $W$ into $L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2)$ is also compact. We can therefore extract a subsequence of $\{\varphi^k\}$ still denoted $\{\varphi^k\}$ such that $\varphi^k \to \bar{\varphi}$ in $L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2])$.

* **Weak lower semi-continuity:** $\|\cdot\|_{L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2])}$ is convex and strongly continuous and so it is weakly lower semi-continuous in $L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2)$ so $\|\frac{\partial \bar{\varphi}}{\partial t}\|_{L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2])} \le \liminf_{k \to +\infty} \|\frac{\partial \varphi^k}{\partial t}\|_{L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2])}$. As, $\varphi^k \to \bar{\varphi}$ in $L^2(0, \bar{T}[\times \Omega, \mathbb{R}^2])$, up to a subsequence, one has pointwise convergence of $\{\varphi^k\}$ to $\bar{\varphi}$ and the dominated convergence theorem enables us to obtain the weak lower semi-continuity of the data fidelity term.

Let $\{\psi_k\}_{k \in \mathbb{N}}$ be a sequence that strongly converges to $\bar{\psi}$ in $L^4(0, \bar{T}; W^{1,4}(\Omega, \mathbb{R}^2))$. Then $\int_{\Omega} \|\bar{\psi}(t) - \psi_k(t)\|_{W^{1,4}(\Omega, \mathbb{R}^2)}^4 dt \underset{k \to +\infty}{\longrightarrow} 0$. By seeing $\|\psi(t) - \psi_k(t)\|_{W^{1,4}(\Omega, \mathbb{R}^2)}^4$ as a real-valued function depending on $t$ defined on $]0, \bar{T}]$ and by applying the reciprocal of the dominated convergence theorem, we get that $\|\bar{\psi}(t) - \psi_k(t)\|_{W^{1,4}(\Omega, \mathbb{R}^2)}^4$ converges to 0 for almost every $t \in ]0, \bar{T}]$ up to a subsequence. So for almost every $t \in ]0, \bar{T}]$, $\psi_k(t)$ strongly converges to $\bar{\psi}(t)$
in $W^{1,4}(\Omega, \mathbb{R}^2)$ and $\det \nabla \psi_k(t) \underset{k \to +\infty}{\to} \det \nabla \bar{\psi}(t)$ in $L^2(\Omega)$. From what was done in the stationary case [15], for almost every $t \in ]0, \bar{T}]$,

$$ \int_{\Omega} QW(\nabla \bar{\psi}(t)) dx \leq \liminf_{k \to +\infty} \int_{\Omega} QW(\nabla \psi_k(t)) dx. $$

So by Fatou's lemma, we get:

$$ \begin{aligned} & \int_0^{\bar{T}} \int_{\Omega} QW(\nabla \bar{\psi}(x, t)) dx dt \\ & \leq \liminf_{k \to +\infty} \int_0^{\bar{T}} \int_{\Omega} QW(\nabla \psi_k(x, t)) dx dt. \end{aligned} $$

Then it is convex and strongly sequentially lower semi-continuous and so it is weakly lower semi-continuous. We finally have

$$ \begin{aligned} & \int_0^{\bar{T}} \int_{\Omega} QW(\nabla \bar{\varphi}(x, t)) dx dt \\ & \leq \liminf_{k \to +\infty} \int_0^{\bar{T}} \int_{\Omega} QW(\nabla \varphi^k(x, t)) dx dt, \end{aligned} $$

which concludes the proof.

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

The well-definedness of $\tilde{\Phi}$ 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 := \tilde{g}(|\nabla R|)$, $g^{1/2}$ and $\nabla g$ are bounded and Lipschitz continuous on $\mathbb{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=0}^N$ be the temporal image sequence. The segmentation evolution problem now reads :

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \tilde{\Phi}}{\partial t} = |\nabla \tilde{\Phi}| \left[ \operatorname{div} \left( \tilde{g}(|\nabla I_n|) \frac{\nabla \tilde{\Phi}}{|\nabla \tilde{\Phi}|} \right) + k((I_n - c_{2,n})^2 - (I_n - c_{1,n})^2) \right] + 4 \frac{\mu'}{d^2} \bar{H}(\tilde{\Phi}(x) + l) \bar{H}(l - \tilde{\Phi}(x)) \\ \displaystyle \int_{\Omega} \left[ \langle x-y, \nabla \tilde{\Phi}(y) \rangle e^{-\frac{|x-y|^2}{d^2}} \bar{H}(\tilde{\Phi}(y)+l) \bar{H}(l-\tilde{\Phi}(y)) \right] dy, \\ c_{1,n} = \frac{\int_{\Omega} \bar{H}(\tilde{\Phi}) I_n(x) dx}{\int_{\Omega} \bar{H}(\tilde{\Phi}) dx}, \quad c_{2,n} = \frac{\int_{\Omega} I_n(x) \bar{H}(-\tilde{\Phi}) dx}{\int_{\Omega} \bar{H}(-\tilde{\Phi}) dx}, \\ \text{for } t \in [t_{n-1}, t_n], N \geq n \geq 1, \\ \tilde{\Phi}(x, 0) = \Phi_0(x), \\ \frac{\partial \tilde{\Phi}}{\partial \bar{n}} = 0, \text{ on } \partial \Omega. \end{array} \right. \quad (6) $$

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 :

$$ \begin{aligned} & \inf_{\varphi_i \in \mathrm{Id} + W_0^{1,4}(\Omega, \mathbb{R}^2)} \int_{\Omega} QW(\nabla \varphi_i) dx \\ & + \frac{\nu}{2} \int_{\Omega} \left( H_{\varepsilon}(\Phi_0 \circ \varphi_1 \circ \cdots \circ \varphi_i) - H_{\varepsilon}(\tilde{\Phi}(\cdot, t_i)) \right)^2 dx, \end{aligned} $$

for $i \in \{1, \dots, N\zeta\}$, $\zeta$ 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 $\varphi_1 \circ \cdots \circ \varphi_{N\zeta}$, and the deformation between the initial frame and the $n$-th frame by $\varphi_1 \circ \cdots \circ \varphi_{n\zeta}$.

**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}(\Omega, \mathbb{R}^2)$ leads to uniform convergence in $\bar{\Omega}$, an extension process as before can be applied on all $\varphi_k$, $k = 1, \dots, i-1$ to ensure the well-definedness of the composition, and the continuous injection $W^{1,4}(\mathbb{R}^2, \mathbb{R}^2) \circ C^0(\mathbb{R}^2, \mathbb{R}^2)$ holds, these three elements combined allowing to handle the fidelity term.

In order to deal with the nonlinearity in $\nabla \varphi$, we propose introducing an auxiliary variable $V_i$ simulating the Jacobian deformation with a quadratic penalty method as in [15]. The decoupled problem becomes :

$$ \begin{aligned} & \inf_{\varphi_i \in \mathrm{Id} + W_0^{1,2}(\Omega, \mathbb{R}^2), V_i \in L^4(\Omega, M_2(\mathbb{R}))} \int_{\Omega} QW(V_i) dx \\ & + \frac{\nu}{2} \int_{\Omega} \left( H_{\varepsilon}(\Phi_0 \circ \varphi_1 \circ \cdots \circ \varphi_i) - H_{\varepsilon}(\tilde{\Phi}(\cdot, t_i)) \right)^2 dx \\ & + \frac{\gamma}{2} \|\nabla \varphi_i - V_i\|_{L^2(\Omega, M_2(\mathbb{R}))}^2, \end{aligned} \quad (7) $$

for $i \in \{1, \dots, N\zeta\}$. 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 o φ) 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 o φ, 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 registration-
segmentation 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 elas-
ticity 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 construc-
tion, 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.

REFERENCES

[1] J.-H. An, Y. Chen, F. Huang, D. Wilson, and E. Geiser. Med. Image Comput. Comput. Assist. Interv., Lecture Notes in Comput. Sci. 3749, chapter A Variational PDE Based Level Set Method for a Simultaneous Segmentation and Non-Rigid Registration, pages 286–293. Springer, Berlin, Heidelberg, 2005.

[2] A. I. Aviles-Rivero, G. Williams, M. J. Graves, and C.-B. Schönlieb. Compressed sensing plus motion (CS+M): A new perspective for improving undersampled MR image reconstruction. Preprint, (https://arxiv.org/abs/1810.10828), 2018.

[3] J. Balogh, D. Victor 3rd, E.H. Asham, S.G. Burroughs, M. Boktour, A. Saharia, X. Li, R.M. Ghobrial, and H.P. Monsour Jr. Hepatocellular carcinoma: a review. *Journal of hepatocellular carcinoma*, 3:41–53, 2016.

[4] G. Barles, P. Cardaliaguet, O. Ley, and A. Monteillet. Existence of weak solutions for general nonlocal and nonlinear second-order parabolic equations. *Nonlinear Analysis: Theory, Methods & Applications*, 71(7-8):2801–2810, 2009.

[5] Y.E. Boink. Combined modelling of optimal transport and segmentation revealing vascular properties. Essay, (http://purl.utwente.nl/essays/70998), 2016.

[6] H. Brezis. Analyse fonctionelle. Théorie et Applications. Dunod, Paris, 2005.

[7] C. Brune, H. Maurer, and M. Wagner. Detection of intensity and motion edges within optical flow via multidimensional control. *SIAM J. Imaging Sci.*, 2(4):1190–1210, 2009.

[8] M. Burger, H. Dirks, L. Frerking, A. Hauptmann, T. Helin, and S. Sil-
tanen. A variational reconstruction method for undersampled dynamic
X-ray tomography based on physical motion models. Inverse Problems,
33(12):124008, 2017.

[9] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic Active Contours. *Int. J. Comput. Vis.*, 22(1):61–87, 1993.

[10] S. Y. Chun and J. A Fessler. Joint Image Reconstruction and Nonrigid Motion Estimation with a Simple Penalty that Encourages Local Invertibility. In *Proceedings of SPIE - The International Society for Optical Engineering*, volume 7258, pages 288–296. SPIE, 2009.

[11] V. Corona, A.I. Aviles-Rivero, N. Debroux, M. Graves, C. Le Guyader, and G. Williams C.-B. Schönlieb and. Multi-tasking to Correct: Motion-Compensated MRI via Joint Reconstruction and Registration. In Scale Space and Variational Methods in Computer Vision, pages 263–274. Springer, Cham, 2019.

[12] M.G. Crandall, H. Ishii, and P.-L. Lions. User's Guide to Viscosity Solutions of Second Order Partial Differential Equations. *Bull. Amer. Math. Soc.*, 27:1–67, 1992.

[13] B. Dacorogna. Direct Methods in the Calculus of Variations, Second Edition. Springer, 2008.

[14] N. Debroux and C. Le Guyader. A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors. *SIAM J. Imaging Sci.*, 11(2):957–990, 2018.

[15] N. Debroux, S. Ozeré, and C. Le Guyader. A non-local topology-preserving segmentation guided registration model. *J. Math. Imaging Vision*, pages 1–24, 2017.

[16] M. Droske and M. Rumpf. Multiscale joint segmentation and registration of image morphology. *IEEE Trans. Pattern Anal. Mach. Intell.*, 29(12):2181–2194, 2007.

[17] A. Gooya, K.M. Pohl, M. Bilello, L. Cirillo, G. Biros, E.R. Melhem, and C. Davatzikos. GLISTR: Glioma Image Segmentation and Registration. *IEEE Trans. Med. Imaging*, 31(10):1941–1954, 2012.

[18] J. Han, B. Berkels, M. Droske, J. Hornegger, M. Rumpf, C. Schaller, J. Scorzin, and H. Urbach. Mumford-Shah Model for One-to-One Edge Matching. *IEEE Trans. Image Process.*, 16(11):2720–2732, 2007.

[19] J. Lamy, T. Pelletier, G. Lienemann, B. Magnin, B. Kerautret, N. Passat, J. Finet, and A. Vacavant. Preliminary results with a new annotation tool for liver volume and inner vessels from DCE-MRI data. In *Virtual Physiological Human*. 2020.
<table><caption>TABLE I
QUANTITATIVE EVALUATION OF THE REGISTRATION IN THE LONGITUDINAL STUDY COMPRISING THE FOUR STEPS OF TWO DCE MRI EXAMS WITH TWO METRICS: DICE COEFFICIENT, MUTUAL INFORMATION AND MIN/MAX OF THE DEFORMATION JACOBIAN DETERMINANT.</caption><thead><tr><th>Considered deformation</th><th>Dice(T,R)</th><th>Dice(R,T○ϕ)</th><th>MI(T,R)</th><th>MI(T○ϕ,R)</th><th>min Det∇ϕ</th><th>max det∇ϕ</th></tr></thead><tbody><tr><td>T (19-02-2013 I) → R (19-02-2013 II)</td><td>0.97687</td><td>0.98084</td><td>0.88958</td><td>0.98288</td><td>0.76866</td><td>1.5092</td></tr><tr><td>T (19-02-2013 I) → R (19-02-2013 III)</td><td>0.91399</td><td>0.98512</td><td>0.57387</td><td>0.95691</td><td>0.59322</td><td>1.9982</td></tr><tr><td>T (19-02-2013 I) → R (19-02-2013 IV)</td><td>0.92772</td><td>0.98417</td><td>0.62573</td><td>0.92490</td><td>0.29685</td><td>2.6096</td></tr><tr><td>T (19-02-2013 I) → R (08-08-2013 I)</td><td>0.92604</td><td>0.99093</td><td>0.64321</td><td>1.1482</td><td>0.44114</td><td>3.6822</td></tr><tr><td>T (19-02-2013 I) → R (08-08-2013 II)</td><td>0.92668</td><td>0.99022</td><td>0.65729</td><td>0.97736</td><td>0.46877</td><td>4.1839</td></tr><tr><td>T (19-02-2013 I) → R (08-08-2013 III)</td><td>0.90736</td><td>0.98647</td><td>0.5604</td><td>0.9057</td><td>0.3663</td><td>5.4014</td></tr><tr><td>T (19-02-2013 I) → R (08-08-2013 IV)</td><td>0.9256</td><td>0.98727</td><td>0.6122</td><td>0.9004</td><td>0.3079</td><td>7.3208</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 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.</figcaption></figure>

[20] C. Le Guyader and L. Vese. A combined segmentation and registration framework with a nonlinear elasticity smoother. *Comput. Vis. Image Underst.*, 115(12):1689-1709, 2011.

[21] N.A. Lord, J. Ho, B.C. Vemuri, and S. Eisenschenk. Simultaneous Registration and Parcellation of Bilateral Hippocampal Surface Pairs for Local Asymmetry Quantification. *IEEE Trans. Med. Imaging*, 26(4):471-478, 2007.

[22] D. Mumford. Optimal approximation by piecewise smooth functions and associated variational problems. *Commun. Pure Appl. Anal.*, pages 577-685, 1989.

[23] S. Ozeré, C. Gout, and C. Le Guyader. Joint Segmentation/Registration Model by Shape Alignment via Weighted Total Variation Minimization and Nonlinear Elasticity. *SIAM J. Imaging Sci.*, 8(3):1981-2020, 2015.

[24] M. Rumpf and B. Wirth. A nonlinear elastic shape averaging approach. *SIAM J. Imaging Sci.*, 2(3):800-833, 2009.

[25] P. Swierczynski, B. W. Papież, J. A. Schnabel, and C. Macdonald. A level-set approach to joint image segmentation and registration with application to CT lung imaging. *Comput. Med. Imaging Graph.*, 65, 2017.

[26] C. Tomasi and T. Kanade. Shape and motion from image streams under orthography: a factorization method. *Int. J. Comput. Vis.*, 9(2):137-154, 1992.

[27] B. Vemuri, J. Ye, Y. Chen, and C. Leonard. Image registration via level-set motion: applications to atlas-based segmentation. *Med. Image Anal.*, 7(1):1-20, 2003.