Variational Multi-Task MRI Reconstruction: Joint Reconstruction, Registration and Super-Resolution

TL;DR

This paper introduces a novel variational multi-task framework for MRI that jointly addresses reconstruction, registration, and super-resolution, significantly improving image quality from undersampled, motion-corrupted data.

Contribution

It presents the first integrated multi-task optimization model combining reconstruction, registration, and super-resolution for MRI, outperforming sequential and bi-task approaches.

Findings

Enhanced image sharpness and detail recovery.

Significant motion artifact reduction.

Outperforms state-of-the-art methods.

Abstract

Motion degradation is a central problem in Magnetic Resonance Imaging (MRI). This work addresses the problem of how to obtain higher quality, super-resolved motion-free, reconstructions from highly undersampled MRI data. In this work, we present for the first time a variational multi-task framework that allows joining three relevant tasks in MRI: reconstruction, registration and super-resolution. Our framework takes a set of multiple undersampled MR acquisitions corrupted by motion into a novel multi-task optimisation model, which is composed of an $L^2$ fidelity term that allows sharing representation between tasks, super-resolution foundations and hyperelastic deformations to model biological tissue behaviors. We demonstrate that this combination yields to significant improvements over sequential models and other bi-task methods. Our results exhibit fine details and compensate for motion producing sharp and highly textured images compared to state of the art methods.

Tables

Table 1. TABLE I: Parameter values using for our model and for all datasets. In this table, ”D” stands for Dataset.

	D1	D2	D3	D4	D5
$a_1$	1	1	1	1	1
$a_2$	50	50	50	100	100
${\gamma }_1$	5	5	5	1	1
${\gamma }_2$	${10}^5$	${10}^5$	${10}^5$	${10}^5$	${10}^5$
${\gamma }_3$	$15$	$15$	$15$	$1$	$1$
$\theta$	$5$	$5$	$5$	$5$	$5$
$\sigma$	$1.5$	$1.5$	$1.5$	$2$	$2$
$k$	$2$	$2$	$2$	$2$	$2$
	$500$	$500$	$500$	$500$	$500$
$\alpha$	$0.01$	$0.01$	$0.01$	$0.001$	$0.001$

Supplemental Material for the Article:
Variational Multi-Task MRI Reconstruction: Joint Reconstruction, Registration and Super-Resolution
Veronica Corona, Angelica I. Aviles-Rivero, Noémie Debroux, Carole Le Guyader, Carola-Bibiane Schönlieb

Equations

$$x=\mathcal{F}u+\eta ,$$

$$u^{\ast }\in \arg \underset{u}{\min }\Vert \mathcal{F}u-x{\Vert }_2^2+\delta \Vert \Phi (u){\Vert }_1,$$

$$\underset{u}{\min }\sum _{i=1}^m\Vert \mathrm{DB}{\mathrm{W}}_{\mathrm{i}}u-f_i{\Vert }_2^2+\lambda \mathrm{Reg}(u),$$

$$W_{Op}(F)=\left\{\begin{array}[]{c}a_{1}\|F\|_{F}^{4}+a_{2}\left(\mbox{det}F-\frac{1}{\mbox{det}F}\right)^{4}\text{ if }\mbox{det}F>0,\\ +\infty\text{ otherwise},\end{array}\right.$$

$$\displaystyle E(u,(\phi_{t})_{t=1,\cdots,T})=\frac{1}{T}\underset{t=1}{\overset{T}{\sum}}\delta TV_{g_{t}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(\mathcal{C}u)}\circ\phi_{t}^{-1})$$

$$+\alpha \mathrm{TV}(u)+\frac{1}{2}\Vert \mathcal{F}((\mathcal{C}u\circ {\varphi }_t^{-1}))-x_t{\Vert }_2^2,$$

$$\underset{u,({\varphi }_t{)}_{t=1,\cdots ,T}}{\inf }G(u,({\varphi }_t{)}_{t=1,\cdots ,T})=E(u,({\varphi }_t{)}_{t=1,\cdots ,T})$$

$$+\frac{1}{T}\underset{t=1}{\sum ^T}{\int }_{\Omega }{W}_{Op}(\nabla {\varphi }_t)dx,$$

$$\iff$$

$$+\delta {\mathrm{TV}}_{g_t}((\mathcal{C}u)\circ {\varphi }_t^{-1})+{\int }_{\Omega }{W}_{Op}(\nabla {\varphi }_t)dx,$$

$$\underset{u,{\varphi }_t,z_t,h_t,f_t}{\min }$$

$$+\frac{{\gamma }_3}{2}\Vert \mathcal{F}{h}_t-x_t{\Vert }_2^2+\alpha \mathrm{TV}(u)$$

$$+\frac{{\gamma }_2}{2}\Vert (h_t-(\mathcal{C}u)\circ {\varphi }_t^{-1})\sqrt{\det \nabla ({\varphi }_t{)}^{-1}}{\Vert }_2^2$$

$$+\frac{1}{2\theta }\Vert f_t-h_t{\Vert }_2^2+{\mathrm{TV}}_{g_t}(f_t).$$

$$z_{11}^{k+1}$$

$$\displaystyle-4a_{2}(1+z_{22}^{k})c_{0}c_{1}+\gamma_{1}\frac{\partial v_{t,1}^{k}}{\partial x}\Bigg{)},$$

$$z_{12}^{k+1}$$

$$\displaystyle+4a_{2}z_{21}^{k}c_{0}c_{1}+\gamma_{1}\frac{\partial v_{t,1}^{k}}{\partial y}\Bigg{)},$$

$$z_{21}^{k+1}$$

$$\displaystyle+4a_{2}z_{12}^{k}c_{0}c_{1}+\gamma\frac{\partial v_{t,2}^{k}}{\partial x}\Bigg{)},$$

$$z_{22}^{k+1}$$

$$\displaystyle-4a_{2}(1+z_{11}^{k})c_{0}c_{1}+\gamma_{1}\frac{\partial v_{t,2}^{k}}{\partial y}\Bigg{)},$$

$$0$$

$$+{\gamma }_2(h_t^k\circ {\varphi }_t^k-\mathcal{C}{u}^k)\nabla h_t^k({\varphi }_t^k),$$

$$\displaystyle h_{t}^{k+1}=\mathcal{A}^{*}\Bigg{\{}(\gamma_{2}\operatorname{det}\nabla(\phi_{t}^{-1})^{k+1}\mbox{Id}+\gamma_{3}\mathcal{S}^{*}\mathcal{S}+\frac{1}{\theta}\mbox{Id})^{-1}$$

$$\displaystyle\bigg{(}\mathcal{A}\big{(}\gamma_{2}\operatorname{det}\nabla(\phi_{t}^{-1})^{k+1}(\mathcal{C}u^{k})\circ(\phi_{t}^{-1})^{k+1}+\frac{f_{t}^{k}}{\theta}\big{)}+\gamma_{3}\mathcal{S}^{*}x_{t}\bigg{)}\Bigg{\}}.$$

$$f_t^{n+1}$$

$$p_t^{n+1}$$

$$\underset{u}{\min }\frac{{\gamma }_2}{2T}\sum _i^T\Vert h_t\circ {\varphi }_t-(\mathcal{C}u){\Vert }_2^2+\alpha \mathrm{TV}(u),$$

$$y^{k+1}$$

$$u^{k+1}$$

$$\displaystyle+\frac{\gamma_{2}}{T}\mathcal{C}^{*}\sum_{t=1}^{T}h_{t}^{k+1}\circ\phi_{t}^{k+1}\Big{)}.$$

$$G(u,({\varphi }_t{)}_{t=1,\cdots ,T})\ge \frac{1}{T}\underset{t=1}{\sum ^T}(a_1\Vert \nabla {\varphi }_t{\Vert }_{L^4(\Omega ,M_2(\mathbb{R}))}^4$$

$$+\frac{a_2}{2}\Vert \text{mbox}det\nabla {\varphi }_t{\Vert }_{L^4(\Omega )}^4+\frac{a_2}{2}\Vert \frac{1}{\text{mbox}det\nabla {\varphi }_t}{\Vert }_{L^4(\Omega )}^4$$

$$+\frac{1}{4}\Vert (\mathcal{C}u)\circ {\varphi }_t^{-1}{\Vert }_{L^2(\Omega )}^2+\delta {\mathrm{TV}}_{g_t}((\mathcal{C}u)\circ {\varphi }_t^{-1})+\alpha \mathrm{TV}(u)$$

$$-\frac{1}{2}\Vert x_t{\Vert }_{L^2({\mathbb{R}}^2)}^2-3a_2\text{mbox}meas(\Omega )).$$

$${\int }_{\Omega }\Vert (\nabla {\varphi }_{t,n}{)}^{-1}(x){\Vert }_F^{\tilde{q}}\text{mbox}det\nabla {\varphi }_{t,n}(x)dx$$

$$={\int }_{\Omega }\Vert \nabla {\varphi }_{t,n}(x){\Vert }_F^{\tilde{q}}(\text{mbox}det\nabla {\varphi }_{t,n}{)}^{1-\tilde{q}}dx,$$

$$\le \Vert \nabla {\varphi }_{t,n}{\Vert }_{L^4(\Omega ,M_2(\mathbb{R}))}^{\tilde{q}}{\Vert \frac{1}{\text{mbox}det\nabla {\varphi }_{t,n}}\Vert }_{L^4(\Omega )}^{\tilde{q}-1},$$

$$\Vert (\mathcal{C}{u}_n)\circ {\varphi }_{t,n}^{-1}\circ {\varphi }_{t,n}-{\alpha }_t\circ {\bar{\varphi }}_t{\Vert }_{L^p(\Omega )}$$

$$\le \Vert (\mathcal{C}{u}_n)\circ {\varphi }_{t,n}^{-1}\circ {\varphi }_{t,n}-{\alpha }_t\circ {\varphi }_{t,n}{\Vert }_{L^p(\Omega )}$$

Full text

I Introduction

Magnetic Resonance Imaging (MRI) is a widely used and non-invasive modality that creates detailed images of the anatomical structures of the human body, including undergoing physiological events. It allows radiologists to examine MRI for diagnosis, treatment monitoring and abnormality/disease detection [1]. However, a central limitation of MRI is the prolonged acquisition period needed to reconstruct an image [2]. This constraint is reputed to be a major contributor to image quality degradation, and therefore, compromising the expert interpretation.

Image degradation appears as motion artefacts including blurring effects and geometric distortions [3, 2]. Therefore, the problem of how to reduce the acquisition time whilst producing high quality images, super-resolved and motion-free, is of a great interest in the community, and it is the problem that we address in this paper.

In particular, in a dynamic MRI setting, acquisitions with low signal-to-noise ratios or small anatomical structures might be severely degraded, affecting the final expert’s outcome [4]. These small structures can appear smeared or blurred, and discerning whether these are artefacts or lesions is very challenging for the expert, leading to potential false positive or negative findings [5]. Moreover, movement distortions are most prominent at contrast edges [6].

Although, it is possible to reduce the artefacts by performing breath-holding techniques, there is still residual motion to be compensated. This is mainly produced because the timescale of physiological motion is shorter than the required time to form an image. Likewise, gating strategies [7, 8, 9, 10, 11, 12], which track either the breathing or cardiac cycles, have been also widely explored. However, they are mainly effective for perpetual breathing motion disregarding all other involuntary physiological motion and therefore only partially accurate. Furthermore, it is challenging to precisely co-register these signals to the corresponding MRI data [13].

As an alternative to the aforementioned techniques, a body of research has developed several algorithmic approaches based on the conceptual definition of Compressed Sensing (CS) which has demonstrated promising results since the seminal paper of Lustig et al. [14]. The main idea of using CS is to reconstruct signals from low-dimensional measurements through iterative optimisation relying on sparsity of the image in a transformed domain. Since then, several promising results have been reported in the body of literature e.g. [15], [16], [17], [18], [19]. However, there is still a need for improving the quality of the MRI reconstruction whilst decreasing the number of measurements.

A commonality of previous techniques is that they perform a single task (just reconstruction). However in most recent years, there has been a great interest for improving medical image reconstruction [20, 21, 22] by using what is called multi-tasking models (also known as joint models). The central idea of this perspective is that by sharing representation between tasks and carefully intertwining them, one can create synergies across challenging problems and reduce error propagation, which results in boosting the accuracy of the outcomes whilst achieving better generalisation capabilities.

Following the multi-task perspective, different works have been presented e.g. [20, 21, 23, 22, 24, 25]. Unlike existing approaches from the literature, and to the best of our knowledge, we are presenting for the first time a model that considers more than two tasks (see Fig. 1). In this work, we introduce a new variational multi-tasking framework that integrates, in a single model, three relevant tasks in MRI: reconstruction, registration and super-resolution. Whilst this is a relevant part of this work, our contributions are:

• •

We propose a computationally tractable and mathematically well-motivated variational multi-task framework for motion correction in MRI, in which our novelties largely rely on:

• –

An original optimisation model that is composed of an $L^{2}$ fidelity term that allows sharing representations between three tasks (reconstruction, super-resolution and registration); a weighted total variation (TV) ensuring robustness of our method to intensity changes; a TV regulariser of the highly resolved reconstruction; and a hyperelasticity-based regulariser. We demonstrate that this combination yields to significant improvements over sequential models and existing multi-task methods.

• –

We show that our optimisation problem can be solved efficiently by using auxiliary variables and then splitting it into sub-problems. We show that this requires lower CPU time than several methods from the body of literature.

• •

We extensively evaluate our approach using five datasets and different acceleration factors. We also compare our multi-task framework against existing approaches. Our experiments are further validated by interpretations of experts.

II Related Work

There have been different attempts to improve motion correction in MRI from undersampled data. Besides motion prevention techniques such as breath-holding, another set of of algorithmic approaches has been devoted to correct for motion using image-based motion tracking, where one needs an explicit estimation of the motion in between scans. The predominant scheme, in this context, is image registration which aims at finding a mapping aligning a moving image to a reference one. Following this perspective, the body of literature can be roughly classified into rigid (translations, rotations) and deformable registration.

In the first category, several approaches have been proposed including [26, 27, 28, 29]. However, physiological motion such as cardiac and respiratory, can hardly be characterised by a simple combination of rotations and translations. To mitigate this limitation, motion correction methods based on deformable registration have been proposed such as [30, 31, 32, 33]. However, in a closer look on the aforementioned approaches, a commonality between them is that the algorithmic approaches are performed sequentially. That is - the motion estimation task is executed only after the image reconstruction is computed (from now we refer to this perspective as sequential model). A clear drawback of using this perspective is that the motion estimation highly depends on the quality of the reconstruction as well as on the selection of the reference image.

More recently, a body of research has solved jointly multiple tasks (the so-called multi-task approach) such as image reconstruction and registration in a unified framework. In particular, in the medical domain and following a variational perspective [34], different works have been reported using multi-task approaches. These include SPECT imaging [35, 36], PET [24] and MRI [20, 21, 22]- to name a few. The works with a closer aim to ours are discussed next.

Authors in [25] and [24] proposed a joint model composed of a motion-aware likelihood function and a smoothing term for a simultaneous image reconstruction and motion estimation for PET data. Schumacher et al. [36] presented an algorithmic approach that combines reconstruction and motion correction for SPECT imaging. The authors proposed a variational approach that includes a regulariser penalising an offset of motion parameter - to favour a mean location of the target object. However, the major limitation is that they only consider rigid motions. In the same spirit, authors of  [37, 38] proposed a generic joint reconstruction/registration framework. That model is based on a penalised-likelihood functional, which uses a weighted least square fidelity term along with a spatial and a motion regulariser.

Odille et al. [23] proposed a joint model for MRI image reconstruction and motion estimation. This approach allows for an estimate of both intra and inter-image motion, meaning that, not only the misalignment problem is addressed but also it allows correcting for blurring/ghosting artefacts. More recently in the context of deep-learning (DL), a number of methods has been investigated for image registration - e.g. [39, 40]. Although, certainly, those approaches deserve attention, their review goes beyond the scope of this paper.

III Proposed method

In this section, we introduce our joint variational framework which addresses simultaneously the following three tasks: MRI reconstruction, registration and super-resolution. We introduce the mathematical formulation as separated tasks and then we show how our novel optimisation model judiciously intertwines them. Finally, we describe the numerical realisation of our approach.

Problem statement. We remark to the reader the focus of this work. Given a set of multiple undersampled MR acquisitions $\{x_{t}\}_{t=1}^{T}$ of low resolution and corrupted by motion, we seek to recover a single high resolved, static and motion-corrected image that represents the true underlying anatomy along with the estimation of the breathing dynamics through deformation maps.

III-A Task 1: CS MRI reconstruction

In particular, in standard dynamic MRI, the acquired data is in a time-spatial-frequency space, i.e. k,t-space, which is composed of $x=(x_{m,t})_{m=1,t=1}^{M,T}\in\mathbb{C}^{M\times T}$ measurements. Therefore, the task of MRI reconstruction from those samples, reads:

[TABLE]

where $\mathcal{F}:\mathbb{R}^{N\times T}\to\mathbb{C}^{M\times T}$ is the undersampled MRI forward operator. More precisely, $\mathcal{F}=\mathcal{S}\mathcal{A}$ where $\mathcal{S}$ is a subsampling operator, $\mathcal{A}$ the Fourier operator, and $\mathcal{F}^{*}:\mathbb{C}^{M\times T}\to\mathbb{R}^{N\times T}$ its adjoint. Moreover, $u\in\mathbb{R}^{N\times T}$ is the stack of reconstructed images, $\eta$ an additive Gaussian noise inherent to the acquisition, and $t$ the temporal coordinate.

The MRI reconstruction task is thus highly ill-posed due to the noise and incomplete measurements. However, (1) can be solved by adding prior information and then casting the problem as a CS-based optimisation problem:

[TABLE]

where the first term, i.e. data fidelity term, ensures consistency with the observed data $x$ whilst $\|.\|_{1}$ enforces sparsity in the transformed domain given by $\Phi$, and $\delta$ is a parameter balancing the influence of each term.

In this work, we focus on the Total Variation (TV) [41] regulariser, which, imposing edge sparsity, leads to piecewise constant reconstructions. It has shown great potential since early developments in MRI reconstruction [14]. However one can easily replace this regulariser by any other one in a plug-and-play fashion.

Although, a large body of literature has shown potential results in the context of undersampled MRI reconstruction using CS or its extended philosophies including [16, 42, 18], there is a still room for improvement, and in particular for the problem of reconstructing a single high quality image that reflects the true underlying anatomy. This motivates the use of two more tasks $-$ image registration and super-resolution, which are described next.

III-B Task 2&3: When Image Registration Meets Image Super-Resolution

In a dynamic MRI setting, there are two tasks that show a natural strong correlation: motion estimation and super-resolution. Therefore, our hypothesis is that by unifying these two tasks, one can create synergies leading to error propagation reduction, and therefore, an increase of the image quality.

In a multi-frame variational framework, super-resolution is the problem of restoring a high-resolution image from several low quality images that are corrupted by motion. From a variational perspective, it can be expressed as:

[TABLE]

where $\operatorname{D}$ and $\operatorname{B}$ are the downsampling and blurring operators correspondingly. Moreover, $\operatorname{W_{i}}$ models the geometric warp existing between the observed images $f_{i}$ and the restored image $u$ to correct for motion. Finally, $\operatorname{Reg}(u)$ is a generic regulariser. In this work, the dowsampling operator is modelled as an averaging window, the blurring kernel is assumed to be Gaussian, and the warping operator is viewed as the deformations from a registration task. Whilst for the regulariser we adopt the TV option, our approach is well-suited for the plug-and-play setting. That is- one can easily replace the TV regulariser with other options.

In particular, for our registration method we have the following. Let $\Omega$ be the image domain, i.e. a connected bounded open subset of $\mathbb{R}^{2}$, and $u:\Omega\rightarrow\mathbb{R}$ be the sought single reconstructed image depicting the true underlying anatomy. We introduce the unknown deformations, between the $t$-th acquisition and the image $u$, as $\phi_{t}:\bar{\Omega}\rightarrow\mathbb{R}^{2}$. We remark that the deformations are smooth mappings with topology preserving and injectivity properties. Moreover, let $v_{t}$ be the associated displacements such that $\phi_{t}=\mbox{Id}+v_{t}$, where Id is the identity function. At the practical level, these deformations should be with values in $\bar{\Omega}$, and Ball’s results [43] guarantee this property theoretically for our model. We also consider $\nabla\phi_{t}:\Omega\rightarrow M_{2}(\mathbb{R})$ to be the gradient of the deformation, where $M_{2}(\mathbb{R})$ is the set of real square matrices of order two.

As MRI images biological soft tissues well-modelled by hyperelastic materials, which allows for large and smooth deformations while keeping an elastic behavior, we propose to view the shapes to be matched in the registration process as isotropic, homogeneous and hyperelastic materials of Ogden type. This is reflected in our formulation as a regularisation on the deformations $\phi_{t}$ based on the stored energy function of such a material.

In two dimensions, the stored energy function of an Ogden material, in its general form, is given by the following expression: $W_{O}(F)=\underset{i=1}{\overset{M}{\sum}}a_{i}\|F\|_{F}^{\gamma_{i}}+\Gamma(\mbox{det}F)$, with $a_{i}>0$, $\gamma_{i}\geq 1$ for all $i=1,\cdots,M$ and $\Gamma\,:\,]0;\infty[\rightarrow\mathbb{R}$ a convex function satisfying $\underset{\delta\rightarrow 0^{+}}{\lim}\Gamma(\delta)=\underset{\delta\rightarrow+\infty}{\lim}\Gamma(\delta)=+\infty$, $\|.\|_{F}$ designating the Frobenius matrix norm.

Following [44], we consider the particular energy:

[TABLE]

with $a_{1}>0$, and $a_{2}>0$. Both changes in length and area are penalised and topology preservation is ensured with this formulation.

III-C Variational Multi-Task Model: Reconstruction, Registration and Super-Resolution

In the body of literature, there have been different attempts of using reconstruction, registration and super-resolution. However, they tackled the tasks separately or jointly but up to two tasks. In this part, we describe, for the first time, how these three tasks can be jointly computed to benefit the final reconstruction. The main idea is to exploit temporal redundancy in the data to compensate for motion artefacts due to breathing and/or involuntary movements whilst increasing the resolution to retrieve finer details in the reconstruction. In particular, we now turn to describe how (2) and (3) can be solved in a multi-task framework.

Our variational multi-task framework takes three key factors: firstly the hyperelastic regulariser (4), secondly a discrepancy measure that joins the reconstruction, super-resolution and the registration tasks, and the TV-based regularisers for reconstruction and super-resolution. Moreover, our model accounts for intensity changes, this, by modifying the CS-classical TV regulariser for the weighted TV to enforce edge alignment (see Definition in Section V of the Supplementary Material). From now $\mathcal{F}$ is acting on one single frame. We thus introduce weights $g_{t}$ as the Canny edge detector applied to $G_{\sigma}*\mathcal{F}^{*}x_{t}$ - for each $t=1,\cdots,T$ - where $G_{\sigma}$ is a Gaussian filter of variance $\sigma$.

We thus consider the following fidelity term and regulariser for our high-resolved image:

[TABLE]

where $\mathcal{C}=DB$ comes from the super-resolution formulation. The first term of $F$ seeks to align the edges of the deformed reconstruction $((\mathcal{C}u)\circ\phi_{t}^{-1})$) with the ones of the different acquisitions, whilst regularising it. The second quantity aims to get $\mathcal{F}((\mathcal{C}u)\circ\phi_{t}^{-1})$ close to the acquisitions $x_{t}$, and thus $\mathcal{F}^{*}(x_{t})$ close to $\mathcal{C}u\circ\phi_{t}^{-1}$ to correct for motion.

Our variational multi-task framework is then defined as a combination of (4) and (5), which leads to the following minimisation problem:

[TABLE]

We now introduce the next theorem to set the well-posedness of our model.

Theorem III.1** (Existence of minimisers).**

Let $\mathcal{F}=\mathcal{S}\mathcal{A}:L^{2}(\mathbb{R}^{2})\rightarrow L^{2}(\mathbb{R}^{2})$, $\mathcal{C}:L^{1}(\Omega^{\prime})\rightarrow L^{p}(\Omega)$, be linear bounded and continuous for the strong topology operators with $p\in]1,\frac{8}{5}[$, and $\Omega\subset\Omega^{\prime}$, $\Omega$ and $\Omega^{\prime}$ connected bounded open subsets of $\mathbb{R}^{2}$ with boundaries of class $\mathcal{C}^{1}$ (verified by the chosen operators). With $\delta,\,\alpha,\,a_{1},\,a_{2}>0$, problem (6) admits minimisers $(\bar{u},(\bar{\phi}_{t})_{t=1,\cdots,T})$ on $\mathcal{U}=\{u\in BV(\Omega^{\prime}),\,\phi_{t}\in\mathcal{W},\,\forall t=1,\cdots,T\,|\,(\mathcal{C}u)\circ\phi_{t}^{-1}\in BV_{g_{t},0}(\Omega),\,\forall t\in\{1,\cdots,T\}\}$, with $\mathcal{W}=\{\psi\in\mbox{Id}+W^{1,4}_{0}(\Omega,\mathbb{R}^{2})\,|\,\mbox{det}\nabla\psi\in L^{4}(\Omega),\,\frac{1}{\mbox{det}\nabla\psi}\in L^{4}(\Omega),\,\mbox{det}\nabla\psi>0\text{ a.e. on }\Omega\}$.

Proof.

The proof can be found in Section V of the supplementary material. ∎

In the next section, we detail how the proposed model (6) can be solved in a computational tractable form.

III-D Optimisation Scheme

The numerical realisation of (6) imposes different challenges due to the nonlinearity and nonconvexity in $\nabla\phi_{t}$ and the composition $(\mathcal{C}u)\circ\phi_{t}^{-1}$ in the fidelity term. In this work, we overcome these difficulties by introducing three auxiliary variables $z_{t},\,h_{t},f_{t}$, this, to mimic $\nabla\phi_{t}$, $(\mathcal{C}u)\circ\phi_{t}^{-1}$ and $h_{t}$. We then relax our problem using quadratic penalty terms. This leads to the following discretised decoupled problem:

[TABLE]

We now can solve our minimisation problem by splitting (7) into five more computational tractable sub-problems. We now turn to give more details on each sub-problem.

Sub-problem 1: Optimisation over $z_{t}$. In practice, $z_{t}=(z_{t,1},z_{t,2})^{T}$ simulates the gradient of the displacements $v_{t}=(v_{t,1},v_{t,2})^{T}$ associated to the deformations $\phi_{t}$. For every $z_{t}$, we have $z_{t}=\begin{pmatrix}z_{11}~{}z_{12}\\ z_{21}~{}z_{22}\end{pmatrix}$. For the sake of readability, we drop here the dependency on $t$. We solve the Euler-Lagrange equation with an $L^{2}$ gradient flow and a semi-implicit finite difference scheme and update $z_{t}$ as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with c_{0}=\bigg{(}\operatorname{det}(I+z_{t}^{k})-\frac{1}{\operatorname{det}(I+z_{t}^{k})}\bigg{)}^{3} and $c_{1}=1+\frac{1}{(\operatorname{det}(I+z_{t}^{k}))^{2}}$.

Sub-problem 2: Optimisation over $\phi_{t}$. We solve the Euler-Lagrange equation in $\phi_{t}$, after making the change of variable $y=\phi_{t}^{-1}(x)$ in the $L^{2}$ penalty term, for all $t$, using an $L^{2}$ gradient flow scheme with a semi-implicit Euler time stepping.

[TABLE]

Sub-problem 3: Optimisation over $h_{t}$. The update in $h_{t}$, for all $t$, has a closed form solution using the subsampling operator $\mathcal{S}$ and the Fourier operator $\mathcal{A}$ along with their adjoints $\mathcal{S}^{*}$ and $\mathcal{A}^{*}=\mathcal{A}^{-1}$:

[TABLE]

Sub-problem 4: Optimisation over $f_{t}$. This is solved via Chambolle projection algorithm [45]. For an inner loop over $n=1,\cdots,M$:

[TABLE]

with $\|.\|$ the Euclidean norm. After enough iterations, we set $f_{t}^{k+1}=f_{t}^{n+1}$.

Sub-problem 5: Optimisation over $u$. Finally, using the same change of variables as in sub-problem 2, the problem in $u$ reads:

[TABLE]

and we solve it with a primal-dual algorithm: [46]:

[TABLE]

We remark that, in this work, we solve the registration problem in $z_{t}$ and $\phi_{t}$ in a multi-scale framework from coarser to finer grids and using a regridding technique [47] (see Supplementary Material Section II).

IV Experimental Results

In this section, we present the experimental results performed to validate our proposed approach.

IV-A Data Description

We evaluate our framework on five publicly available datasets.

• •

Dataset 1,2 & 3111https://zenodo.org/record/55345#.XBOkvi2cbUZ. These datasets are 2D T1-weighted data [48] acquired during free breathing of the entire thorax. It was acquired with a 3T Philips Achieva system with matrix size = $215\times 173$, slice thickness=8mm, TR=3.1ms and TE=1.9ms. We remark that each dataset refers to three different patients.

• •

Dataset 4 & 5222http://www.vision.ee.ethz.ch/~organmot/chapter_download.shtml. The datasets are 4DMRI data acquired during free-breathing of the right liver lobe [49]. It was acquired on a 1.5T Philips Achieva system, TR=3.1 ms, coils=4, slices=25, matrix size = $195\times 166$, over roughly one hour on 22 to 30 sagittal slices and a temporal resolution of $2.6-2.8$ Hz.

IV-B Evaluation Protocol

To validate our theory, we expensively evaluate our model as follows.

Comparison against Sequential Models. For the first part of our evaluation, we compared our variational multi-task approach against two well-known models, rigid (RIGID), and hyperelastic (HYPER), for deformations. To run this comparison, we solve the CS reconstruction model with TV, and then register all the frames to a reference frame used as initialisation in our proposed approach. For this, we use the well-established FAIR toolbox [50], where we select rigid and hyperelastic transformations. Finally, we perform the super-resolution task with TV.

Comparison against other Multi-task Approaches. As to the best of our knowledge, this is the first variational model joining three tasks, we compare our model against two models that only joints two tasks- reconstruction and the motion estimation. More precisely, we compared our method against DC-CS [20] and GW-CS [21]. To show robustness and generalisation capabilities of our approach, we ran the comparisons using fully sampled data and acceleration factors = $\{2,4,5,6,8\}$.

Metrics Evaluation. As we seek to recover a single high resolved and motion corrected image, there is not ground truth for this task. Therefore, our evaluation is based on a standard protocol for evaluating MRI reconstruction, that is a user-study (expert scoring). For this, we design a a three-point Likert rating scale in which experts were asked to indicate the level of agreement, ranging from best reconstruction to worst reconstruction. The study is also supported by a nonparametric statistical test. Detailed protocol can be found in Section III of the Supplementary Material. Moreover, to further support our multi-task model, we also offer CPU time comparison against all the compared approaches.

The experiments reported in this section were run under the same conditions in a CPU-based Matlab implementation. We used an Intel core i7 with 4GHz and a 16GB RAM.

IV-C Parameter Selection

In our experiments, we set the parameter of our approach and the compared ones as described next. For our experiments, we set the parameters as displayed in Table I. Whilst for the sequential approaches, based on the FAIR implementation, we set the hyperelastic regularisation parameter for Datasets 1,2 & 3 $=1$ and for Datasets 4 & 5 $=0.1$.

IV-D Results and Discussion

We evaluate our proposed approach following the scheme described in Section IV-B.

$\triangleright$** Is our Multi-tasking Approach Better than a Sequential one?** We start evaluating our approach against two sequential models. We remark to the reader that sequential means to execute tasks (reconstruction registration and super-resolution) one after another. In particular, we compared our approach against two well-known models for deformations: rigid (RIGID) and hyperelastic (HYPER). Results of this comparison are displayed in Figs. 4 and 5, and using different acceleration factors.

In a closer look at those Figs., one can see that our reconstructions have better sharp edges and retrieve fine details, in the heart and below the lung areas, than the sequential approaches. Particularly, the rigid transformation is not able to compensate for the whole motion and thus blurring effects are visible, especially under the lungs, for all acceleration factors. Hyperelastic deformations however, have more degrees of freedom and are capable of better compensating for motion which is manifested as sharp edges in the HYPER reconstruction. Moreover, the darker structure, at the center bottom of the heart, disappears or is much less visible in the HYPER reconstructions than in our approach. This effect is observed for all acceleration factors.

Besides, as the acceleration factor increases, the HYPER reconstruction loses the initial contrast, which is particularly visible for the acceleration factor of 8. In contrast our multi-task framework is able to preserve it nicely. This shows the robustness of our method to noise and corrupted data. The benefits of our multi-task framework is prevalent to all datasets (see Section IV of the Supplementary Material for Dataset 2, 4 and 5).

Figs. 4 and 5 show that hyperelastic deformations are better suited to deal with complex physiological motions, as the RIGID reconstructions exhibits strong blurring artefacts, this, due to residual movements amplified as the acceleration factor increases. Also, our method is able to preserve small structures in the kidney and the white blood vessels in the liver even for large acceleration factors contrary to the sequential HYPER approach. For the acceleration factor of 8, the HYPER reconstruction suffers more from staircaising effects than our approach and loses the initial contrast.

Overall, we can show that sharing representation between tasks (i.e. our multi-task approach) leads to better MRI reconstructions than if one performs the task separately. This is strongly supported by two factors, the computational time and the expert agreement. Following common protocol for MRI evaluation, we performed a user-study, in which we asked to twelve experts (radiologist trainees and experienced) to evaluate reconstructions with all acceleration factors in Figs. 4 and 5 (see Supplementary material for Datasets 2, 4 and 5).

The outcome is displayed in Fig. 3. At left side of this figure, one can see that overall (i.e. for all reconstruction/all acceleration factors) our approach was ranked best, with a $44.29\%$ of agreement, in comparison with the output from the other methods. We also ran the nonparametric Friedman test, per acceleration and therefore accounting for FDR, and we found that there is significant statistical difference- that is, our approach offered the best reconstructions.

To further support our model performance, we also analyse the difference maps to assess the quality of our registration, and therefore, its motion correction potential. To do this, we inspect the uncorrected average of the difference image, between a reference frame and each individual one, which is displayed at the left side of Fig. 8. From this column, we can observe that the motion is significant in both the datasets. However, when we inspect the mean difference between our reconstruction and the individual registered acquisitions ($h_{t}\circ\phi_{t}$), at the middle and right sides of Fig. 8, one can see that the structures are very well-aligned resulting in a much smaller range in difference maps. Overall, our approach successfully corrects for motion even at low undersampling rates, and this effect is preserved for all datasets.

$\triangleright$** Is it Three-Task Better Than Two-Task Framework?** In a multi-task framework a key factor is to assess if the tasks are not affecting negatively the final MRI reconstruction. To evaluate this factor, we ran a set of experiments of our approach against two multi-task frameworks DC-CS [20] and GW-CS [21]. These approaches perform only two-tasks (reconstruction and motion estimation), these approaches are our baselines as, to the best of our knowledge, there exists no approaches that join three tasks.

The MRI reconstruction from our model against DC-CS and GW-CS can be seen in Figs. 6 and 7 (see Supplementary Material Section IV for further visualisation with the remaining datasets). In a closer look at these figures, one can observe very blurred reconstructions from DC-CS, which can be interpreted as a failure of the model to capture the complex intrinsic nature of physiological motions. In contrast, GW-CS and our reconstructions are sharper even for very low undersampling factors and compensate well for motion.

However, our method is more robust to noise and outliers (as displayed in the compared reconstructions). Although, the GW-CS reconstructions preserve fine textures and small structures, they are noisier than ours. That is, our approach preserve improves in terms of preserving information whilst removing noise in comparison with GW-CS. This effect is elevated even more as the acceleration factor increases. For example, for an acceleration factor of 8, artefacts and noise are visible in the heart and under the lungs in the GW-CS reconstruction whereas ours is clearer.

Another example of the good performance of our approach can be seen in Fig.7, in which we are able to retrieve more clinically useful texture and fine details than the GW-CS technique. This is particularly visible in the central part of the heart where noise is visible in the GW-CS reconstructions especially as the acceleration factor increases.

To further support our results, we display, at the right side of Fig. 3, the overall outcome of the user-study. From this plot, we can see that the majority of the expert agreed that our reconstructions are better than the compared approaches. Although, the second best ranked is GW-CS, it fails to correct for noise which compromises the readability of the underlying texture. Moreover, as soon as the acceleration factor increases, the noise level jumps, reducing drastically the readability and interpretability of the GW-CS reconstructions whereas our method retrieves relevant small structures and denoises the reconstruction.

$\triangleright$** The CPU Cost of Our Multi-Tasking Approach - Does It Pay Off?** From previous sections, we demonstrated that our approach achieves a better reconstruction in comparison with other approaches, however, does this improvement come to pay off in terms of computational time? Therefore, in this section, we highlight the computational advantages of our model. We remark to the reader that all comparisons were run under the same conditions.

The CPU time, for all approaches, is displayed at Fig. 2 for several acceleration factors (the remaining can be found in Section IV of the Supplementary Material). Firstly, we observe that, in terms of sequential models, our model outperforms RIGID and HYPER reporting the lower CPU time. However and in terms of the other multi-task methods, the computational time for GW-CS is much longer compared to our method (and only performing two-tasks). Whilst the DC-CS approach readily competes with our approach, from the image quality standpoint our method offers by far better results in terms of reconstruction. We emphasise that the CPU times for our model and DC-CS are still on the same range, but the proposed method is performing three-tasks instead of two like the DC-CS. These advantages highlight our optimisation scheme that allows computing a complex problem in a very computational tractable form.

V Conclusion

In this work, we proposed a novel variational multi-task framework to achieve higher quality and super-resolved reconstructions. Our method compensates for motion due to breathing in undersampled data. To the best of our knowledge, it is the first variational framework that allows computing three tasks jointly.

In particular, our multi-task framework is composed of four major components: an $L^{2}$ fidelity term intertwining MRI reconstruction, super-resolution and registration; a weighted TV ensuring robustness of our method to intensity changes by promoting edge alignment; a TV regulariser of the super-resolved reconstruction; and a hyperelasticity-based regulariser modelling biological tissue behavior and allowing for large and smooth deformations. We exploit the temporal redundancy to correct for blurring artefacts and increase image quality. As a result, we obtain a single highly resolved and clear image reconstruction representing the true underlying anatomy.

The advantages of our model is that we guarantee preservation of anatomical structures whilst keeping fine details and less blurry and noise artefacts in the final reconstructions. We extensively evaluated our method against sequential and another multi-task methods from the body of literature. We demonstrated that our method achieves the best results whilst demanding low CPU time. Our method was further supported by a user-study (experts).

Future Work. This multi-task framework is indeed very well-suited for the plug-and-play setting when one (or more) imaging tasks could be replaced by different algorithms. For instance, the modelling of the regularisation functional for the high resolution image reconstruction could be replaced in a plug-and-play fashion. It opens the door to hybrid methods as deep-learning can be used for the super-resolution task.

Acknowledgment

VC acknowledges the financial support of the Cambridge Cancer Centre and the Cambridge Research UK. Support from the Centre for Mathematical Imaging in Healthcare (CMIH), Cantab Capital Institute for the Mathematics of Information (CCIMI) and the National Physical Laboratory (NPL) are greatly acknowledged.

I Outline

This supplementary material extends the experimental results, practicalities and theoretical analysis to further support our proposed multi-task framework. This document is organised as follows.

• •

Section II. In this section, we offer a detailed explanation about the influence of the parameters of out variational multi-task model.

• •

Section III. We further describe the user-study protocol that we follow to support our outcomes and comparisons.

---

• •

Section IV. We display further visual comparisons of our model against some works from the state-of-the-art. This section further validate our claim in terms of our model performance.

• •

Section V. We provide further mathematical details for the paper. First, we recall the definition of the weighted total variation, then we introduce the proof, of Theorem III.1 regarding the existence of minimisers result.

II Parameters Reasoning of Our Multi-Task Framework

In this section, we discuss the influence of each parameter. $a_{1}$ and $a_{2}$ control the regularisation of the deformations. Whilst the former acts on the smoothness of the deformations, the latter can be seen as a measure of rigidity. That is- the bigger $a_{2}$ is, the more rigid the deformations are and the less accurate the registration becomes. It thus behaves as a trade-off between the ability to handle large and nonlinear deformations, and topology preservation.

Moreover, $\gamma_{1}$ and $\gamma_{2}$ are chosen big, to ensure the closeness of the auxiliary variables to the original ones, and $\theta$ is set small for the same purpose. $\gamma_{3}$ weights the fidelity term joining the three tasks, and it is often chosen to be close to 1. Finally, $\alpha$ offers a balance between regularity and fidelity to the data for the super-resolved reconstructed image $u$.

Finally, we also explicitly define he regridding algorithm that we used in Section III.D from the main paper. This is displayed at Algorithm 1.

III Further Details on the User-Study

To further support our results, we performed a user-study following standard protocol in clinical settings. That is- to ask experts and trainees to evaluate the reconstructions based on a scoring system. For this study, we had twelve experts and trainees from the area of radiology.

Scoring Procedure. We create an electronic survey in which, after giving the instructions to the users, they were provided with two-part evaluation. For the first part, they evaluated the reconstructions related to sequential models. To do this, they were provided with the reconstructions of the five dataset, for acceleration factor ={fully sampled, 2x, 4x, 6x, 8x}, and for the approaches ={RIGID, HYPER, OURS} (see main paper for description on these models). With the purpose of capturing they scores, we design a three-point Likert rating scale in which experts were asked to indicate the level of agreement, ranging from best reconstruction to worst reconstruction. Similarly, this protocol was followed to evaluate other multi-task methods = {DC-CD, GW-CS, OURS}. Fig. 9 displays the results in terms of % the worst ranked for all compared approaches.

Statistical Analysis. The circle plots displayed in the main paper and in Fig. 9 of this supplementary material, reflects the averaged results of all scoring (i.e. for all acceleration and all reconstructions). However, to test if there was statistical significant difference between the approaches we took into account the the scoring per each reconstruction, and ran the nonparametric Friedman test to compare the three approaches (for both sequential and multi-task frameworks). We also applied corrections for multiple comparisons during the statistical analysis. Overall, we found statistical significant difference between our approach and the compared ones. This supports our previous discussions on our model offering the best reconstruction in comparison to the compared approached for both sequential and joint models.

IV Further Visual Experiments

In this section, we extend the reconstruction comparisons from the main paper. In particular, for the sequential model we now included Dataset 2, 4 and 5. The results are displayed in Figs. 12, 13 and 14 in which we observe that our joint approach successfully corrects for motion and aligns the different acquisitions resulting in a sharp reconstruction that contains very fine details, such as preserving the texture in the heart and vessels in the liver and kidney. It further support the finding from the main paper.

Moreover, we also provide further evaluation of our framework against other multi-task approaches (DC-CS [20] and GW-CS [21]). The reconstruction results are displayed at Figs. 15, 16 and 17. In a closer look at those reconstructions, one can observe that our approach produces very sharp reconstructions while preserving the texture in the relevant anatomical areas. This is validated for different acceleration factors. In contrast, DC-CS fails to carefully align the acquisition as we can see from the blurring artefacts around the edges. Additionally all the small details get blurred and smeared out. The GW-CS seems to perform reasonably well for fully-sampled data, however failing to remove the noise. This is accentuated at higher reduction factors, where the reconstructions contain significant noise which compromises the clinical interpretation of the images.

Moreover, we present for two datasets the estimated motion $\phi_{t}$ and its inverse $\phi_{t}^{-1}$ for a given time frame in Fig. 11. We can see that our proposed approach produces a reasonable estimation of the motion, where the motion fields are visualised by a deformation grid. Additionally, we show the corresponding Jacobian determinant maps for the deformations. In these plots, we can see that the determinants remain positive meaning that our model ensures topology preservation both from a mathematical and practical point of view. The values are interpreted as follows: small deformations when values are closer to 1, big expansions when values are greater than 1, and big contractions when values are smaller than 1. Moreover, one can observe that the determinants remain positive in all cases, that is to say, our estimated deformations are physically meaningful and preserve the topology as required in a registration framework.

Finally, we also report in Fig. 10, the CPU times for the remaining accelaretion factors. Again, we can see that our approach allows computing a complex problem in a very computationally tractable way. It is indeed much faster than sequential approaches and than the GW-CS. It is comparable to the DC-CS approach, which, however, only computes two tasks.

V Further Mathematical Details

We recall the definition of the weighted BV-space and the associated weighted TV.

Definition 1** ([51, Definition 2]).**

Let $w$ be a weight function satisfying some properties (defined in [51] and fulfilled by $g_{t}$). We denote by $BV_{w}(\Omega)$ the set of functions $f\in L^{1}(\Omega,w)$, which are integrable with respect to the measure $w(x)dx$, such that: $\sup\left\{\int_{\Omega}f\mbox{div}(\varphi)\,dx\,:\,|\varphi|\leq w\text{ everywhere},\,\\ \varphi\in Lip_{0}(\Omega,\mathbb{R}^{2})\right\}$ $<+\infty,$ where $Lip_{0}(\Omega,\mathbb{R}^{2})$ is the space of Lipschitz continuous functions with compact support. We denote by $TV_{w}$ the previous quantity.

V-A Detailed Proof of Theorem III.1

We report the detailed proof of the Theorem III.1, from the main paper , regarding the existence of minimisers. From Section III.C from the main paper.

Proof.

The proof is based on arguments coming from the calculus of variations and is divided into three steps. We recall the assumptions on $g_{t}$: $g_{t}:\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$, $g_{t}(0)=1$, $g_{t}$ is strictly decreasing, $\underset{r\rightarrow+\infty}{\lim}\,g(r)=0$, and there exists $c>0$ such that $c\leq g_{t}\leq 1$ everywhere. We also have $\phi_{t}:\bar{\Omega}\rightarrow\bar{\Omega}$ thanks to Ball’s results [43] as seen later.

Coercivity inequality: We first have that $G(0,(\mbox{Id})_{t=1,\cdots,T})=4a_{1}\mbox{meas}(\Omega)+\frac{1}{2T}\underset{t=1}{\overset{T}{\sum}}\|x_{t}\|_{L^{2}(\mathbb{R}^{2})}^{2}<+\infty$. Let $u\in BV(\Omega^{\prime})$, $\phi_{t}\in\mathcal{W},\,\forall t\in\{1,\cdots,T\}$ such that $(\mathcal{C}u)\circ\phi_{t}\in BV_{g_{t},0}(\Omega),\forall t\in\{1,\cdots,T\}$, we then derive a coercivity inequality:

[TABLE]

Indeed, $((\mathcal{C}u)\circ\phi_{t}^{-1})\in BV_{g_{t},0}(\Omega)\subset BV(\Omega)\subset L^{2}(\Omega)$ ([51]) with continuous embedding, and $((\mathcal{C}u)\circ\phi_{t}^{-1})_{e}$ is the extension by [math] outside the domain $\Omega$ of $(\mathcal{C}u)\circ\phi_{t}^{-1}$, then $\|((\mathcal{C}u)\circ\phi_{t}^{-1})_{e}\|_{L^{2}(\mathbb{R}^{2})}=\|(\mathcal{C}u)\circ\phi_{t}^{-1}\|_{L^{2}(\Omega)}<+\infty$, and finally Plancherel’s theorem gives us $\|\mathcal{F}((\mathcal{C}u)\circ\phi_{t}^{-1})_{e}\|_{L^{2}(\mathbb{R}^{2})}^{2}\leq\|\mathcal{A}((\mathcal{C}u)\circ\phi_{t}^{-1})_{e}\|_{L^{2}(\mathbb{R}^{2})}^{2}=\|((\mathcal{C}u)\circ\phi_{t}^{-1})_{e}\|_{L^{2}(\mathbb{R}^{2})}^{2}=\|(\mathcal{C}u)\circ\phi_{t}^{-1}\|_{L^{2}(\Omega)}^{2}$.

The infimum is thus finite. Let $(u_{n},\,(\phi_{t,n})_{t=1,\cdots,T})_{n}\in\{u\in BV(\Omega^{\prime}),\,\phi_{t}\in\mathcal{W},\,\forall t=1,\cdots,T\,|\,(\mathcal{C}u)\circ\phi_{t}^{-1}\in BV_{g_{t},0}(\Omega),\,\forall t\in\{1,\cdots,T\}\}$ be a minimizing sequence such that $\underset{n\rightarrow+\infty}{\lim}G(u_{n},(\phi_{t,n})_{t=1,\cdots,N})=\underset{(u,(\phi_{t})_{t=1,\cdots,T})\in\mathcal{U}}{\inf}G(u,(\phi_{t})_{t=1,\cdots,T})$.

Extraction of converging subsequences: From the previous coercivity inequality and the finiteness of the infimum we deduce that:

• •

$(\phi_{t,n})_{n}$ is uniformly bounded according to $n$ in $W^{1,4}(\Omega,\mathbb{R}^{2})$ for each $t=1,\cdots,T$ by using the generalized Poincaré’s inequality and the fact that $\phi_{t,n}$ is equal to the identity on the boundary $\partial\Omega$.

• •

$(\mbox{det}\nabla\phi_{t,n})_{n}$ is uniformly bounded according to $n$ in $L^{4}(\Omega)$ for each $t=1,\cdots,T$.

• •

$((\mathcal{C}u_{n})\circ\phi_{t,n})_{n}$ is uniformly bounded according to $n$ in $BV_{g_{t}}(\Omega)$ and thus in $BV(\Omega)$ since $c\leq g_{t}\leq 1$ everywhere for each $t=1,\cdots,T$.

Therefore, we can extract subsequences (but for the sake of conciseness we keep the same notations) such that:

• •

$\phi_{t,n}\underset{n\rightarrow+\infty}{\rightharpoonup}\bar{\phi}_{t}$ in $W^{1,4}(\Omega,\mathbb{R}^{2})$ for each $t=1,\cdots,T$ and by continuity of the trace operator, we deduce that $\bar{\phi}_{t}\in\mbox{Id}+W^{1,4}_{0}(\Omega,\mathbb{R}^{2})$.

• •

$\mbox{det}\nabla\phi_{t,n}\underset{n\rightarrow+\infty}{\rightharpoonup}\delta_{t}$ in $L^{4}(\Omega)$ for all $t=1,\cdots,T$.

• •

$((\mathcal{C}u_{n})\circ\phi_{t,n}^{-1})\underset{n\rightarrow+\infty}{\longrightarrow}\alpha_{t}$ in $L^{q}(\Omega)$ for $q\in[1,2[$ and thus in $L^{1}(\Omega,g_{t})$ with $\alpha_{t}\in BV(\Omega)\subset BV_{g_{t}}(\Omega)$. By continuity of the trace operator, we deduce that $\alpha_{t}\in BV_{g_{t},0}(\Omega)$.

• •

$((\mathcal{C}u_{n})\circ\phi_{t,n}^{-1})\underset{n\rightarrow+\infty}{\rightharpoonup}\alpha_{t}$ in $L^{2}(\Omega)$ for each $t=1,\cdots,T$ by uniqueness of the weak limit in $L^{1}(\Omega)$ and the continuous embedding of $L^{2}(\Omega)\subset L^{1}(\Omega)$. Also, since $\alpha_{t}\in BV_{g_{t},0}(\Omega)$, we can extend it by [math] outside the domain $\Omega$ and we have $((\mathcal{C}u_{n})\circ\phi_{t,n}^{-1})_{e}\underset{n\rightarrow+\infty}{\rightharpoonup}(\alpha_{t})_{e}$ in $L^{2}(\mathbb{R}^{2})$, for each $t=1,\cdots,T$.

Also, from [52, Theorem 8.20], we have that $\mbox{det}\nabla\phi_{t,n}\underset{n\rightarrow+\infty}{\rightharpoonup}\mbox{det}\nabla\bar{\phi}_{t}$ in $L^{2}(\Omega)$ for each $t=1,\cdots,T$, and by uniqueness of the weak limit in $L^{2}(\Omega)$ and the continuous embedding $L^{4}(\Omega)\subset L^{2}(\Omega)$, we deduce that $\delta_{t}=\mbox{det}\nabla\bar{\phi}_{t}$ for each $t=1,\cdots,T$. Also, since we have for $\tilde{q}=2+\frac{1}{2}>2$ :

[TABLE]

using Hölder’s inequality with $\tilde{p}=\frac{4}{\tilde{q}}$, $r=\frac{4}{\tilde{q}-1}$ and noticing that $\frac{4(\tilde{q}-1)}{4-\tilde{q}}=4$. This quantity is uniformly bounded according to $n$ from what precedes and we deduce from [53, Theorem 1 and 2] that $(\phi_{t,n})$ are bi-Hölder’s homeomorphisms and therefore $\phi_{t,n}^{-1}$ exists.

Weak lower semi-continuity: $W_{Op}$ is convex and continuous. If $\psi_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\bar{\psi}$ in $W^{1,4}(\Omega,\mathbb{R}^{2})$, then $\nabla\psi_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\nabla\bar{\psi}$ in $L^{4}(\Omega,M_{2}(\mathbb{R}))$ and we can extract a subsequence still denoted $(\nabla\psi_{n})$ such that $\nabla\psi_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\nabla\bar{\psi}$ almost everywhere in $\Omega$. If $\delta_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\bar{\delta}$ in $L^{4}(\Omega)$, then we can extract a subsequence still denoted $(\delta_{n})$ such that $\delta_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\bar{\delta}$ almost everywhere in $\Omega$. Then by continuity of $W_{Op}$, we get that $W_{Op}(\nabla\psi_{n}(x),\delta_{n}(x))\underset{n\rightarrow+\infty}{\longrightarrow}W_{Op}(\nabla\bar{\psi}(x),\bar{\delta}(x))$ almost everywhere in $\Omega$. Then by applying Fatou’s lemma, we have that $\underset{n\rightarrow+\infty}{\lim\inf}\int_{\Omega}W_{Op}(\nabla\psi_{n}(x),\delta_{n}(x))\,dx\geq\int_{\Omega}W_{Op}(\nabla\bar{\psi}(x),\bar{\delta}(x))\,dx$. Since $W_{Op}$ is convex, so is $\int_{\Omega}W_{Op}(\xi(x),\delta(x))\,dx$, and we can apply [54, Corollaire III.8] to get that $\int_{\Omega}W_{Op}(\xi(x),\delta(x))\,dx$ is lower semicontinuous in $L^{4}(\Omega,M_{2}(\mathbb{R}))\times L^{4}(\Omega)$. We deduce that $\underset{n\rightarrow+\infty}{\lim\inf}\,\int_{\Omega}W_{Op}(\nabla\phi_{t,n}(x),\mbox{det}\nabla\phi_{t,n}(x))\,dx\geq\int_{\Omega}W_{Op}(\nabla\bar{\phi}_{t}(x),\mbox{det}\nabla\bar{\phi}_{t}(x))\,dx$.

Also since $(\mathcal{C}u_{n})\circ\phi_{t,n}^{-1}\underset{n\rightarrow+\infty}{\longrightarrow}\alpha_{t}$ in $L^{1}(\Omega)$ and so in $L^{1}(\Omega,g_{t})$, then by the semi-continuity theorem from [51, Theorem 3.2], we conclude that $\operatorname{TV}_{g_{t}}(\alpha_{t})\leq\underset{n\rightarrow+\infty}{\lim\inf}\,\operatorname{TV}_{g_{t}}((\mathcal{C}u_{n})\circ\phi_{t,n}^{-1})$ for all $t=1,\cdots,T$.

$\mathcal{F}$ is a linear operator and continuous for the strong topology from $L^{2}(\mathbb{R}^{2})$ to $L^{2}(\mathbb{R}^{2})$. Therefore, by applying [54, Theorem III.9], $\mathcal{F}$ is continuous from the weak topology of $L^{2}(\mathbb{R}^{2})$ to the weak topology of $L^{2}(\mathbb{R}^{2})$. As $((\mathcal{C}u_{n})\circ\phi_{t,n})_{e}\underset{n\rightarrow+\infty}{\rightharpoonup}(\alpha_{t})_{e}$ in $L^{2}(\mathbb{R}^{2})$, we deduce that $\mathcal{F}((\mathcal{C}u_{n})\circ\phi_{t,n})_{e}\underset{n\rightarrow+\infty}{\rightharpoonup}\mathcal{F}(\alpha_{t})_{e}$ and thus $\mathcal{F}((\mathcal{C}u_{n})\circ\phi_{t,n})_{e}-x_{t}\underset{n\rightarrow+\infty}{\rightharpoonup}\mathcal{F}(\alpha_{t})_{e}-x_{t}$ in $L^{2}(\mathbb{R}^{2})$. By the lower semi-continuity of the norm, we deduce that $\underset{n\rightarrow+\infty}{\lim\inf}\|\mathcal{F}((\mathcal{C}u_{n})\circ\phi_{t,n})_{e}-x_{t}\|_{L^{2}(\mathbb{R}^{2})}^{2}\geq\|\mathcal{F}(\alpha_{t})_{e}-x_{t}\|_{L^{2}(\mathbb{R}^{2})}^{2}$.

We now need to prove that $(\mathcal{C}u_{n})\circ\phi_{t,n}^{-1}\circ\phi_{t,n}=\mathcal{C}u_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\alpha_{t}\circ\bar{\phi}_{t}=\bar{U}$ in $L^{p}(\Omega)$ for all $t=1,\cdots,T$. We first have :

[TABLE]

We now focus on the first term and make the change of variable $y=\phi_{t,n}(x)\Leftrightarrow x=\phi_{t,n}^{-1}(y)$ and $dy=\mbox{det}\nabla\phi_{t,n}(x)\,dx,\,dx=\frac{1}{\mbox{det}\nabla\phi_{t,n}(\phi_{t,n}^{-1}(y))}\,dy$ :

[TABLE]

using Hölder’s inequality and making another change of variable. Since $(\mathcal{C}u_{n})\circ\phi_{t,n}^{-1}\underset{n\rightarrow+\infty}{\longrightarrow}\alpha_{t}$ in $L^{\frac{5p}{4}}(\Omega)$, as $\frac{5p}{4}<2$, and $\|\frac{1}{\mbox{det}\nabla\phi_{t,n}}\|_{L^{4}(\Omega)}$ is uniformly bounded, we deduce that $\int_{\Omega}|(\mathcal{C}u_{n})\circ\phi_{t,n}^{-1}\circ\phi_{t,n}-\alpha_{t}\circ\phi_{t,n}|\,dx\underset{n\rightarrow+\infty}{\longrightarrow}0$.

According to [55, Theorem 6.70], there exists a sequence $(\xi_{t}^{k})_{k}$ of functions in $\mathcal{C}^{\infty}_{c}(\Omega)$ such that $\|\alpha_{t}-\xi_{t}^{k}\|_{L^{1}(\Omega)}\underset{k\rightarrow+\infty}{\longrightarrow}0$ and $\int_{\Omega}|\nabla\xi_{t}^{k}|\underset{k\rightarrow+\infty}{\longrightarrow}\int_{\Omega}|\nabla\alpha_{t}|+\int_{\partial\Omega}|\alpha_{t}|\,dx$ with $\int_{\partial\Omega}|\alpha_{t}|\,dx=0$ since $\alpha_{t}=0$ on $\partial\Omega$, for all $t=1,\cdots,T$. Thus $(\xi_{t}^{k})$ is uniformly bounded according to $k$ in $BV(\Omega)\underset{c}{\subset}L^{q}(\Omega)$, for $q\in[1,2[$, and therefore $\xi_{t}^{k}\underset{k\rightarrow+\infty}{\longrightarrow}\alpha_{t}$ in $L^{q}(\Omega)$ for $q\in[1,2[$. Let $\varepsilon>0$. Thus we fix $N\in\mathbb{N}^{*}$ such that $\|\xi_{t}^{N}-\alpha_{t}\|_{L^{\frac{5p}{4}}(\Omega)}\leq\frac{\varepsilon}{3}$. We now have :

[TABLE]

with $L_{\varepsilon}$ the Lipschitz constant of $\xi_{t}^{N}$ since it belongs to $\mathcal{C}^{\infty}_{c}(\Omega)$ and so is Lipschitz continuous, and using a change of variable. As $\phi_{t,n}\underset{n\rightarrow+\infty}{\rightharpoonup}\bar{\phi}_{t}$ in $W^{1,4}(\Omega)\underset{c}{\subset}L^{p}(\Omega,\mathbb{R}^{2})$, there exists $K\in\mathbb{N}^{*}$ such that for any $n\geq K$, $\|\phi_{t,n}-\bar{\phi}_{t}\|_{L^{p}(\Omega,\mathbb{R}^{2})}\leq\frac{\varepsilon}{3L_{\varepsilon}}$. From now on, we assume $n$ satisfies $n\geq K$, and we use Hölder’s inequality :

[TABLE]

with $\|\frac{1}{\mbox{det}\nabla\phi_{t,n}}\|_{L^{4}(\Omega)}$ uniformly bounded according to $n$ and $\|\frac{1}{\mbox{det}\nabla\bar{\phi}_{t}}\|_{L^{4}(\Omega)}$ bounded from what precedes. So by letting $\varepsilon$ tend to [math], we obtain that $\int_{\Omega}|\alpha_{t}\circ\phi_{t,n}-\alpha_{t}\circ\bar{\phi}_{t}|^{p}\,dx\underset{n\rightarrow+\infty}{\longrightarrow}0$ and consequently $(\mathcal{C}u_{n})\circ\phi_{t,n}^{-1}\circ\phi_{t,n}=\mathcal{C}u_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\alpha_{t}\circ\bar{\phi}_{t}$ for all $t=1,\cdots,N$. By uniqueness of the limit, we have that $\bar{U}=\alpha_{t}\circ\bar{\phi}_{t}\Leftrightarrow\alpha_{t}=\bar{U}\circ\bar{\phi}_{t}^{-1}$ for all $t=1,\cdots,T$ in $L^{p}(\Omega)$, with $\bar{U}\in L^{p}(\Omega)$ and $\bar{U}\circ\bar{\phi}_{t}^{-1}\in BV_{g_{t},0}(\Omega)$.

We now set $\tilde{u}_{n}=\frac{1}{|\Omega^{\prime}|}\int_{\Omega^{\prime}}u_{n}\,dx$, and $u_{0,n}=u_{n}-\tilde{u}_{n}$. We clearly have $\int_{\Omega^{\prime}}u_{0,n}\,dx=0$ for all $n$, and $\operatorname{TV}(u_{0,n})=\operatorname{TV}(u_{n})$ is uniformly bounded thanks to the coercivity inequality. We denote this uniform bound by $\nu$. Thus using Poincaré-Wirtinger’s inequality we obtain

[TABLE]

with $c_{1}>0$. We now need a bound for $\|\mathcal{C}\tilde{u}_{n}\|_{L^{p}(\Omega)}$:

[TABLE]

as $\mathcal{C}u_{n}$ strongly converges in $L^{p}(\Omega)$. Besides, we have $\|\mathcal{C}\|_{\mathcal{L}(L^{1}(\Omega^{\prime}),L^{p}(\Omega))}\|u_{0,n}\|_{L^{1}(\Omega^{\prime})}\leq\|\mathcal{C}\|_{\mathcal{L}(L^{1}(\Omega^{\prime}),L^{p}(\Omega))}c_{1}\nu=c_{3}<+\infty$ and thus $\|\mathcal{C}\tilde{u}_{n}\|_{L^{p}(\Omega)}\leq 2\sqrt{c_{3}^{2}+c_{2}^{2}}$. But we also know that $\|\mathcal{C}\tilde{u}_{n}\|_{L^{p}(\Omega)}=\frac{1}{|\Omega^{\prime}|}|\int_{\Omega^{\prime}}u_{n}\,dx|\|\mathcal{C}\mathbb{1}\|_{L^{p}(\Omega)}\leq 2\sqrt{c_{2}^{2}+c_{3}^{2}}$. Since $\mathcal{C}\mathbb{1}\neq 0$, we have $\frac{1}{|\Omega^{\prime}|}|\int_{\Omega^{\prime}}u_{n}\,dx|\leq\frac{2\sqrt{c_{2}^{2}+c_{3}^{2}}}{\|\mathcal{C}\mathbb{1}\|_{L^{p}(\Omega)}}=c_{4}<+\infty$. We therefore have

[TABLE]

Thus $u_{n}$ is uniformly bounded according to $n$ in $BV(\Omega^{\prime})$ and there exists a subsequence still denoted $(u_{n})$ such that $u_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\bar{u}$ in $L^{1}(\Omega^{\prime})$ with $\bar{u}\in BV(\Omega^{\prime})$. By continuity of the operator $\mathcal{C}$ and the uniqueness of the limit, we deduce that $\mathcal{C}u_{n}\underset{n\rightarrow+\infty}{\longrightarrow}\mathcal{C}\bar{u}=\bar{U}$ in $L^{p}(\Omega)$. By the semi-continuity theorem, we get $\operatorname{TV}(\bar{u})\leq\underset{n\rightarrow+\infty}{\lim\inf}\,\operatorname{TV}(u_{n})$.

By combining all the results, we obtain that $+\infty>\underset{(u,(\phi_{t})_{t=1,\cdots,T})\in\mathcal{U}}{\inf}G(u,(\phi_{t})_{t=1,\cdots,T})=\underset{n\rightarrow+\infty}{\lim\inf}\,G(u_{n},(\phi_{t,n})_{t=1,\cdots,T})\geq\frac{1}{T}\underset{t=1}{\overset{T}{\sum}}\delta\operatorname{TV}_{g_{t}}((\mathcal{C}\bar{u})\circ\bar{\phi}_{t}^{-1})+\frac{1}{2}\|\mathcal{F}((\mathcal{C}\bar{u})\circ\bar{\phi}_{t}^{-1})_{e}-x_{t}\|_{L^{2}(\mathbb{R}^{2})}^{2}+\int_{\Omega}W_{Op}(\nabla\bar{\phi}_{t})\,dx+\alpha\operatorname{TV}(\bar{u})$. It thus means that $\mbox{det}\nabla\bar{\phi}_{t}\in L^{4}(\Omega)$, $\frac{1}{\mbox{det}\nabla\bar{\phi}_{t}}\in L^{4}(\Omega)$, and $\mbox{det}\nabla\bar{\phi}_{t}>0$ almost everywhere in $\Omega$ for all $t=1,\cdots,T$. Indeed, since $W_{Op}(\nabla\bar{\phi}_{t}(x),\mbox{det}\nabla\bar{\phi}_{t}(x))=+\infty$ when $\mbox{det}\nabla\bar{\phi}_{t}(x)\leq 0$, it means that the set on which it happens must be of null measure otherwise we would have $\int_{\Omega}W_{Op}(\nabla\bar{\phi}_{t},\mbox{det}\nabla\bar{\phi}_{t})\,dx=+\infty$. Also, by applying the same reasoning for each $\phi_{t,n}$, we prove that $\bar{\phi}_{t}$ is a bi-Hölder homeomorphism and have that $\bar{\phi}_{t}\in\mathcal{W}$ for each $t=1,\cdots,T$.

We thus have proved the existence of minimisers for our problem (4).

∎