# Learning optical flow for fast MRI reconstruction

Timothée Schmoderer, Angelica I Aviles-Rivero, Veronica Corona, Noémie

Debroux, Carola-Bibiane Schönlieb

► To cite this version:

Timothée Schmoderer, Angelica I Aviles-Rivero, Veronica Corona, Noémie Debroux, Carola-Bibiane Schönlieb. Learning optical flow for fast MRI reconstruction. Inverse Problems, 2021, 37 (9), pp.095007. 10.1088/1361-6420/ac164a. hal-03329929

HAL Id: hal-03329929

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

Submitted on 31 Aug 2021

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

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

<figure><img src="image_1.png" /></figure>
Learning Optical Flow for Fast MRI Reconstruction

Timothée Schmoderer¹, Angelica I Aviles-Rivero²,
Veronica Corona³, Noémie Debroux⁴ and
Carola-Bibiane Schönlieb³

¹ Normandie Univ, Institut National des Sciences Appliquées de Rouen, Laboratory of Mathematics, 76000 Rouen, FR

E-mail: timothee.schmoderer@insa-rouen.fr

² Department of Pure Mathematics and Mathematical Statistics, University of Cambridge, UK

E-mail: ai323@cam.ac.uk

³ Department of Applied Mathematics and Theoretical Physics (DAMTP), Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 OWA, UK

E-mail: vc324@cam.ac.uk, cbs31@cam.ac.uk

⁴ Université Clermont Auvergne, CNRS, SIGMA Clermont, Institut Pascal, F-63000 Clermont-Ferrand, FR

E-mail: noemie.debroux@uca.fr

June 2021

**Abstract.** Reconstructing high-quality magnetic resonance images (MRI) from undersampled raw data is of great interest from both technical and clinical point of views. To this date, however, it is still a mathematically and computationally challenging problem due to its severe ill-posedness, resulting from the highly undersampled data. Whilst a number of techniques have been presented to improve image reconstruction, they only account for spatio-temporal regularisation, which shows its limitations in several relevant scenarios including dynamic data. In this work, we propose a new mathematical model for the reconstruction of high-quality medical MRI from few measurements. Our proposed approach combines - *in a multi-task and hybrid model* - the traditional compressed sensing formulation for the reconstruction of dynamic MRI with motion compensation by learning an optical flow approximation. More precisely, we propose to encode the dynamics in the form of an optical flow model that is sparsely represented over a learned dictionary. This has the advantage that ground truth data is not required in the training of the optical flow term. Furthermore, we present an efficient optimisation scheme to tackle the non-convex problem based on an alternating splitting method. We demonstrate the potentials of our approach through an extensive set of numerical results using different datasets and acceleration factors. Our combined approach reaches and outperforms several state-of-the-art techniques for multi-tasking reconstruction and other classic variational reconstruction schemes. Finally, we show the ability of our technique to transfer phantom based knowledge to real datasets.

Keywords: MRI Reconstruction, Optical Flow, Dictionary Learning, Multi-Task Model
Submitted to: *Inverse Problems*

# 1. Introduction

The problem of reconstructing high quality images, from dynamic Magnetic Resonance (MR) measurements, whilst reducing the inherent involuntary motion during the acquisition is a central topic in the community. Although it has been widely explored, it is still considered an open problem that we address in this work. This central topic in MRI comes as a consequence of the long acquisition time, which makes the image formation highly sensitive to motion leading to image degradation and compromising the expert interpretation [67, 88].

There are several attempts in the body of literature, to improve the reconstruction under inherent motion. A set of solutions to reduce motion artefacts includes breath-holding techniques [59, 88, 23] and gating strategies [70, 60, 35, 25, 33, 88]. The former rely on the patients' ability to hold their breath. However, the time needed to form an image is much longer than the average time a human being is capable of holding its breath [59]. The latter aim at tracking either the breathing or the cardiac cycles via external sensors. However, the co-registration of these signals to the image ones is still challenging [29, 88] and these techniques are effective only for perpetual motions disregarding any other involuntary motion or arrhythmia [88] and thus are only partially accurate.

The idea of accelerating the MRI acquisition by using under-sampling strategies comes as a complement or alternative to the aforementioned techniques. This option has gained great interest since the seminal paper [43] in which Compressed Sensing techniques [15, 11] were successfully adapted to the MRI reconstruction problem and adopted in clinical practice. CS relies on the property of MR images to be sparse in a transformed domain. Examples are [45, 55, 44, 40, 37, 21, 34, 39], in which authors applied several known or learnt transformations alone or in combination with parallel imaging.

Extending the idea of CS, another set of works has been devoted to low-rank matrix completion by taking advantage of the highly redundant information present in dynamic MRI sequences. For example, the work of Liang [37], in which a singular value decomposition (SVD) is used to compute the temporal basis functions, has paved the way to techniques relying on Principal Component Analysis (PCA) to reconstruct from under-sampled measurements [57, 18, 80, 1]. Other works have proposed to incorporate local low rank constraints in small regions to reduce the computational load [77, 90, 49], or a combination of sparse and low-rank constraints e.g. [39, 48, 91]. However they come at the expense of block artefacts [68].

Previous techniques were devoted to address a single task to improve MRI reconstruction. Another set of techniques were devoted to design models that seek to take advantages of the shared information among several tasks. Following this perspective, different works have been reported using so-called joint models, in which
the main idea is to intertwine several related tasks into a single unified functional. The potential in terms of image quality improvement whilst reducing error propagation has been demonstrated in [3, 42, 12, 92]. There are several works that follow this idea either for reconstructing a single image e.g.[52, 14, 83, 61, 28] or a full image sequence e.g. [38, 50, 51, 9, 89, 92]. Some of them rely on prior and external knowledge to generate motion models from sensors, navigators, etc. [50] or driven by previously computed eigenmodes of cardiac and respiratory motions [51]. Iterative schemes, with two coupled inverse problems but not derived from the same original optimisation problem, also aim at computing simultaneously the motion map and the reconstructed images e.g. [65, 66, 19, 2, 63]. These models alternate between the resolution of an independent motion estimation problem, and the resolution of a motion informed reconstruction of MR sequences.

Motivated by the success of joint models, in this work, our proposed approach follows the same perspective as the works of [3, 92, 9], in which Compressed Sensing and Optical Flow (OF) principles are combined. We emphasise that the main difference stands in the optical flow formulation. That is – we propose a new optimisation model, in which the flow estimation is encoded in a sparse representation over a learnt dictionary. More precisely, we propose a new multi-task framework for reconstructing high quality MR images from reduced measurements. Our approach seeks to improve the reconstruction quality through a better optical flow estimation using a learning-based approach. More precisely, we plug-in to the reconstruction energy a learnt sparse optical flow loss. Our optical flow model is computed as a sparse linear combination of basis functions from a learnt dictionary. Whilst these are important part of our work, our contributions are as follows:

* We introduce a unified joint and hybrid (i.e. a model that combines variational and learning strategies) model for fast MRI reconstruction, which we called **MRI Reconstruction with Dictionary Learnt Motion Compensation (MRIR-DLMC)**, in which we highlight:

- We propose a mathematical model that simultaneously performs the MRI reconstruction using Compressed Sensing principles and motion estimation via a TV-$L^1$-based dictionary learning approach for Optical Flow estimation.

- We propose a tractable and efficient numerical scheme to solve the initial optimisation problem.

- We show that our proposed technique is able to do transfer knowledge from phantom to realistic data.

* We extensively evaluate our approach using several datasets and acceleration factors. From the results,

- We show that integrating motion information in the Compressed Sensing MRI reconstruction model improves the overall quality of the reconstructed sequence in comparison to mathematical variational state of the art MRI reconstruction techniques.
- We demonstrate that the quality of the motion model is linked to the quality of the reconstruction, and that learned parametrisation for the motion as the one proposed in this paper outperforms simple motion models with handcrafted regularisation on the velocity field.

- We demonstrate generalisation properties by showing improved performance over four datasets and using several undersampling factors.

## 2. Related Work

The body of literature has reported several techniques for Fast MRI reconstruction using either single or joint models. In particular, the focus of our technique is joint models along with learning optical flow techniques. In this section, we review the literature in turn.

**Joint Models for Fast MRI reconstruction.** The literature has reported promising techniques for fast MRI reconstruction, which mainly focuses on exploiting solely compressed sensing principles e.g. [45, 55, 44, 40, 37, 21, 34, 39]. Another set of researchers have explored different alternatives the so-called joint models. The core idea of this perspective is to share knowledge from different tasks to improve performance. A challenging part when designing joint models is to select which tasks and how they will be sharing information. In this work, we are particularly interested in the setting where one has dynamic MRI sequences. To this purpose, we include motion estimation as proxy for improving MRI reconstruction. Following this perspective, different works have been reported using so-called joint models, in which the main idea is to intertwine several related tasks into a single unified functional. The potential in terms of image quality improvement whilst reducing error propagation has been shown in [3, 42, 12, 92]. There are several works that follow this idea either for reconstructing a single image e.g. [52, 14, 83, 61, 28] or a full image sequence e.g. [6, 79, 10, 51].

Authors of that [38] proposed the self-contained Deformation Corrected-Compressed Sensing (DC-CS) model relying on CS principles along with a Demons algorithm to estimate the motion. More recently, Zhang et al. in [89] extended this model to the concept of Blind Compressed Sensing, in which the transformation leading to sparsity is also estimated. In particular, our proposed approach follows the same perspective as the works of [3, 92, 9], in which CS and OF principles are combined. We emphasise that the main difference stands in the optical flow formulation. That is- we propose a new optimisation model, in which the flow estimation is encoded in a sparse representation over a learnt dictionary.

**Motion Estimation.** There are two popular strategies for motion estimation based either on registration [72] or optical flow [20]. In this work, we focus on the latter. Optical flow relies on the brightness constancy assumption, and consists in estimating the velocities of movement of brightness patterns in an image. Optical flow is a very ill-posed problem due to the aperture problem. In the body of literature, there are two main categories of techniques when designing optical flow: sparse and dense optical flow
techniques e.g. [86, 26, 73, 87, 41]. Whilst sparse optical flow seeks to estimate the
relative motion of some selected features of interest, dense optical flow outputs perpixel
estimation of the motion in the scene. In this work, we focus on dense estimation - that
is, we seek to compute a pixelwise flow estimation.

Variational approaches were first introduced to solve the OF problem, starting with
the seminal work of Horn and Schunck [30]. In that work, authors introduced an $L^2$
fidelity term coming from the linearised brightness constancy assumption and an $L^2$
regularity term on the flow fields to ensure smoothness and make the problem well-
posed. Since then, there has been a large body of works trying to improve the quality
and the accuracy of the obtained flows in a variational setting [20, 58].

More recently, several authors have explored the optical flow problem using tools
drawn from machine learning. Astonishing results have been reported using deep
learning either for normal scene conditions e.g.[75, 16, 31] or extreme scenarios e.g.[36].
We remark that our purpose in this work is different from those works. Whilst they aim
to improve the optical flow estimation itself, *our goal is instead to use OF as proxy to
improve the MRI reconstruction*. In particular, we are interested in learning the optical
flow using sparse coding. Our motivation to follow this philosophy, is given by the
promising results reported in image restoration e.g [17, 46]. The central idea of sparse
coding is to find a dictionary such that each measurement can be well-approximated by
a sparse linear combination of atoms (i.e. elements of a dictionary).

There have been several attempts to learn a sparse optical flow. The works of
that [76, 71] used known basis functions whilst authors in [32] learned an overcomplete
dictionary coming either from ground truth set or from classic variational estimation
[27, 86]. The major drawback of these techniques is the lack of convergence guarantee
but this has been overcome by the work of Bao [5]. In that work, authors introduced
a globally convergent algorithm to solve generic dictionary learning and sparse coding
problems. Whilst sparse coding has been widely explored for natural images, sparse
optical flow models have yet to be tackled in the medical domain and only few works
have been reported e.g. [56].

3. MRIR-DLMC for Fast MRI: Proposed Model

This section is devoted to the depiction of our joint and hybrid model for simultaneous
MRI reconstruction and motion estimation called **MRI Reconstruction** with **Dictionary
Learnt Motion Compensation MRIR-DLMC**. The section is organised in three key
parts: i) the CS MRI reconstruction scheme ii) our proposed approach for learning
optical flow; and iii) our combined MRI reconstruction and motion estimation model.

3.1. *Compressed Sensing MRI Reconstruction*

We consider a dynamic MRI setting, where we denote by **f** the time sequence of raw
data collected in a spatial-frequency space (k, t-space). The associated forward model –
relating the MR measurements **f** with **m** the magnetisation image showing the anatomy of the patient in the image domain Ω, being a connected bounded open subset of R² – reads:

$$ \mathbf{f}(\mathbf{k}, t) = \int_{\Omega} \mathbf{m}(\mathbf{x}, t) \exp(-j\mathbf{k}^T \cdot \mathbf{x}) d\mathbf{x} + \eta(\mathbf{k}, t), \quad (1) $$

where **x** is the spatial coordinate and **k** the frequency variable. Moreover, *t* is the temporal coordinate varying from 0 (initial acquisition time) to *T* (final acquisition time), *η* is the inherent acquisition noise which can be modelled as Gaussian noise for the MRI application and **m** is the image sequence of a moving part of the body.

For fast MRI reconstruction, one seeks to reconstruct the images from a highly reduced number of acquired measurements, which translates into $\mathbf{f}(\mathbf{k}, t)$ being available only for a small amount of $\mathbf{k} \in \mathbb{R}^2$. To do this, we denote by $K$ the undersampled Fourier operator, that is to say the composition of a sampling mask also encoding coil sensitivities (i.e. sensitivity profile of a phased array coil element) with the Fourier transform. Nevertheless, the resulting inverse reconstruction problem becomes highly ill-posed and can only be solved by introducing a prior knowledge on the nature of the reconstructed images.

The MRI reconstruction from highly undersampled data, using compressed sensing principles, was first investigated in the seminal paper of Lustig [43]. The main idea of CS in MRI is to take advantage of the natural sparsity of images in some transform domains, and especially the sparsity of MRI in the Wavelet domain. One can use the $L^1$ norm as a sparsity measure for which the unconstrained CS reconstruction minimisation problem reads:

$$ \inf_{\mathbf{m}} E(\mathbf{m}) = \int_{0}^{T} \frac{1}{2} \|K\mathbf{m} - \mathbf{f}\|_{L^2(\mathbb{R}^2)}^2 + \lambda_1 TV(\mathbf{m}) + \lambda_2 \|\Psi\mathbf{m}\|_{L^1(\mathbb{R}^4)} dt, \quad (2) $$

where $TV(\mathbf{m})$ denotes the spatial total variation to enforce sparse edges and $\Psi$ indicates the Wavelet transform. Moreover, $\lambda_1$ and $\lambda_2$ are non-negative tuning parameters. The first term is a matching criterion to enforce closeness between the reconstructed images **m** and the MR measurements **f**, while the remaining terms regularise the images by adding prior knowledge and making the problem well-posed. In this work, we focus on the total variation combined with the Wavelet transform favouring better medical image reconstruction with sharp edges, but this can be replaced by any other transformation easily such as the Fourier transform or the discrete cosine transform.

Although the body of literature has shown the potentials of using the model described in (2), it is not time dependent and therefore is not performing well enough in dynamic settings, where motion and temporal redundancy are involved. With the purpose to improve the reconstruction quality, one can account for the inherent motion in the scene via motion estimation that can be computed using Optical Flow (OF). In what follows, we describe our new approach to improve OF approximation via a sparse representation in a learnt dictionary, and show its connection with MRI reconstruction.
## 3.2. Motion Estimation: Learning Optical Flow

Optical flow (OF) consists in estimating the vector field of apparent velocities of brightness patterns between consecutive frames. Assuming small displacements, the essential brightness constancy assumption can be linearised and reads:

$$ \frac{\partial \mathbf{m}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{m} = 0, \qquad (3) $$

where $\mathbf{m} = \mathbf{m}(\mathbf{x}, t)$ denotes the dynamic sequence of images with $\frac{\partial \mathbf{m}}{\partial t}$ being the temporal derivative of the image sequence, $\nabla \mathbf{m}$ the spatial gradient, and $\mathbf{u} = [u_x, u_y]^T$ the unknown motion field. To deal with the aperture problem, one can embed (3) in a variational formulation which reads, in its generic form:

$$ \inf_{\mathbf{u}} E(\mathbf{u}) = \int_{0}^{T} E_{d}(\mathbf{u}) + \lambda E_{r}(\mathbf{u}) dt, \qquad (4) $$

where $E_d$ is a fidelity term of the form $E_d(\mathbf{u}) = \phi(\frac{\partial \mathbf{m}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{m})$, and $E_r$ is a regularisation of the flow field to make the problem well-posed. In their seminal paper [30], Horn and Schunck set $\phi(\cdot) = \| \cdot \|_{L^2(\Omega)}^2$ and $E_r(\cdot) = \| \nabla \cdot \|_{L^2(\Omega)}^2$. Moreover, in the work of [58] authors proposed $\phi(\cdot) = \| \cdot \|_{L^1(\Omega)}$ and $E_r(\cdot) = TV(\cdot)$, which is known as the $TV - L^1$ formulation. They show that this formulation leads to better and more accurate results than the ones in [30]. In the remainder of the paper, we build upon the $TV - L^1$ model.

**Remark.** The brightness constancy assumption usually is not perfectly fulfilled for real-world problems, e.g. [74, 84], and in particular, with medical data due to the rapid changes in the scene (organs). However, we show that this assumption still enables us to get a good approximation of the observed motion, for our application at hand, to improve the MRI reconstructions. A relaxation of this assumption [24, 7], for further correcting the motion estimation, to obtain a better reconstruction quality will be the object of future work. Some examples include extending the brightness constancy to high-order constancy [7], integrating local neighbourhood information [8] and using high-order random field [64].

Whilst the variational formulations for the OF estimation is limited by its own design, one can further improve the approximation of the motion scene by expressing the flow field as a sparse linear combination of basis functions in either off-the-shelf dictionaries [76, 71] or learnt ones [32, 86, 27]. Inspired by these works, we seek to increase the accuracy of the OF estimation by encoding it in a sparse representation over a learnt over-complete dictionary. From now on, to comply with the dictionary learning and sparse representation philosophies, we consider a discrete spatial setting. The overcomplete flow dictionary $D$ is thus learnt from a partition $\mathcal{P}$ of an initial flow estimation decomposed into small overlapping patches. In here, we define the patch extracting operator as $R_p$ at pixel $p$. We follow the strategy from [32], in which the flow dictionary $D$ is simplified by decoupling the horizontal and vertical motions as
$D = \begin{bmatrix} D_x & 0 \\ 0 & D_y \end{bmatrix}$, where both $D_x$ and $D_y$ are composed of $N_d$ elements such that $N_d \times |\mathcal{P}|$ is larger than the image dimension. Formally, this sparse representation over a learnt dictionary is incorporated in the variational setting through this additional term:

$$E_{\text{sparse}}(\mathbf{u}, D, \mathbf{a}) = \sum_{\mathbf{p} \in \mathcal{P}} \|R_{\mathbf{p}}\mathbf{u} - D\mathbf{a}_{\mathbf{p}}\|_{F}^{2} + \tau \|\mathbf{a}_{\mathbf{p}}\|_{1} \quad (5) \\ \text{s.t. } \|D_{x,i}\|_{2} \le 1, \|D_{y,i}\|_{2} \le 1, 1 \le i \le N_d,$$

where $\| \cdot \|_F$ is the Frobenius norm and $\| \cdot \|_1$ is the $l^1$ norm. The new variable $\mathbf{a}_p = ((a_x, a_y)_p)^T$ represents the sparse coefficients of the patch $R_p u$ in $D$. In the model (5), the first term ensures the flow field is given by a linear combination of basis functions in the dictionary whilst the second one enforces sparsity of the coefficients. At the learning level, the model (5) can be solved using either ground truth when available [32], or from an initial approximation of the flows using a variational model [86, 27]. The assumption of having ground truth is not realistic in the medical domain, and therefore we follow the second strategy. Motivated by the results in [86, 27] showing better results with a dictionary learnt from a classic variational formulation than the initial ones, we propose to build our dictionary upon the $TV-L^1$ model, and thus consider the following discrete optical flow formulation:

$$E_{\text{OF}}(\mathbf{u}, D, \mathbf{a}) = \int_{0}^{T} \lambda_{3} \left\| \frac{\partial \mathbf{m}}{\partial t} + \nabla \mathbf{m} \cdot \mathbf{u} \right\|_{1} + \lambda_{4} TV(\mathbf{u}) + \sum_{\mathbf{p} \in \mathcal{P}} \lambda_{5} \|R_{\mathbf{p}}\mathbf{u} - D\mathbf{a}_{\mathbf{p}}\|_{F}^{2} + \lambda_{6} \|\mathbf{a}_{\mathbf{p}}\|_{1} \quad (6) \\ \text{s.t. } \|D_{x,i}\|_{2} \le 1, \|D_{y,i}\|_{2} \le 1, 1 \le i \le n_d,$$

where $\lambda_3, \lambda_4, \lambda_5$ and $\lambda_6 > 0$ are weighting parameters to balance the influence of each term. The constraints on $D$ ensure uniqueness of the solution. Given a reference optical flow $\mathbf{u}^{\text{ref}}$, the dictionary is learnt as the result of the following minimisation problem under constraints:

$$E(D, \mathbf{a}) = \int_{0}^{T} \|R_{\mathbf{p}}\mathbf{u}_{\mathbf{p}}^{\text{ref}} - D\mathbf{a}_{\mathbf{p}}\|_{F}^{2} dt, \quad (7) \\ \|D_{x,j}\|_{2} \le 1, \|D_{y,j}\|_{2} \le 1, \quad 1 \le j \le N_{d}, \\ \|\mathbf{a}_{\mathbf{p}}\|_{0} \le k_{0} \quad \forall \mathbf{p} \in \mathcal{P},$$

where the pseudo norm $l^0$ counts the non-zero elements of $\mathbf{a}_p$. The model in (7) uses a constraint on $\mathbf{a}$ rather than the $l^1$ relaxation of the $l^0$ pseudo-norm as it shows better numerical performances.

### 3.3. Joint Model: Learning Optical Flow for MRI Reconstruction

The potentials of CS MRI reconstruction and optical flow for motion estimation have been demonstrated separately in the previous subsections. We now demonstrate
the strong correlation between both tasks – as motion estimation depends on the reconstruction quality, and the reconstruction accuracy can be substantially improved by a good approximation of the motion scene – and show their mutual benefit by intertwining them. This statement has been shown in recent works such as [3, 4, 9] that have inspired our work, in which one seeks to solve a single optimisation model unifying both tasks in a pure variational setting. More precisely, the two formulations for CS MRI reconstruction and TV-L¹ Optical Flow estimation can be cast in a unified optimisation problem as:

$$ \inf_{\mathbf{m},\mathbf{u}} \int_0^T \frac{1}{2} \|K\mathbf{m} - \mathbf{f}\|_{L^2(\mathbb{R}^2)}^2 + \lambda_1 \mathrm{TV}(\mathbf{m})^p + \lambda_2 \|\psi\mathbf{m}\|_{L^1(\mathbb{R}^4)}^p + \lambda_3 \left\| \frac{\partial \mathbf{m}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{m} \right\|_{L^1(\Omega)} \quad (8) \\ + \lambda_4 \mathrm{TV}(\mathbf{u})^q dt, $$

with $p > 1$ and $q \ge 1$. For now, we have just added the sparsity constraint on the Wavelet domain of the reconstruction sequence in a purely variational setting compared to [9, 3, 4]. The existence theorem given in [9] still holds.

**Theorem 1** (Existence of minimisers). Let the assumptions of [9, Theorem 3.1] be satisfied (the MRI reconstruction operator fulfilling the required properties). Then there exists a minimiser of the full continuous model (8) in the space

$$ \{(\mathbf{m}, \mathbf{u}) : \mathbf{m} \in L^p(0, T; \mathrm{BV}(\Omega)), \mathbf{u} \in L^q(0, T; \mathrm{BV}(\Omega))^2, \nabla \cdot \mathbf{u} \in \Theta\}, $$

$\Theta$ being a Lebesgue space including an upper bound constraint.

*Proof.* The proof is adapted from the one of [9, Theorem 3.1] and we just have to prove the weak-* lower semi-continuity property of the additional term $\|\psi\mathbf{m}\|_{L^p(0,T;L^1(\mathbb{R}^4))}$. Since the Wavelet transform is a linear operator in each direction defined by a convolution, then it is a strongly continuous operator from $L^2(\mathbb{R}^2)$ to $L^2(\mathbb{R}^4)$. Also, since the spatial dimension is two, we have the compact embedding of $BV(\Omega)$ into $L^2(\Omega)$ and using the continuous embedding of $L^2(\mathbb{R}^4)$ in $L^1(\mathbb{R}^4)$, one can prove the weak-* lower semi-continuity of the additional term $\|\psi\mathbf{m}\|_{L^p(0,T;L^1(\mathbb{R}^4))}$. $\square$

We remark that unlike the models where authors applied directly the TV-L¹ optical flow model, we seek to further improve that formulation via a learnt dictionary. To do this, we propose to encode the flow estimation in a sparse representation over a learnt dictionary. We therefore consider this new hybrid model mixing continuous and discrete variables and terms:

$$ \inf_{\mathbf{m},\mathbf{u},D,\mathbf{a}} E(\mathbf{m}, \mathbf{u}, D, \mathbf{a}) = \int_0^T \frac{1}{2} \|K\mathbf{m} - \mathbf{f}\|_F^2 + \lambda_1 \mathrm{TV}(\mathbf{m})^p + \lambda_2 \|\psi\mathbf{m}\|_1^p \quad (9) \\ + \lambda_3 \left\| \frac{\partial \mathbf{m}}{\partial t} + \nabla \mathbf{m} \cdot \mathbf{u} \right\|_1 + \lambda_4 \mathrm{TV}(\mathbf{u})^q dt \\ + \sum_{k=1}^{N_t-1} \sum_{\mathbf{p} \in \mathcal{P}} \lambda_5 \|R_{\mathbf{p},k}\mathbf{u} - D\mathbf{a}_{\mathbf{p},k}\|_F^2 + \lambda_6 \|\mathbf{a}_{\mathbf{p},k}\|_1, $$
where $N_t$ is the discrete total number of frames. We still have the existence of minimisers theorem.

**Theorem 2 (Existence of minimisers).** Let the assumptions of [9, Theorem 3.1] be satisfied (the MRI reconstruction operator fulfilling the required properties). Then there exists a minimiser of (9) in the space

$$ \{ (\mathbf{m}, \mathbf{u}, D, \mathbf{a}) : \mathbf{m} \in L^p(0, T; BV(\Omega)), \mathbf{u} \in L^q(0, T; BV(\Omega))^2, D \in \mathbb{R}^{P_s^2 \times N_d}, \\ \mathbf{a} \in \mathbb{R}^{N_d \times P_n \times (N_t-1)}, \forall 1 \le j \le N_d, \|D_j\| \le 1, \nabla \cdot \mathbf{u} \in \Theta \}, $$

Θ being a Lebesgue space including an upper bound constraint, $R_{\mathbf{p},k}$ : $L^q(0, T; BV(\Omega))^2 \to \mathbb{R}^{P_s^2}$ being the discretisation (L²-mean over the spatio-temporal pixel domain) patch extractor at pixel $\mathbf{p} \in \mathcal{P}$ with $P_s$ the patch size and at discrete time $k$, $\mathbf{A}_{\mathbf{p},k}$ being a vector of size $N_d$ corresponding to the sparse coefficient for the patch centered at pixel $\mathbf{p} \in \mathcal{P}$ and for the time $k$.

*Proof.* The proof follows the same arguments as previously. Indeed, the coercivity as well as the weak and weak-⋆ lower semi-continuity of the continuous terms remain true since all the added terms are non-negative. We therefore focus only on the discrete terms and variables. Let $(\mathbf{m}_n, \mathbf{u}_n, D_n, \mathbf{a}_n)$ be a minimising sequence in the space described in the theorem. We thus have that $\|D_n\|$ is uniformly bounded with respect to $n$ by definition of the functional space. We also have this coercivity inequality :

$$ \begin{aligned} E(\mathbf{m}_n, \mathbf{u}_n, D_n, \mathbf{a}_n) = & \int_0^T \frac{1}{2} \|K\mathbf{m}_n - \mathbf{f}\|_F^2 + \lambda_1 TV(\mathbf{m}_n)^p + \lambda_2 \|\psi\mathbf{m}_n\|_1^p \\ & + \lambda_3 \left\| \frac{\partial \mathbf{m}_n}{\partial t} + \nabla \mathbf{m}_n \cdot \mathbf{u}_n \right\|_1 + \lambda_4 TV(\mathbf{u}_n)^q dt \\ & + \sum_{k=1}^{N_t-1} \sum_{\mathbf{p} \in \mathcal{P}} \lambda_5 \|R_{\mathbf{p},k} \mathbf{u}_n - D \mathbf{a}_{\mathbf{p},k,n}\|_F^2 + \lambda_6 \|\mathbf{a}_{\mathbf{p},k,n}\|_1, \\ & \ge \sum_{k=1}^{N_t-1} \lambda_6 \|\mathbf{a}_{\mathbf{p},k,n}\|_1. \end{aligned} $$

$E$ is proper (by taking $\mathbf{u} = 0, \mathbf{m} = 0, D = \text{Id}, \mathbf{a} = 0$), $E(\mathbf{u}, \mathbf{m}, D, \mathbf{a}) = \int_0^T \|f\|_{L^2(\Omega)}^2 dt$), positive with a finite infimum. Thus for $n$ large enough, we get $E(\mathbf{m}_n, \mathbf{u}_n, D_n, \mathbf{a}_n) \le \inf_{\mathbf{u},\mathbf{m},D,\mathbf{a}} E(\mathbf{u}, \mathbf{m}, D, \mathbf{a}) + 1 < +\infty$. Therefore $(\mathbf{a}_n)$ is uniformly bounded according to $n$. According to Bolzano-Wierstrass theorem we can thus extract subsequences from $(\mathbf{a}_n)$ and $(D_n)$ that are convergent. We name $\bar{D}$ and $\bar{a}$ the limits of these subsequences still denoted $(D_n)$ and $(\mathbf{a}_n)$ for the sake of conciseness.

Thanks to [9] and the previous proof, it only remains to prove the weak-⋆ lower semi-continuity of $\lambda_5 \|R_{\mathbf{p},k} \mathbf{u}_n - D \mathbf{a}_{\mathbf{p},k,n}\|_F^2$ with respect to $\mathbf{u}_n$ and its lower semi-continuity with respect to $D_n$ and $\mathbf{a}_n$. The continuity of the last term is obvious by continuity of the $L^1$ norm in the discrete case. By definition of the operator $R_{\mathbf{p},k}$, it is based on the $L^2$ norm which is weak-⋆ lower semi-continuous (see [9] for more details). The matrix
product and the norm are continuous which leads to the weak-⋆ lower semi-continuity of $\lambda_5 \|R_{\mathbf{p},k}\mathbf{u}_n - D\mathbf{a}_{\mathbf{p},k,n}\|_F^2$ with respect to $\mathbf{u}_n$ and its lower semi-continuity with respect to $D_n$ and $\mathbf{a}_n$. This concludes the proof. $\square$

We then follow the discretise-then-optimise strategy and derive the discrete formulation of *our proposed joint and hybrid model* after the learning step of the dictionary reads:

$$ 
\begin{aligned}
\inf_{\mathbf{m},\mathbf{u},D,\mathbf{a}} E(\mathbf{m}, \mathbf{u}, D, \mathbf{a}) = & \sum_{k=1}^{N_t-1} \frac{1}{2} \|K\mathbf{m}_k - \mathbf{f}_k\|_F^2 + \lambda_1 TV(\mathbf{m}_k) + \lambda_2 \|\psi\mathbf{m}_k\|_1 \\
& + \lambda_3 \|\frac{\partial \mathbf{m}_k}{\partial t}\|_1 + \nabla \mathbf{m}_k \cdot \mathbf{u}_k\|_1 + \lambda_4 TV(\mathbf{u}_k) \\
& + \sum_{\mathbf{p} \in \mathcal{P}} \lambda_5 \|R_{\mathbf{p},k}\mathbf{u}_k - D\mathbf{a}_{\mathbf{p},k}\|_F^2 + \lambda_6 \|\mathbf{a}_{\mathbf{p},k}\|_1 dt
\end{aligned}
\quad (10) $$

Since all the norms are equivalent in the discrete setting, we have removed the $p$ and $q$ coefficients. We now turn to the numerical resolution of this challenging problem.

## 4. MRIR-DLMC for Fast MRI: Numerical Realisation

Due to the non-convexity of the energy coming from the brightness constancy assumption, its non-differentiability with the TV and $L^1$ terms, and complex linear operators applied to the unknowns ($K$ and $\psi$), the problem eq. (10) is very challenging to solve numerically. We thus designed an alternating minimisation scheme, to solve our new hybrid model, in which we fix all the variables except one and solve the simplified sub-problem for each unknown in turn. The algorithm is therefore divided into four steps which will be described in the following sub-sections, each one relying on the primal dual Chambolle and Pock algorithm [13]. We now recall the general setting of this method. The Chambolle and Pock procedure aims at solving the nonlinear primal problem $\min_{x \in X} F(Cx) + G(x)$ with the following primal dual formulation:

$$ \min_{x \in X} \max_{y \in Y} \langle Cx, y \rangle + G(x) - F^*(y). \quad (11) $$

The hypothesis are: $X$ and $Y$ are finite dimensional real vector spaces, the map $C: X \to Y$ is a continuous linear operator, $G: X \to [0, +\infty[$ and $F^*: Y \to [0, +\infty[$ are proper, convex, lower-semicontinuous functions and $F^*$ is the Legendre-Fenchel conjugate of $F$. The following (1) converges towards a saddle-point of (11). The resolvent operators defined through

$$ x = (I + \tau \partial F)^{-1}(y) = \arg \min_x \left\{ \frac{\|x - y\|^2}{2\tau} + F(x) \right\}. \quad (13) $$

Our reconstructed sequence is initialised as the zero filling inverse Fourier reconstruction from the sub-sampled measurements. Zero filling is a common choice for initialisation as it does not add any substantial computational time. A discussion on other possible
**Algorithm 1:** Chambolle and Pock iteration

Choose $\tau, \sigma > 0$ such that $\tau\sigma\|C\|^2 < 1$, $\theta \in [0, 1]$, $(x^0, y^0) \in X \times Y$ and set $\bar{x}^0 = x^0$

Update $x^n, y^n$ and $\bar{x}^n$ as follows:

$$ 
\left\{ 
\begin{aligned}
y^{n+1} &= (I + \sigma \partial F^*)^{-1}(y^n + \sigma C \bar{x}^n) \\
x^{n+1} &= (I + \tau \partial G)^{-1}(x^n - \tau C^* y^{n+1}) \\
\bar{x}^{n+1} &= x^{n+1} + \theta(x^{n+1} - x^n)
\end{aligned} 
\right. 
\qquad (12) 
$$

initialisations can be seen in Appendix D. The first approximation of the optical flow $\mathbf{u}^0$ is computed on this image sequence by the $TV - L^1$ method. We consider a sequence of $N_t$ frames with size $N_x \times N_y$ pixels. Each frame will be decomposed in overlapping square patches of size $P_s \times P_s$ and we choose $N_d$, the number of atoms in the dictionary such that $N_d(P_s)^2 > N_x N_y$. We denote by $P_n$ the number of patches for a given frame. Then the main loop of our algorithm is described by algorithm 2.

**Algorithm 2:** Main Loop

Given a threshold $\epsilon > 0$.

Let $\mathbf{m}^0$ be a first approximation of the MRI sequence (ZF).

Let $\mathbf{u}^0$ be a first approximation of the optical flow ($TV - L^1$).

Learning the dictionary $D$ (if needed).

Let $\mathbf{a}^0$ be a first approximation of the sparse decomposition of $\mathbf{u}^0$ in $D$ (from the learning part).

repeat

Compute $\mathbf{m}^{n+1} = \arg\min_{\mathbf{m}} E(\mathbf{m}, \mathbf{u}^n, \mathbf{a}^n)$

Compute $\mathbf{u}^{n+1} = \arg\min_{\mathbf{u}} E(\mathbf{m}^{n+1}, \mathbf{u}, \mathbf{a}^n)$

Compute $\mathbf{a}^{n+1} = \arg\min_{\mathbf{a}} E(\mathbf{m}^{n+1}, \mathbf{u}^{n+1}, \mathbf{a})$

until $\|\mathbf{m}^{n+1} - \mathbf{m}^n\| < \epsilon \|\mathbf{m}^n\|$

return $\mathbf{m}^{n+1}$.

In what follows, we describe each of the four required algorithms, namely: learning the dictionary, image reconstruction, optical flow estimation, and sparse coding.

## 4.1. Learning the dictionary

In this section, we present the process of learning the dictionary for the optical flow consisting in the minimisation of eq. (7) following the strategy of [53] which we recall here for the sake of completeness. Given a discrete reference optical flow $\mathbf{u}^{ref} = (u_x^{ref}, u_y^{ref})^T \in [\mathbb{R}^{N_x \times N_y \times (N_t-1)}]^2$, dictionary learning aims at finding the best dictionary $D_{\mathbf{u}} = \begin{pmatrix} D_x & 0 \\ 0 & D_y \end{pmatrix}$ with $D_x$ and $D_y$ having $N_d$ atoms such that all the data
can be expressed as a linear combination of these atoms. Let $R_{\mathbf{p}} = \begin{pmatrix} R_{x,\mathbf{p}} & 0 \\ 0 & R_{y,\mathbf{p}} \end{pmatrix}$ be the extracting patch operator centred at pixel $\mathbf{p}$. Let $\mathcal{P}$ be the set of all patches, we set $P_n = |\mathcal{P}|$ and $R = (R_{\mathbf{p}})_{\mathbf{p}\in\mathcal{P}}$ so that the extracting patch operator $R$ goes from $[\mathbb{R}^{N_x \times N_y \times (N_t-1)}]^2$ onto $[\mathbb{R}^{P_s^2 \times P_n \times (N_t-1)}]^2$. Let $\mathbf{a} = \begin{pmatrix} a_x \\ a_y \end{pmatrix} \in [\mathbb{R}^{N_d \times P_n \times (N_t-1)}]^2$ be the sparse coefficients of the decomposition of $\begin{pmatrix} u_x \\ u_y \end{pmatrix}$ in $D_{\mathbf{u}}$.

From now on, we keep the notation $u$ for both $u_x$ and $u_y$, the notation $R$ for both $R_x$ and $R_y$, the notation $D$ for both $D_x$ and $D_y$ and the notation $a$ for both $a_x$ and $a_y$ as the learning process is decoupled for both directions and can be performed in parallel computing.

Let $\bar{u} = Ru^{ref} \in \mathbb{R}^{P_s^2 \times P_n \times (N_t-1)}$ be such that $\bar{u}_k \in \mathbb{R}^{P_s^2 \times P_n}$ is the matrix of all patch extracted from $\bar{u}$ at time $k$. And let $a = (a_k)_{k=1}^{N_t-1}$ defined by

$$ \min_{\|a\|_{l0} \le \epsilon} \sum_{k=1}^{N_t-1} \|\bar{u}_k - D a_k\|_2^2. \qquad (14) $$

The pseudo-norm $l^0$ of $a \in \mathbb{R}^{N_d \times P_n \times (N_t-1)}$ is given by $\|a\|_{l^0} = \min_{1 \le k \le N_t-1} |\{(a_k)_{i,j} \ne 0\}|$ and $\epsilon > 0$ controls the amount of sparsity. Dictionary learning performs an optimisation both on the dictionary $D$ and on the set of coefficients $a$. This joint optimisation reads

$$ \min_{D \in \mathcal{D}, a \in \mathcal{A}} E(D, a) = \sum_{k=1}^{N_t-1} \frac{1}{2} \| \bar{u}_k - D a_k \|_F^2 \qquad (15) $$

The constraint set on $D$ is $\mathcal{D} = \{D \in \mathbb{R}^{P_s^2 \times N_d}, \forall 1 \le j \le N_d, \|D_j\| \le 1\}$, the columns of the dictionary are unit normalised, the sparsity constraint set on $a$ reads $\mathcal{A} = \{a \in \mathbb{R}^{D_n \times P_n}, \forall 1 \le j \le P_n, |a_j|_0 \le \epsilon\}$.

**Remark.** These constraints ensure uniqueness of the solution, indeed if they were not here, we could always find a transformation $P$ such that $Pa$ has smaller $l^0$ norm and the transformed dictionary $DP^{-1}$ has bigger norm.

Following standard arguments, we propose a block coordinate descent method to minimise the energy:

$$ D^{n+1} \in \arg \min_{D \in \mathcal{D}} E(D, a^n), \quad a^{n+1} \in \arg \min_{a \in \mathcal{A}} E(D^{n+1}, a). \qquad (16) $$

Convergence on such alternating minimisation scheme can be shown, see for instance [78].

For the first sub-problem, we assume that $a$ is fixed, then it becomes a quadratic optimisation problem under constraints with respect to $D$:

$$ \min_{D \in \mathcal{D}} \left\| \bar{u}_k \left( \sum_{k=1}^{N_t-1} a_k^T \right) \sqrt{\sum_{k=1}^{N_t-1} a_k a_k^T} - D \sqrt{\sum_{k=1}^{N_t-1} a_k a_k^T} \right\|_F^2. \qquad (17) $$
The minimisation problem can be exactly solved as:

$$D^{n+1} = \mathrm{proj}_{\mathcal{D}} \left( D^n - \tau_1 \left( D^n \sum_{k=1}^{N_t-1} a_k a_k^T - \sum_{k=1}^{N_t-1} \bar{u}_k a_k^T \right) \right) \quad (18)$$

where the projection step is $\tau_1 < \frac{2}{\|\sum_{k=1}^{N_t-1} a_k a_k^T\|}$ and $\mathrm{proj}_{\mathcal{D}}(D) = \frac{D}{\|D\|}$.

Then we consider $D$ fixed and minimise with respect to $a$. This again leads to a constrained problem which can be solved by projections for each time step on the constraint set,

$$a_k^{n+1} = \mathrm{proj}_{\mathcal{A}} (a_k^n - \tau_2 D^T (D a_k^n - U_k)) \quad (19)$$

where this time, $\tau_2 < \frac{2}{\|DD^T\|}$. For $1 \le j \le P_n$ we denote $|\bar{a}_{k,j}(1)| \le \dots \le |\bar{a}_{k,j}(N_d)|$ the order of magnitude of the vector $a_{k,j} \in \mathbb{R}^{N_d}$. Then the projection operator on $\mathcal{A}$ reads:

$$\forall 1 \le i \le N_d, \quad \tilde{a}_{k,j}(i) = \begin{cases} a_{k,j}(i) & \text{if } |a_{k,j}(i)| \ge |\bar{a}_{k,j}(\epsilon)| \\ 0 & \text{otherwise.} \end{cases} \quad (20)$$

**Remark.** Two levels of parallel computing can be performed for this process. First we can learn the dictionary $D_x$ for the horizontal component $u$ of the optical flow and the dictionary $D_y$ for the vertical component simultaneously. Additionally, for each component we can minimise $a_k$ in parallel.

## 4.2. Image reconstruction

This section describes the process of MRI reconstruction which is the core of the method. While in [43] a smooth regularisation of the $L^1$ norm was used, in this work, we propose to solve the non-smooth optimisation problem using algorithm 1. Assume $\mathbf{u}$ is known, the problem then reads

$$\min_{\mathbf{m}} E(\mathbf{m}) = \int_{0}^{T} \frac{1}{2} \|K\mathbf{m} - \mathbf{f}\|^2 + \lambda_1 \|\nabla \mathbf{m}\|_{2,1} + \lambda_2 \|\Psi \mathbf{m}\|_{1} + \lambda_3 \left\| \frac{\partial \mathbf{m}}{\partial t} + \nabla \mathbf{m} \cdot \mathbf{u} \right\|_{1} \, dt \quad (21)$$

with $\|\nabla \cdot \mathbf{m}\|_{2,1}$ representing the discrete spatial isotropic TV.

We proceed as in [9]. Let $C = [K, \nabla, \Psi, D_t + \mathbf{u} \cdot \nabla]^T$ be the linear operator acting on $\mathbf{m}$ so that the minimisation is seen as $\min_{\mathbf{m}} E(C\mathbf{m})$. And let $\mathbf{y} = [y_0, y_1, y_2, y_3]^T$ be the collection of dual variables for each operator. From eq. (11) we have $F := E$ and $G = 0$ and the Legendre-Fenchel conjugate of the energy is given by

$$E^*(\mathbf{y}) = \int_{0}^{T} \frac{1}{2} \|y_0\|_2^2 + \langle y_0, \mathbf{f} \rangle + \delta_{\{y: \|y\|_{2,\infty} \le 1\}} \left( \frac{y_1}{\lambda_1} \right) + \delta_{\{y: \|y\|_{\infty} \le 1\}} \left( \frac{y_2}{\lambda_2} \right) \\ + \delta_{\{y: \|y\|_{\infty} \le 1\}} \left( \frac{y_3}{\lambda_3} \right) dt \quad (22)$$
with $\delta_I(y) = 0$ if $y \in I$ and $+\infty$ otherwise. The primal-dual algorithm requires the proximity operator $(I + \sigma \partial E^*)^{-1}(\mathbf{y}^n + \sigma C \bar{\mathbf{m}}^n)$. The components $y_i^{n+1}$ for $i = 0, 1, 3$ are computed exactly as in [9] and we give only

$$y_2^{n+1} = \arg \min_y \left\{ \int_0^T \frac{1}{2} \|y - y_2^n - \sigma \Psi (\mathbf{m}^n + \theta(\mathbf{m}^n - \mathbf{m}^{n-1}))\|_2^2 + \delta_{\{y: \|y\|_{\infty} \le 1\}} \left(\frac{y}{\lambda_2}\right) dt \right\}.$$

This subproblem admits a closed form solution recalled here,

$$\tilde{\mathbf{y}}^{n+1} = \mathbf{y}^n + \sigma C (\mathbf{m}^n + \theta(\mathbf{m}^n - \mathbf{m}^{n-1}))$$

$$\mathbf{y}^{n+1} = (I + \sigma \partial E^*)^{-1}(\mathbf{y}^n) = \begin{cases} y_0^{n+1} = \frac{y_0^n - \sigma \mathbf{f}}{\sigma+1}, \\ y_1^{n+1} = \pi_{\lambda_1}(\tilde{y}_1^n), \\ y_2^{n+1} = \pi_{\lambda_2}(\tilde{y}_2^n), \\ y_3^{n+1} = \pi_{\lambda_3}(\tilde{y}_3^n) \end{cases}.$$

The projection onto the unit ball is given by $\pi_\lambda(y) = \frac{y}{\max(1, \frac{\|y\|_2}{\lambda})}$. The parameters $\tau, \sigma$ refer to step sizes, and have to fulfil the following constraints: $\tau\sigma\|C\|^2 \le 1$ with $\|C\| \approx 2 + \sqrt{8} + \sqrt{2}(1 + \|\mathbf{u}\|)$.

### 4.3. Sparse flow representation and optical flow approximation

The minimisation over the variable $\mathbf{a}$ can be split for each component $\mathbf{a}_\mathbf{p}$ which is given by $\tilde{E}(\mathbf{a}_\mathbf{p}) = \lambda_5 \|R_\mathbf{p}\mathbf{u} - D\mathbf{a}_\mathbf{p}\|_F^2 + \lambda_6 \|\mathbf{a}_\mathbf{p}\|_1$. From algorithm 1 we take $F = \lambda_6 \|\cdot\|_1$, the proximal operator of its Legendre-Fenchel conjugate is given by projection operator:

$$(I + \sigma \partial F^*)^{-1}(\mathbf{y}) = \pi_{\lambda_6}(\mathbf{y}) \qquad (23)$$

Moreover since $G = \lambda_4 \|R_\mathbf{p}\mathbf{u} - D\mathbf{u}\|_2^2$ is smooth, its resolvent reduces to the gradient,

$$\tau \nabla G(\mathbf{a}_\mathbf{p}) = 2\tau\lambda_5 D^T (D\mathbf{a}_\mathbf{p} - R_\mathbf{p}\mathbf{u}). \qquad (24)$$

Finally we optimise on the optical flow $\mathbf{u}$, given fixed $\mathbf{m}$ and $\mathbf{a} = (\mathbf{a}_\mathbf{p})_{\mathbf{p} \in \mathcal{P}}$. Let $E(\mathbf{u}) = F(\nabla\mathbf{u}) + G(\mathbf{u})$ with $F(\cdot) = \lambda_4 \|\cdot\|_{2,1}$ and $G = \lambda_3 \|\frac{\partial \mathbf{m}}{\partial t} + \nabla\mathbf{m}\|_1 + \lambda_5 \sum_\mathbf{p} \|R_\mathbf{p} \cdot -D\mathbf{a}_\mathbf{p}\|_F^2$. For simplicity of the notations we introduce the operators $\mathcal{A} = I + 2\tau\lambda_5 \sum R_\mathbf{p}^T R_\mathbf{p}$, $\rho(\mathbf{u}) = \frac{\partial \mathbf{m}}{\partial t} + \mathbf{u} \cdot \nabla\mathbf{m}$ and the notation $\tilde{\mathbf{u}} = \mathbf{u} + 2\tau\lambda_5 \sum R_\mathbf{p}^T D\mathbf{a}$. Then the proximal operator of $G$ reads,

$$\mathcal{A}(I + \tau \partial G)^{-1}(\mathbf{u}) = \tilde{\mathbf{u}} + \begin{cases} -\lambda_3 \nabla \mathbf{m} & \text{if } \rho(\mathcal{A}^{-1} \tilde{\mathbf{u}}) > \tau \lambda_3 \mathcal{A}^{-1} \|\nabla \mathbf{m}\|^2, \\ \lambda_3 \nabla \mathbf{m} & \text{if } \rho(\mathcal{A}^{-1} \tilde{\mathbf{u}}) < -\tau \lambda_3 \mathcal{A}^{-1} \|\nabla \mathbf{m}\|^2, \\ -\rho(\tilde{\mathbf{u}}) \frac{\nabla \mathbf{m}}{\|\nabla \mathbf{m}\|^2} & \text{else,} \end{cases} \qquad (25)$$

and the $F^*$ one is

$$(I + \sigma \partial F^*)^{-1}(\mathbf{y}) = \pi_{\lambda_6}(\mathbf{y}). \qquad (26)$$

Details of the computation of eq. (25) can be found in appendix Appendix C and details about the discretization of the operators can be found in appendix Appendix A. We now turn to the numerical experiments.
# 5. Numerical Experiments

In this section we show the results of our numerical experiments. All tests and comparisons were run under the same CPU-based implementation (Intel Core i7-8550U) with MATLAB 2018b. Compared techniques were either used as provided by the author or implemented when not available.

## 5.1. Data Description and Evaluation Protocol

Our approach is evaluated by using data coming from: cine cardiac, right lobe liver and phantom generated cardiac cine. The datasets were saved as fully sampled raw data and then were retrospectively undersampled using a variable-density random sampling pattern, and in particular, a cartesian pattern. ‡ Figure 1 as suggested by Lustig [43]. Datasets characteristics are as follows:

* Dataset I: A cine cardiac dataset acquired from a healthy volunteer, from [69]. It is a retrospectively gated bSSFP short-axis view of a health volunteer measured on 3T scanner(Skyra, Siemens Healthcare, Erlangen, Germany) with the following scan parameters: FOV=274.62 mm×340 mm, matrix size=208×168, 25 cardiac phases with a temporal resolution of 42.72 ms, TR/TE/FA=3.56 ms/1.78 ms/40°, 6 mm slice-thickness, TA=16 s.

* Dataset II and III: Realistic cardiac cine simulation generated using the MRXCAT phantom framework [85]. Whilst the Dataset II was simulated during breath holding, the III was set with free respiratory motion. Matrix size - 128 × 128 and 24 frames.

* Dataset IV: 4DMRI dataset [81] acquired from the right liver lobe during free-breathing. It was acquired using a 1.5T Philips Achieva system with: SENSE factor 1.7 and half scan, flip angle 70°, TR=3.1ms, slices=25 and matrix size of 195 × 166, coils=4, and a temporal resolution of 2.6 – 2.8Hz.

We remark that in this work, we followed common protocol by different learning paradigms when treating multi-coils. That is, we collapsed the multi-coil data (such as Dataset IV) into a single view before performing the reconstruction.

We have followed common protocol for evaluating our MRIR-DLMC using two metrics over the full reconstructions [43, 38, 45, 47]: the peak signal-to-noise ratio (PSNR) and the structural similarity (SSIM). During the following sections we use MRIR-DLMC and OURS to refer to our technique. The metric-wise computations were performed between the gold-standard (i.e. fully sampled case) and the reconstructed MR images. The explicit definition of these metrics can be found in appendix Appendix B. To demonstrate the potentials of our approach, we designed a two-part evaluation protocol. The first one consists in comparing our technique versus three top reference techniques

<footnote>‡ Raw data and method to create the undersampled datasets can be downloaded here https://github.com/tschmoderer/2019_mri_reconstruction/blob/master/src/data/data_mri.tar.xz</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: (Top column) Visual samples of the datasets used in our experiments. (Bottom row) visualisation of some undersampling patterns used in our experiments using acceleration factor={2x,4x,6x,8x}.</figcaption></figure>

for MRI reconstructions: i) *zero-filling* (ZF) approach i.e. an inverse Fourier transform on the undersampled data ii) the classic *compressed sensing* (CS) algorithm [43] and iii) the *low rank and sparsity* (L+S) technique of [54]. The second comparison scheme is based on techniques that perform also two tasks simultaneously: *compressed sensing and motion* (CS+M) [3] and the recent method of (JPDAL) [92]. For visualisation purposes, we display a cropped version of the output reconstruction to inspect better the details.

## 5.2. Parameter Selection & Dictionary Setup

The parameters were set based on an empirical testing in a coarse to fine search strategy by maximising the SSIM metric on the dataset CINE and then we kept the same parameters for all experiments in order to show the inherent strength of the method and not its dependency to parameters. These parameters are fixed during the optimisation iterations and are stable for all the experiments conducted. The parameters $τ$ and $σ$ in the Chambolle & Pock iteration are chosen to fulfil the constraint $τσ||C||^2$ at each step. This is done by exactly computing the norm of the operator $C$ when possible, or by finding the largest set of $τ$ and $σ$ satisfying the constraint by numerical trials. From this, we set the parameter for our experiments as displayed in Table 1.

For each dataset, several experiments are performed depending on the utilisation of the dictionary. First, for each dataset and acceleration factor, we train a dictionary based on a $TV - L^1$ approximation of the optical flow computed on the inverse Fourier transform of undersampled data. In a second step, to show the transfer learning capability of our model, we train the dictionary on a $TV - L^1$ approximation of the optical flow from fully sampled phantom data, and use it for real cine cardiac applications.
<table><caption>Table 1: Parameter values used in our experiments.</caption><thead><tr><th>WEIGHING<br>PARAMETERS eq. (10)</th><th>IMAGE<br>RECONSTRUCTION<br>(algorithm 1 for eq. (21))</th><th>SPARSE<br>APPROXIMATION<br>(algorithm 1 with eqs. (23) and (24))</th><th>OPTICAL FLOW<br>APPROXIMATION<br>eqs. (25) and (26)</th><th>DICTIONARY<br>LEARNING<br>eq. (15)</th></tr></thead><tbody><tr><td>λ<sub>1</sub> = 0.003</td><td></td><td></td><td></td><td></td></tr><tr><td>λ<sub>2</sub> = 0.0001</td><td></td><td></td><td></td><td></td></tr><tr><td>λ<sub>3</sub> = 0.001</td><td>σ = 0.05</td><td>σ = 0.99</td><td>σ = 0.5</td><td>ε = 0.7N<sub>d</sub></td></tr><tr><td>λ<sub>4</sub> = 0.001</td><td>τ = 0.05</td><td>τ = 0.99</td><td>τ = 0.25</td><td>τ<sub>1</sub> in (18)</td></tr><tr><td>λ<sub>5</sub> = 0.001</td><td>θ = 1</td><td>θ = 1</td><td>θ = 1</td><td>τ<sub>2</sub> in (19)</td></tr><tr><td>λ<sub>6</sub> = 0.0001</td><td></td><td></td><td></td><td></td></tr><tr><td>ε = 0.001</td><td></td><td></td><td></td><td></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: (A) Samples extracted from our learned dictionary with 1024 atoms and patches of size 16 × 16. (B) Evolution of the energy during the learning process.</figcaption></figure>

We remark that one of the strengths of our approach is to avoid the need to have a ground truth to construct the dictionary. This is a highly desirable property as in the medical domain, the assumption of having ground truth is strong. We use over-complete patch-based dictionary using the TV − L¹ OF approximation. The dictionary training process involves the computation of the sparse decomposition of the reference flow, so we use this procedure to initialise our sequence a⁰. Our dictionary is composed of N_d = 1024 atoms, but for visualisation purposes we display 64 atoms illustrated on the left of fig. 2. The samples are coloured according to the standard visualisation of optical flow. The plot on the right of fig. 2 displays the evolution of the energy through the learning process to illustrate the speed and stability of convergence. We remark a fast convergence towards a dictionary D and a sparse decomposition a. This translates into a small computational time required to get the optimal dictionary and sparse decomposition.
<figure><img src="image_4.png" /><figcaption>Figure 3: Visual comparison of image quality of our model (MRIR-DLMC) against variational state of the art reconstruction techniques namely: zero filling (ZF), compressed sensing [43] (CS), and low rank and sparse decomposition [54], along with the gold standard fully-sampled reconstruction for two real datasets I and IV. For each dataset, the results are shown for three acceleration factors 2x, 4x and 8x.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Visual comparison of both reconstruction and optical flow estimation quality of our model (OURS) against two joint models namely the purely variational compressed sensing and optical flow model [9, 4] (CS+M) and the technique introduced in [92] (JPDAL), along with the gold standard fully-sampled reconstruction for datasets II (cardiac phantom) and IV (liver MRI). For each dataset, the results are shown for three acceleration factors 2x, 4x and 8x.</figcaption></figure>

## 5.3. Results & Discussion

In this section, we extensively evaluate our approach following the evaluation protocol stated in Section 4.1.

**Comparison against single task MRI reconstruction methods.** In fig. 3 and following common protocol for MRI reconstruction evaluation, we report the reconstructed images obtained for two real datasets I and IV and three acceleration factors 2x, 4x, and 8x, for a visual comparison of our model against variational MRI reconstruction techniques: zero-filling (ZF), compressed sensing [43] (CS) and low-rank and sparse decomposition [22] (L+S). We remark that the comparison between our
<figure><img src="image_6.png" /><figcaption>Figure 5: Reconstruction outputs under several acceleration counts along with time-lines (see white vertical and horizontal lines) references. Bottom part of each reconstruction displays the reconstruction outputs for the selected time-lines.</figcaption></figure>

technique and no motion can be observed when comparing the results from the CS reconstructions against ours. Also, the impact of the dictionary learning term can be measured by comparing the CS+M results and ours. As expected, the zero-filling method performs the worst in all cases and exhibit blurring and wave-like artefacts. We remark that our method outperforms the CS and L+S techniques at all acceleration factors. For an acceleration factor of 2x, this is particularly visible in the region beneath the ventricular chambers of the heart in which both CS and L+S reconstructions still exhibit wave-like and blurring artefacts contrary to ours. As the acceleration factor increases, the difference is even more visible as our method reduces the blurring and wave-like artefacts, better preserves fine details and exhibits sharper edges. Especially, the white blood vessels of the liver are better recovered with our technique than with the other presented ones. The edges of the right and left ventricular chambers of the heart as well as of the liver are sharper with our algorithm than with purely reconstruction methods specially for an acceleration factor of 8x. The blurring and wave-like artefacts in the liver and the myocardium are reduced with our technique for acceleration factors 4x and 8x compared to other reconstruction methods. All these statements are in favour of incorporating motion information in the reconstruction process of dynamic sequences to improve the overall image quality. To further support our technique, in Figure 5 we included time-line visualisations, for different reconstruction counts vs the gold-standard, for two different line directions corresponding to the white lines. The time-lines show the systolic and diastolic cardiac phases over the ventricles. The visualisations help to identify subtle missing changes and/or introduced artefacts.

**Comparison against joint MRI reconstruction and motion estimation models.** We now compare our method with other existing joint models addressing simultaneously MRI reconstruction and motion estimation namely a purely variational compressed sensing and optical flow model [9, 4], named CS+M, to evidence the improvement induced by the dictionary learning process, and the new compressed sensing and affine optical flow model introduced in [92] called JPDAL. Both the
reconstruction and the motion estimation are reported in fig. 4 for datasets II and IV and for three acceleration factors 2x, 4x and 8x, along with the gold standard fully-sampled reconstruction. Regarding the motion maps, we observe that the JPDAL retrieves small displacements mainly located on the boundary of the organs whereas our model is able to estimate bigger displacements located in the whole organ which seems to be closer to the physiological motion involved. Also, by visual comparison, adding the dictionary learning process tends to smooth and to better identify small moving regions compared to the purely variational model. The motion estimation improvement is reflected in higher quality and visually more pleasing reconstructed images especially for large acceleration factors. For acceleration factors 2x and 4x, the reconstruction for all the methods are visually comparable and competitive. For an acceleration factor of 8x, the CS+M reconstructions exhibit blurring artefacts especially in the liver and the left ventricular chamber of the heart which are less visible in our results. This is probably due to residual movements and less accurate motion estimation and shows the benefits of adding the dictionary learning process in the model to improve image quality. Small details such as the white blood vessels are also better recovered with our technique than with the CS+M one. The JPDAL reconstructions for an acceleration factor of 8x are also more blurred than ours and the homogeneous regions of the left ventricular chamber and the above area are better recovered with our algorithm than with the JPDAL model.

**Global Evaluation Performance.** We begin by reporting visual outputs of our proposed algorithmic technique versus single and other joint models along with the gold-standard (fully-sampled case). The results are displayed in fig. 6 using a high undersampling factor 6x. By visual inspection, we observe that the outputs generated by the single task techniques (i.e. ZF, CS and LS) are the ones displaying more artefacts than ours, and in particular, as expected, the worst performance is reported by ZF. This in terms of blurring artefacts and loss of preservation of relevant anatomy parts - for example in the cardiac datasets this can be observed around the papillary muscle of the heart whilst in the last dataset in the interior of the liver lobe. These observations are further supported by the corresponding difference maps, in which our technique reports closer reconstruction to the gold-standard than the compared single task techniques. For further support of our proposed technique, we also display selected outputs of our approach versus other techniques with similar philosophy as ours (i.e. joint models). A closer look at these outputs we can see that while CS+M and JDPAL perform better than the single task techniques, our proposed approach readily competes with these techniques, and in several cases, outperforms them. To further support our technique, we ran a nonparametric statistical test (Wilcoxon test) for two group comparison. From the test we found that our technique is statistically significantly different than other techniques including JPDAL (significance level = 0.05; p-value= 0.024). Moreover, we also reported a comparison in terms of computational time (in seconds), the results are displayed in Figure 6. These plots compare the computational time for all datasets and different undersampling factors for our technique and other multi-task models. From an
<table><caption>Table 2: Numerical comparison of our technique vs single and joint technique for different acceleration factors. The results are reported as the average of the corresponding metric over all the corresponding dataset.</caption><thead><tr><th rowspan="2">DATASET</th><th rowspan="2">RECONSTRUCTION SCHEME</th><th colspan="2">2x</th><th colspan="2">4x</th><th colspan="2">6x</th><th colspan="2">8x</th></tr><tr><th>PSNR</th><th>SSIM</th><th>PSNR</th><th>SSIM</th><th>PSNR</th><th>SSIM</th><th>PSNR</th><th>SSIM</th></tr></thead><tbody><tr><td rowspan="6">Dataset I</td><td>Zero-Filling</td><td>31.40</td><td>90.19</td><td>25.14</td><td>81.40</td><td>25</td><td>76.97</td><td>22.99</td><td>72.71</td></tr><tr><td>CS</td><td>32.85</td><td>93.58</td><td>31.57</td><td>88.57</td><td>27.98</td><td>81.48</td><td>22.97</td><td>72.76</td></tr><tr><td>L+S</td><td>34.21</td><td>92.77</td><td>31.09</td><td>86.20</td><td>27.74</td><td>80.10</td><td>22.98</td><td>72.68</td></tr><tr><td>CS+M</td><td>36.72</td><td>96.23</td><td>31.53</td><td>90.26</td><td><strong>28.39</strong></td><td>80.17</td><td><strong>24.9</strong></td><td>72.68</td></tr><tr><td>MC+JPDAL</td><td>36.7</td><td><strong>97.85</strong></td><td><strong>32.72</strong></td><td>92.06</td><td>27.80</td><td>84.29</td><td>23.15</td><td>75.51</td></tr><tr><td>MRIR-DLMC</td><td><strong>38.01</strong></td><td>97.33</td><td>32.35</td><td><strong>92.26</strong></td><td>27.65</td><td><strong>84.76</strong></td><td>23.12</td><td><strong>76.03</strong></td></tr><tr><td rowspan="6">Dataset IV</td><td>Zero-Filling</td><td>22.43</td><td>72.95</td><td>17.84</td><td>57.66</td><td>18.49</td><td>50.80</td><td>16.84</td><td>44.98</td></tr><tr><td>CS</td><td>26.61</td><td>83.48</td><td>22.7</td><td>69.31</td><td>19.36</td><td>54.15</td><td>17.03</td><td>45.86</td></tr><tr><td>L+S</td><td>24.26</td><td>77.68</td><td>21.71</td><td>63.39</td><td>19.19</td><td>52.46</td><td>16.82</td><td>44.96</td></tr><tr><td>CS+M</td><td>31.91</td><td>91.73</td><td>25.86</td><td>76.32</td><td>21.14</td><td>58.37</td><td>20.01</td><td>49.17</td></tr><tr><td>MC+JPDAL</td><td>32.81</td><td>93.07</td><td>27.3</td><td>82.01</td><td>22.34</td><td>65.62</td><td>19.39</td><td>54.1</td></tr><tr><td>MRIR-DLMC</td><td><strong>34.21</strong></td><td><strong>94.16</strong></td><td><strong>28.28</strong></td><td><strong>84.26</strong></td><td><strong>22.65</strong></td><td><strong>67.17</strong></td><td><strong>19.57</strong></td><td><strong>55.02</strong></td></tr><tr><td rowspan="6">Dataset II</td><td>Zero-Filling</td><td>28.27</td><td>78.87</td><td>23.97</td><td>68.08</td><td>20.53</td><td>61.39</td><td>22.57</td><td>61.09</td></tr><tr><td>CS</td><td>31.91</td><td>83.11</td><td>28.58</td><td>76.35</td><td>24.83</td><td>67.01</td><td>22.58</td><td>61.12</td></tr><tr><td>L+S</td><td>30.59</td><td>82.17</td><td>27.74</td><td>73.7</td><td>24.63</td><td>65.35</td><td>22.58</td><td>60.91</td></tr><tr><td>CS+M</td><td>32.65</td><td>87.1</td><td>30.5</td><td>81.91</td><td>24.84</td><td>70.30</td><td>21.94</td><td>63.06</td></tr><tr><td>MC+JPDAL</td><td><strong>36.5</strong></td><td><strong>93</strong></td><td><strong>32.45</strong></td><td><strong>86.44</strong></td><td><strong>28.08</strong></td><td><strong>76.92</strong></td><td>24.46</td><td>69.45</td></tr><tr><td>MRIR-DLMC</td><td>33.19</td><td>88.46</td><td>30.5</td><td>82.77</td><td>27.76</td><td>76.18</td><td><strong>25.14</strong></td><td><strong>71.9</strong></td></tr><tr><td rowspan="6">Dataset III</td><td>Zero-Filling</td><td>28.84</td><td>79.3</td><td>24.31</td><td>67.84</td><td>20.57</td><td>60.48</td><td>22.46</td><td>59.61</td></tr><tr><td>CS</td><td>31.19</td><td>82.51</td><td>28.66</td><td>74.89</td><td>25.26</td><td>65.87</td><td>22.48</td><td>59.61</td></tr><tr><td>L+S</td><td>30.59</td><td>82.8</td><td>27.99</td><td>72.93</td><td>24.9</td><td>64.5</td><td>22.46</td><td>59.38</td></tr><tr><td>CS+M</td><td>32.11</td><td>86.13</td><td>30.02</td><td><strong>89.16</strong></td><td>24.28</td><td>67.74</td><td>21.55</td><td>60.59</td></tr><tr><td>MC+JPDAL</td><td><strong>35.37</strong></td><td><strong>91.86</strong></td><td><strong>31.24</strong></td><td>83.85</td><td><strong>27.3</strong></td><td><strong>74.16</strong></td><td><strong>24.14</strong></td><td><strong>66.66</strong></td></tr><tr><td>MRIR-DLMC</td><td>32.76</td><td>87.48</td><td>29.41</td><td>80.12</td><td>26.9</td><td>73.43</td><td>24.05</td><td>66.5</td></tr></tbody></table>

inspection to these plots, one can see that our method not only improves performance wise but also demands lower computational time that other existing techniques following similar philosophy.

For a more detailed quantitative analysis, we report the global results in Table 2 for all compared techniques and for acceleration factor of 2x, 4x, 6x and 8x. The displayed results are the average of the image metrics, i.e. PSNR and SSIM, across the full corresponding dataset. With respect to the realistic datasets, we reported the overall best performance whilst for the simulated datasets our method performs similarly well, and places second only behind JPDAL. We further support our technique by computing the computational time (in seconds) of our technique and other multi-task techniques: CS+M and MC+JPDAL. The results are reported in Figure 7. In a closer inspection at this figure, one can observe that our technique required for all cases less computational
<figure><img src="image_7.png" /><figcaption>Figure 6: Visual comparison of our technique versus single and other joint models. The results are displayed using 4x and 6x acceleration factor. A closer look at the difference maps, one can see that our technique reports reconstructions closer to the gold-standard with respect to single task techniques whilst readily compete with, and several times outperforms, the other compared joint models.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Computational time comparison of our approach versus other multi-task techniques reported in seconds. The results are reported from several acceleration factors (2x, 4x, 6x and 8x) and for all datasets.</figcaption></figure>

<table><caption>Table 3: Numerical comparison of our technique vs other reconstruction methods. The numerical values are computed as the averages of the similarity metrics over the complete corresponding dataset. w/TL denotes the transfer learning capability of our technique, that is- the results are from training our dictionary with phantom datasets and applied to the real cardiac cine.</caption><thead><tr><th rowspan="2">DATASET</th><th rowspan="2">RECONSTRUCTION SCHEME</th><th colspan="2">2x</th><th colspan="2">4x</th><th colspan="2">6x</th><th colspan="2">8x</th></tr><tr><th>PSNR</th><th>SSIM</th><th>PSNR</th><th>SSIM</th><th>PSNR</th><th>SSIM</th><th>PSNR</th><th>SSIM</th></tr></thead><tbody><tr><td rowspan="6">CARDIAC CINE</td><td>Zero-Filling</td><td>31.40</td><td>90.19</td><td>25.14</td><td>81.40</td><td>25</td><td>76.97</td><td>22.99</td><td>72.71</td></tr><tr><td>CS</td><td>32.85</td><td>93.58</td><td>31.57</td><td>88.57</td><td>27.98</td><td>81.48</td><td>22.97</td><td>72.76</td></tr><tr><td>L+S</td><td>34.21</td><td>92.77</td><td>31.09</td><td>86.20</td><td>27.74</td><td>80.10</td><td>22.98</td><td>72.68</td></tr><tr><td>CS+M</td><td>36.72</td><td>96.23</td><td>31.53</td><td>90.26</td><td><strong>28.39</strong></td><td>80.17</td><td><strong>24.9</strong></td><td>72.68</td></tr><tr><td>MC+JPDAL</td><td>36.7</td><td><strong>97.85</strong></td><td><strong>32.72</strong></td><td>92.06</td><td>27.80</td><td>84.29</td><td>23.15</td><td>75.51</td></tr><tr><td>MRIR-DLMC</td><td><strong>38.01</strong></td><td>97.33</td><td>32.35</td><td><strong>92.26</strong></td><td>27.65</td><td>84.76</td><td>23.12</td><td>76.03</td></tr><tr><td></td><td>RIR-DLMC w/TL</td><td>37</td><td>97.05</td><td>31.79</td><td>90.71</td><td>27.70</td><td><strong>84.85</strong></td><td>24.80</td><td><strong>81.01</strong></td></tr></tbody></table>

resources than the compared techniques.

## 5.4. Transfer Learning Capabilities of our Technique

We have shown in the previous section the advantages of our technique in terms of quality and computational performance. Besides these benefits, robustness of our training process when no ground-truth is available has been evidenced which is a highly desirable property. In this section, we highlight another strength of our technique, which is the ability of learning over phantom datasets and transfer this knowledge to real datasets. The ability to transfer knowledge from phantom to real data is very challenging yet highly desirable for any algorithmic technique. As pointed out in [62], this capability is not straightforward and such generalisation might be difficult to achieve due to the overparametrisation of the model given by the fundamental gap between synthetic and real data. To prove this capability, we show a numerical comparison in table 3 and
denote the new results as 'MRIR-DLMC w/TL'. From the results, one can see that firstly, our approach is capable to perform similarly well, when trained with phantoms datasets, than JPDAL and ours (trained with realistic dataset). Our 'MRIR-DLMC w/TL' readily compete numerically with the best scores and in some cases outperforms the compared techniques. This capability opens several opportunities in the medical domain as we avoid the need to retrain our dictionary.

Overall, our proposed model considers local and global neighbourhood (patch based process) flow relations, we encode the regularisation of the motions and discontinuities through a sparse representation rather than in a separately manner, and by learning through dictionaries we can characterise complex patterns and structures. Our technique offers three major advantages over the compared techniques i) it improves performance wise over the compared techniques, and in particular, over JPDAL with a statistical difference, ii) it requires less computational time than the other techniques and iii) our technique can narrow the gap between synthetic and real data- as our technique performs well when trained on phantoms and then applied to real data.

## 6. Conclusion

This work addresses the problem of reconstructing high-quality images from dynamic under-sampled MRI measurements. This is a fundamental problem in MRI due to the long acquisition time required to form an image. During scanning, involuntary and breathing motion often lead to image degradation, compromising clinical interpretation and relevance for diagnosis.

To tackle this problem, we propose a novel mathematical model to improve the reconstruction quality by estimating motion in dynamic MRI. The underlying idea is to exploit the strong correlation between the two tasks of image reconstruction and motion estimation. The main motivation of considering motion as a second task is because, in a dynamic setting, inherent motion is contained in the scene and therefore it strongly depends on the reconstruction quality. Our research hypothesis is that the reconstruction accuracy can be substantially improved by a good approximation of the motion scene. Our proposed approach combines, in a multi-task and hybrid model, the traditional compressed sensing formulation for the reconstruction of dynamic MRI and motion compensation by learnt optical flow approximation. First, we introduce the classical compressed sensing reconstruction formulation along with motion estimation algorithms depicted in [3]. Second, we show how the optical flow can be learnt from reference samples following the work introduced in [32]. More precisely, our optical flow model is computed as a sparse linear combination of basis functions from a learnt dictionary. Then, we embed in a single functional the reconstruction process and the optical flow estimation by dictionary learning. We present in details our efficient and tractable optimisation framework adopted to solve the non-convex problem based on an alternating splitting scheme. Finally, we extensively demonstrate the potential of our proposed model in the context of dynamic MRI through a set of numerical experiments.
In particular, we present comparisons with single task MRI reconstructions as well as with joint MRI reconstruction and motion estimation models, demonstrating that our approach reaches and outperforms state-of-the-art among mathematical variational models and joint models.

Our multi-task hybrid model has shown enormous potential for improving the quality of dynamic MRI reconstruction from highly undersampled data by carefully intertwining two imaging tasks that are traditionally performed separately. Furthermore, our motion estimation is based on a learnt approach that relies on a sparse representation in a dictionary of reference samples, with the great advantage that ground truth data is not required. In addition, we explore the potential of our proposed method in the context of transfer learning, that is we learn our dictionary on phantom data and apply it to real dataset. Our results are promising especially in the medical domain, where many times we deal with a severe lack of ground truth data. A perspective for future research is to embed the proposed algorithm in a coarse to fine pyramidal approach. We would also add to the optical flow constraint another term coding the possible brightness variations in the sequence. Finally, we will take our findings and apply them to a more clinical study. Overall, our technique has three major advantages over other techniques. Firstly, our technique statistically improves reconstruction quality over the compared techniques. Secondly, our model requires less computational resources than other multi-task models. Finally, we show the ability of our technique to transfer knowledge from phantom to real data.

## Appendix A. Description of the data structures and the operators

We assume that the space-time discrete grid consists of the following set of points:

$$ \{(i, j, t) : i = 1, \dots, N_x, j = 1, \dots, N_y, t = 1, \dots, N_t\}. $$

Therefore, the reconstructed image sequence $\mathbf{m}$ takes values in $\mathbb{R}^{N_x N_y N_t}$, whereas the components $u_x$, $u_y$ of the optical flow take values in $\mathbb{R}^{N_x N_y (N_t - 1)}$. We now discuss the discretisation of the operators. Assume we have a sequence $\mathbf{m}$ of $N_t$ frames of size $N_x \times N_y$ pixels. We use a forward discretisation for the temporal derivative and a central discretisation for the spatial derivative. Neumann boundary conditions are applied:

$$ \mathbf{m}_t(i, j, t) = \begin{cases} \mathbf{m}(i, j, t+1) - \mathbf{m}(i, j, t) & \text{if } t < N_t, \\ 0 & \text{else,} \end{cases} $$

$$ \mathbf{m}_x(i, j, t) = \frac{1}{2} \begin{cases} \mathbf{m}(i+1, j, t) - \mathbf{m}(i-1, j, t) & \text{if } i > 1 \text{ and } i < N_x \text{ and } t < N_t, \\ 0 & \text{else,} \end{cases} $$

$$ \mathbf{m}_y(i, j, t) = \frac{1}{2} \begin{cases} \mathbf{m}(i, j+1, t) - \mathbf{m}(i, j-1, t) & \text{if } j > 1 \text{ and } j < N_y \text{ and } t < N_t, \\ 0 & \text{else.} \end{cases} $$
The adjoint operator then yields:

$$
\begin{align*}
\mathbf{y}_t(i, j, t) &= - \begin{cases} \mathbf{y}(i, j, t) & \text{if } t = 1, \\ \mathbf{y}(i, j, t) - \mathbf{y}(i, j, t-1) & \text{if } t > 1 \text{ and } t < N_t, \\ -\mathbf{y}(i, j, t-1) & \text{if } t = N_t, \end{cases} \\
\mathbf{y}_x(i, j, t) &= \frac{-1}{2} \begin{cases} \mathbf{y}(i+1, j, y) & \text{if } i \le 2 \text{ and } t < N_t \\ \mathbf{y}(i-1, j, y) - \mathbf{y}(i+1, j, y) & \text{if } i > 2 \text{ and } i < N_x - 1 \text{ and } t < N_t \\ -\mathbf{y}(i-1, j, y) & \text{if } i \ge N_x - 1 \text{ and } t < N_t \\ 0 & \text{else,} \end{cases} \\
\mathbf{y}_y(i, j, t) &= \frac{-1}{2} \begin{cases} \mathbf{y}(i, j+1, y) & \text{if } j \le 2 \text{ and } t < N_t \\ \mathbf{y}(i, j-1, y) - \mathbf{y}(i, j+1, y) & \text{if } j > 2 \text{ and } j < N_y - 1 \text{ and } t < N_t \\ -\mathbf{y}(i, j-1, y) & \text{if } j \ge N_y - 1 \text{ and } t < N_t \\ 0 & \text{else.} \end{cases}
\end{align*}
$$

The discrete derivatives for the optical flow vector field are calculated using forward
differences and Neumann boundary conditions. The corresponding adjoint operator
consists of backward differences with Dirichlet boundary conditions and is applied to
the dual variable **y**. In the following, *u* denotes either the horizontal component **u**<sub>x</sub> or
the vertical **u**<sub>y</sub>.

$$
\begin{align*}
u_x(i, j, t) &= \begin{cases} u(i+1, j, t) - u(i, j, t) & \text{if } i < N_x \\ 0 & \text{if } i = N_x, \end{cases} \\
u_y(i, j, t) &= \begin{cases} u(i, j+1, t) - u(i, j, t) & \text{if } j < N_y \\ 0 & \text{if } j = N_y, \end{cases} \\
\nabla \cdot \mathbf{y}(i, j, t) &= \begin{cases} \mathbf{y}^1(i, j, t) & \text{if } i = 1 \\ \mathbf{y}^1(i, j, t) - \mathbf{y}^1(i-1, j, t) & \text{if } i > 1 \text{ and } i < N_x \\ -\mathbf{y}^1(i-1, j, t) & \text{if } i = N_x \end{cases} \\
&\quad + \begin{cases} \mathbf{y}^2(i, j, t) & \text{if } j = 1 \\ \mathbf{y}^2(i, j, t) - \mathbf{y}^2(i, j-1, t) & \text{if } j > 1 \text{ and } j < N_y \\ -\mathbf{y}^2(i, j-1, t) & \text{if } j = N_y. \end{cases}
\end{align*}
$$

**Appendix B. Error measures**

We evaluate the performance of our model in terms of quality for the reconstructed
images. The main metric thus used is the structural similarity index (SSIM) [82] which
measures the differences in luminance, contrast and structure of the ground truth image
m and the reconstruction m<sub>r</sub> as follows,

$$
SSIM := \frac{(2\mu_m\mu_{m_r} + C_1)(2\sigma_{m,m_r} + C_2)}{(\mu_m^2 + \mu_{m_r}^2 + C_1)(\sigma_m^2 + \sigma_{m_r}^2 + C_2)}
$$

where $\mu_m, \mu_{m_r}, \sigma_m, \sigma_{m_r}$ and $\sigma_{m,m_r}$ are local means, standard deviations and cross covariance for ground truth image $m$ and reconstruction $m_r$ respectively. The constants
are fixed to $C_1 = 0.01^2$ and $C_2 = 0.03^2$. The SSIM takes values between -1 and 1 where 1 stands for perfect similarity. Moreover we calculate the peak signal-to-noise ratio (PSNR) between the ground truth and the reconstruction given by

$$PSNR := 10 \log_{10} \left( \frac{\max(m^2)}{\text{mean}((m - m_r)^2)} \right). \qquad (B.1)$$

## Appendix C. Proximal operators

In this section we present the computation of the proximal operator of the optical flow minimisation for the Chambolle and Pock algorithm [13]. We recall the unconstrained problem,

$$\min_{\mathbf{u}} \int_{0}^{T} \lambda_{3} \left\| \frac{\partial \mathbf{m}}{\partial t} + \nabla \mathbf{m} \cdot \mathbf{u} \right\|_{1} + \lambda_{4} \|\nabla \mathbf{u}\|_{1} + \lambda_{5} \sum_{\mathbf{p} \in \mathcal{P}} \|R_{\mathbf{p}}\mathbf{u} - D\mathbf{a}_{\mathbf{p}}\|_{F}^{2} \, dt.$$

To simplify the notation, we introduce the affine operator $\rho(\mathbf{u}) = \frac{\partial \mathbf{m}}{\partial t} + \nabla \mathbf{m} \cdot \mathbf{u}$. The challenging part is to compute the resolvent operator of $G(\mathbf{u}) = \lambda_3 \|\rho(\mathbf{u})\|_1 + \lambda_5 \sum_{\mathbf{p} \in \mathcal{P}} \|R_{\mathbf{p}}\mathbf{u} - D\mathbf{a}_{\mathbf{p}}\|_F^2$. We have,

$$ (I + \tau \partial G)^{-1}(\mathbf{u}) = \arg \min_{\mathbf{v}} \left\{ \frac{\|\mathbf{v} - \mathbf{u}\|^2}{2\tau} + G(\mathbf{v}) \right\} \qquad (C.1) \\ = \arg \min_{\mathbf{v}} \left\{ \frac{\|\mathbf{v} - \mathbf{u}\|^2}{2\tau} + \lambda_3 \|\rho(\mathbf{v})\|_1 + \lambda_5 \sum_{\mathbf{p} \in \mathcal{P}} \|R_{\mathbf{p}}\mathbf{v} - D\mathbf{a}_{\mathbf{p}}\|_F^2 \right\}. $$

The problem is decomposed in three cases: $\rho(\mathbf{u}) > 0$, $\rho(\mathbf{u}) < 0$ and $\rho(\mathbf{u}) = 0$, the previous constraint are intended component-wise. In the first two cases the $l^1$-norm is differentiable so they can be treated in a similar fashion:

$$ \rho(\mathbf{v}) > 0 \implies \frac{\mathbf{v} - \mathbf{u}}{\tau} + \lambda_3 \nabla \mathbf{m} + 2\lambda_5 \sum R_{\mathbf{p}}^T (R_{\mathbf{p}} \mathbf{v} - D\mathbf{a}_{\mathbf{p}}) = 0 \\ \implies (I + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T R_{\mathbf{p}}) \mathbf{v} = \mathbf{u} - \tau\lambda_3 \nabla \mathbf{m} + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T D\mathbf{a}_{\mathbf{p}}. \qquad (C.2) $$

The equality eq. (C.2) gives the solution of eq. (C.1). Like in eq. (25) we note the operator $(I + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T R_{\mathbf{p}})$ by $\mathcal{A}$. We then plug the solution eq. (C.2) in the condition $\rho(\mathbf{v}) > 0$ to obtain a necessary condition for the realisation of this solution:

$$ \rho \left( \mathcal{A}^{-1} (\mathbf{u} + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T D\mathbf{a}_{\mathbf{p}}) \right) - \tau\lambda_3 \mathcal{A}^{-1} \|\nabla \mathbf{m}\|^2 > 0. \qquad (C.3) $$

When $\rho(\mathbf{v}) = 0$ the $l^1$-norm is no longer differentiable so we have to consider its sub-differential.

$$ \rho(\mathbf{v}) = 0 \implies \frac{\mathbf{v} - \mathbf{u}}{\tau} + \lambda_3 [-\nabla \mathbf{m}; \nabla \mathbf{m}] + 2\lambda_5 \sum R_{\mathbf{p}}^T (R_{\mathbf{p}} \mathbf{v} - D\mathbf{a}_{\mathbf{p}}) = 0 \\ \implies (I + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T R_{\mathbf{p}}) \mathbf{v} = \mathbf{u} - \tau\lambda_3 [-\nabla \mathbf{m}; \nabla \mathbf{m}] + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T D\mathbf{a}_{\mathbf{p}}. \qquad (C.4) $$
In the previous eq. (C.4) we note by $[-\nabla \mathbf{m}; \nabla \mathbf{m}]$ any element of the sub-differential of the $l^1$-norm, then the condition $\rho(\mathbf{v}) = 0$ gives us the required element:

$$ \rho \left( \mathcal{A}^{-1}(\mathbf{u} + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T D\mathbf{a}_{\mathbf{p}}) \right) - \tau\lambda_3[-\nabla \mathbf{m}; \nabla \mathbf{m}] \cdot \nabla \mathbf{m} = 0 \\ \implies \tau\lambda_3[-\nabla \mathbf{m}; \nabla \mathbf{m}] = \rho \left( \mathcal{A}^{-1}(\mathbf{u} + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T D\mathbf{a}_{\mathbf{p}}) \right) \frac{\nabla \mathbf{m}}{\|\nabla \mathbf{m}\|^2}. \quad (C.5) $$

Finally,

$$ \mathcal{A}(I + \tau \partial G)^{-1}(\mathbf{u}) = \tilde{\mathbf{u}} + \begin{cases} -\tau\lambda_3 \nabla \mathbf{m} & \text{if } \rho(\mathcal{A}^{-1}\tilde{\mathbf{u}}) > \tau\lambda_3 \mathcal{A}^{-1} \|\nabla \mathbf{m}\|^2 \\ \tau\lambda_3 \nabla \mathbf{m} & \text{if } \rho(\mathcal{A}^{-1}\tilde{\mathbf{u}}) < -\tau\lambda_3 \mathcal{A}^{-1} \|\nabla \mathbf{m}\|^2 \\ -\tau\lambda_3 \rho(\mathcal{A}^{-1}\tilde{\mathbf{u}}) \frac{\nabla \mathbf{m}}{\|\nabla \mathbf{m}\|^2} & \text{otherwise} \end{cases} $$

where $\tilde{\mathbf{u}} = \mathbf{u} + 2\tau\lambda_5 \sum R_{\mathbf{p}}^T D\mathbf{a}_{\mathbf{p}}$.

## Appendix D. Supplementary Experiments

In this section, we extend the comparison results of Table 2 from the main paper. Table D1 shows a performance comparison of our technique under different initialisations (ZF, CS and CS+M). A closer look at the results shows that as better the initialisation as higher the reconstruction quality. Whilst using initialisation such as CS or CS+M indeed improves the performance, it also implies several drawbacks such as substantial computational load as additional optimisation process need to be solved, and biased reconstructions closer to the initialisation. Moreover, the problem at hand might be seen as being solved twice. For these reasons, we follow common protocol for initialisation and, in the main paper, we use zero filling as our initialisation.

<table><caption>Table D1: Performance comparison of our technique under different initialisations and different undersampling factors. The results are reported as the average of the corresponding metric on Dataset I.</caption><thead><tr><th>ACCELERATION FACTOR</th><th>INITIALISATION</th><th>SSIM</th><th>PSNR</th><th>INITIALISATION</th><th>SSIM</th><th>PSNR</th><th>INITIALISATION</th><th>SSIM</th><th>PSNR</th></tr></thead><tbody><tr><td>2x</td><td rowspan="4">ZF</td><td>97.33</td><td>38.01</td><td rowspan="4">CS</td><td>97.70</td><td>37.47</td><td rowspan="4">CS+M</td><td>98.34</td><td>37.75</td></tr><tr><td>4x</td><td>92.26</td><td>32.35</td><td>92.47</td><td>32.37</td><td>92.90</td><td>33.24</td></tr><tr><td>6x</td><td>84.76</td><td>27.65</td><td>84.68</td><td>28.07</td><td>84.92</td><td>29.05</td></tr><tr><td>8x</td><td>76.03</td><td>23.12</td><td>76.19</td><td>23.66</td><td>76.61</td><td>24.40</td></tr></tbody></table>

## Acknowledgments

AIAR gratefully acknowledges the financial support of the CMIH and CCIMI University of Cambridge. VC acknowledges the financial support of the Cancer Research UK. CBS acknowledges support from the Leverhulme Trust project on 'Breaking the non-convexity barrier', EPSRC grant Nr. EP/M00483X/1, EP/S026045/1, the EPSRC
Centre Nr. EP/N014588/1, Wellcome Trust project "All-in-one cancer imaging", the RISE projects CHiPS and NoMADS, the Cantab Capital Institute for the Mathematics of Information and the Alan Turing Institute.

[1] Omar Arif, Hammad Afzal, Haider Abbas, Muhammad Faisal Amjad, Jiafu Wan, and Raheel Nawaz. Accelerated Dynamic MRI Using Kernel-Based Low Rank Constraint. *Journal of medical systems*, 43(8):271, 2019.

[2] M Salman Asif, Lei Hamilton, Marijn Brummer, and Justin Romberg. Motion-adaptive spatio-temporal regularization for accelerated dynamic MRI. *Magnetic Resonance in Medicine*, 70(3):800–812, 2013.

[3] Angelica I Aviles-Rivero, Noémie Debroux, Guy Williams, Martin J Graves, and Carola-Bibiane Schönlieb. Compressed sensing plus motion (cs+ m): A new perspective for improving undersampled mr image reconstruction. *Medical Image Analysis*, 68:101933, 2021.

[4] Graves M Aviles-Rivero AI, Williams G and Schönlieb CB. CS+M: A simultaneous reconstruction and motion estimation approach for improving undersampled MRI reconstruction. In *In Proceedings of the 26th Annual Meeting ISMRM*, 2018.

[5] Chenglong Bao, Hui Ji, Yuhui Quan, and Zuowei Shen. Dictionary learning for sparse coding: Algorithms and convergence analysis. *IEEE transactions on pattern analysis and machine intelligence*, 38(7):1356–1369, 2015.

[6] Çağdaş Bilen, Ivan Selesnick, Yao Wang, Ricardo Otazo, and Daniel K Sodickson. A motion compensating prior for dynamic mri reconstruction using combination of compressed sensing and parallel imaging. In *2011 IEEE Signal Processing in Medicine and Biology Symposium (SPMB)*, pages 1–6. IEEE, 2011.

[7] Thomas Brox, Andrés Bruhn, Nils Papenberg, and Joachim Weickert. High accuracy optical flow estimation based on a theory for warping. In *European conference on computer vision*, pages 25–36. Springer, 2004.

[8] Andrés Bruhn, Joachim Weickert, and Christoph Schnörr. Lucas/kanade meets horn/schunck: Combining local and global optic flow methods. *International journal of computer vision*, 61(3):211–231, 2005.

[9] M. Burger, H. Dirks, and C. Schönlieb. A variational model for joint motion estimation and image reconstruction. *SIAM Journal on Imaging Sciences*, 11(1):94–128, 2018.

[10] Aurélien Bustin, Anne Menini, Martin A Janich, Darius Burschka, Jacques Felblinger, Anja CS Brau, and Freddy Odille. Motion Estimated-Compensated Reconstruction with Preserved-Features in Free-Breathing Cardiac MRI. In *Reconstruction, Segmentation, and Analysis of Medical Images*, pages 70–80. Springer, 2016.

[11] Emmanuel J Candes, Justin K Romberg, and Terence Tao. Stable signal recovery from incomplete and inaccurate measurements. *Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences*, 59(8):1207–1223, 2006.

[12] Rich Caruana. Multitask learning. *Machine learning*, 28(1):41–75, 1997.

[13] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. *Journal of Mathematical Imaging and Vision*, pages 120–145, 2011.

[14] V. Corona, A.I. Aviles-Rivero, N. Debroux, M.Grave, C. Le Guyader, C.-B. Schönlieb, and G. Williams. Multi-tasking to correct: motion-compensated MRI via joint reconstruction and registration. *Scale Space and Variational Methods in Computer Vision LNCS conference proceedings*, Springer., 2019.

[15] David L Donoho. Compressed sensing. *IEEE Transactions on information theory*, 52(4):1289–1306, 2006.

[16] Alexey Dosovitskiy, Philipp Fischer, Eddy Ilg, Philip Hausser, Caner Hazirbas, Vladimir Golkov, Patrick Van Der Smagt, Daniel Cremers, and Thomas Brox. Flownet: Learning optical flow with convolutional networks. In *IEEE international conference on computer vision (CVPR)*,
pages 2758–2766, 2015.

[17] Michael Elad and Michal Aharon. Image denoising via sparse and redundant representations over learned dictionaries. *IEEE Transactions on Image Processing*, 15(12):3736–3745, 2006.

[18] Li Feng, Monvadi B Srichai, Ruth P Lim, Alexis Harrison, Wilson King, Ganesh Adluru, Edward VR Dibella, Daniel K Sodickson, Ricardo Otazo, and Daniel Kim. Highly accelerated real-time cardiac cine MRI using k-t SPARSE-SENSE. *Magnetic Resonance in Medicine*, pages 64–74, 2013.

[19] M Filipovic, P-A Vuissoz, A Codreanu, M Claudon, and J Felblinger. Motion compensated generalized reconstruction for free-breathing dynamic contrast-enhanced MRI. *Magnetic Resonance in Medicine*, 65(3):812–822, 2011.

[20] Denis Fortun, Patrick Bouthemy, and Charles Kervrann. Optical flow modeling and computation: A survey. *Computer Vision and Image Understanding*, 134:1–21, 2015.

[21] Urs Gamper, Peter Boesiger, and Sebastian Kozerke. Compressed sensing in dynamic MRI. *Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine*, 59(2):365–373, 2008.

[22] Hao Gao, Stanislas Rapacchi, Da Wang, John Moriarty, Conor Meehan, James Sayre, Gerhard Laub, Paul Finn, and Peng Hu. Compressed sensing using prior rank, intensity and sparsity model (PRISM): applications in cardiac cine MRI. In *Proceedings of the Annual Meeting of ISMRM*, 2012.

[23] Nadine Gdaniec, Holger Eggers, Peter Börnert, Mariya Doneva, and Alfred Mertins. Robust abdominal imaging with incomplete breath-holds. *Magnetic Resonance in Medicine*, 71(5):1733–1742, 2014.

[24] Michael A Gennert and Shahriar Negahdaripour. Relaxing the brightness constancy assumption in computing optical flow. Technical report, MASSACHUSETTS INST OF TECH CAMBRIDGE ARTIFICIAL INTELLIGENCE LAB, 1987.

[25] Rohini George, Theodore D. Chung, Sastry S. Vedam, Viswanathan Ramakrishnan, Radhe Mohan, Elisabeth Weiss, and Paul J. Keall. Audio-visual biofeedback for respiratory-gated radiotherapy: Impact of audio instruction and audio-visual biofeedback on respiratory-gated radiotherapy. *International Journal of Radiation Oncology*Biology*Physics, 65(3):924–933, 2006.

[26] David Gibson and Michael Spann. Robust optical flow estimation based on a sparse motion trajectory set. *IEEE Transactions on Image Processing*, 12(4):431–445, 2003.

[27] Joel Gibson and Oge Marques. Sparsity in optical flow and trajectories. *Signal, Image and Video Processing*, 10(3):487–494, 2016.

[28] Melissa W Haskell, Stephen F Cauley, and Lawrence L Wald. TARgeted Motion Estimation and Reduction (TAMER): data consistency based motion mitigation for MRI using a reduced model joint optimization. *IEEE transactions on medical imaging*, 37(5):1253–1265, 2018.

[29] Jeremy D.P. Hoisak, Katharina E. Sixel, Romeo Tirona, Patrick C.F. Cheung, and Jean-Philippe Pignol. Prediction of lung tumour position based on spirometry and on abdominal displacement: Accuracy and reproducibility. *Radiotherapy and Oncology*, 78(3):339–346, 2006.

[30] Berthold KP Horn and Brian G Schunck. Determining optical flow. In *Techniques and Applications of Image Understanding*, volume 281, pages 319–331. International Society for Optics and Photonics, 1981.

[31] Eddy Ilg, Nikolaus Mayer, Tonmoy Saikia, Margret Keuper, Alexey Dosovitskiy, and Thomas Brox. Flownet 2.0: Evolution of optical flow estimation with deep networks. In *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 2462–2470, 2017.

[32] Kui Jia, Xiaogang Wang, and Xiaoou Tang. Optical flow estimation using learned sparse model. In *2011 International Conference on Computer Vision*, pages 2391–2398. IEEE, 2011.

[33] Steve B. Jiang. Technical aspects of image-guided respiration-gated radiation therapy. *Medical Dosimetry*, 31(2):141–151, 2006.

[34] Hong Jung, Kyunghyun Sung, Krishna S Nayak, Eung Yeop Kim, and Jong Chul Ye. k-t FOCUSS: a general compressed sensing framework for high resolution dynamic MRI. *Magnetic*
*Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine*, 61(1):103–116, 2009.

[35] Shuichiro Kaji, Philip C Yang, Adam B Kerr, WH Wilson Tang, Craig H Meyer, Albert Macovski, John M Pauly, Dwight G Nishimura, and Bob S Hu. Rapid evaluation of left ventricular volume and mass without breath-holding using real-time interactive cardiac magnetic resonance imaging system. *Journal of the American College of Cardiology*, 38(2):527–533, 2001.

[36] Ruoteng Li, Robby T Tan, Loong-Fah Cheong, Angelica I Aviles-Rivero, Qingnan Fan, and Carola-Bibiane Schönlieb. RainFlow: Optical flow under rain streaks and rain veiling effect. In *IEEE International Conference on Computer Vision (ICCV)*, pages 7304–7313, 2019.

[37] Zhi-Pei Liang. Spatiotemporal imaging with partially separable functions. In *IEEE International Symposium on Biomedical Imaging (ISBI)*, pages 988–991, 2007.

[38] Sajan Goud Lingala, Edward DiBella, and Mathews Jacob. Deformation corrected compressed sensing (DC-CS): a novel framework for accelerated dynamic MRI. *IEEE transactions on medical imaging*, 34(1):72–85, 2014.

[39] Sajan Goud Lingala, Yue Hu, Edward DiBella, and Mathews Jacob. Accelerated dynamic MRI exploiting sparsity and low-rank structure: kt SLR. *IEEE Transactions on Medical Imaging*, pages 1042–1054, 2011.

[40] Sajan Goud Lingala and Mathews Jacob. Blind compressive sensing dynamic mri. *IEEE transactions on medical imaging*, 32(6):1132–1145, 2013.

[41] Ce Liu, Jenny Yuen, and Antonio Torralba. Sift flow: Dense correspondence across scenes and its applications. *IEEE transactions on pattern analysis and machine intelligence*, 33(5):978–994, 2010.

[42] Jiulong Liu, Angelica I Aviles-Rivero, Hui Ji, and Carola-Bibiane Schönlieb. Rethinking medical image reconstruction via shape prior, going deeper and faster: Deep joint indirect registration and reconstruction. *arXiv preprint arXiv:1912.07648*, 2019.

[43] Michael Lustig, David Donoho, and John M. Pauly. Sparse MRI: The application of compressed sensing for rapid MR imaging. *Magnetic Resonance in Medicine*, 58(6):1182–1195, 2007.

[44] Michael Lustig and John M Pauly. SPIRIT: iterative self-consistent parallel imaging reconstruction from arbitrary k-space. *Magnetic Resonance in Medicine*, pages 457–471, 2010.

[45] Michael Lustig, Juan M Santos, David L Donoho, and John M Pauly. k-t SPARSE: High frame rate dynamic MRI exploiting spatio-temporal sparsity. In *Proceedings of the Annual Meeting of ISMRM*, volume 2420, 2006.

[46] Julien Mairal, Francis Bach, Jean Ponce, Guillermo Sapiro, and Andrew Zisserman. Non-local sparse models for image restoration. In *2009 IEEE 12th International Conference on Computer Vision*, pages 2272–2279. IEEE, 2009.

[47] Angshul Majumdar. *Compressed sensing for magnetic resonance image reconstruction*. Cambridge University Press, 2015.

[48] Angshul Majumdar and Rabab K Ward. Exploiting rank deficiency and transform domain sparsity for MR image reconstruction. *Magnetic Resonance Imaging*, pages 9–18, 2012.

[49] Xin Miao, Sajan Goud Lingala, Yi Guo, Terrence Jao, Muhammad Usman, Claudia Prieto, and Krishna S Nayak. Accelerated cardiac cine MRI using locally low rank and finite difference constraints. *Magnetic Resonance Imaging*, 2016.

[50] Freddy Odille, Anne Menini, Jean-Marie Escanyé, Pierre-André Vuissoz, Pierre-Yves Marie, Marine Beaumont, and Jacques Felblinger. Joint reconstruction of multiple images and motion in MRI: Application to free-breathing myocardial t₂ quantification. *IEEE transactions on medical imaging*, 35(1):197–207, 2015.

[51] Freddy Odille, Sergio Uribe, Philip G. Batchelor, Claudia Prieto, Tobias Schaeffter, and David Atkinson. Model-based reconstruction for cardiac cine MRI without ECG or breath holding. *Magnetic Resonance in Medicine*, 63(5):1247–1257, 2010.

[52] Freddy Odille, Pierre-André Vuissoz, Pierre-Yves Marie, and Jacques Felblinger. Generalized reconstruction by inversion of coupled systems (GRICS) applied to free-breathing MRI.
*Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine*, 60(1):146-157, 2008.

[53] Bruno A Olshausen and David J Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. *Nature*, 381(6583):607-609, 1996.

[54] Ricardo Otazo, Emmanuel Candès, and Daniel K Sodickson. Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components. *Magnetic Resonance in Medicine*, pages 1125-1136, 2015.

[55] Ricardo Otazo, Daniel Kim, Leon Axel, and Daniel K Sodickson. Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI. *Magnetic Resonance in Medicine*, pages 767-776, 2010.

[56] Nora Ouzir, Adrian Basarab, Olivier Lairez, and Jean-Yves Tourneret. Robust optical flow estimation in cardiac ultrasound images using a sparse representation. *IEEE Transactions on Medical Imaging*, pages 741-752, 2018.

[57] Henrik Pedersen, Sebastian Kozerke, Steffen Ringgaard, Kay Nehrke, and Won Yong Kim. k-t PCA: temporally constrained k-t BLAST reconstruction using principal component analysis. *Magnetic Resonance in Medicine*, pages 706-716, 2009.

[58] Javier Sánchez Pérez, Enric Meinhardt-Llopis, and Gabriele Facciolo. TV-L1 optical flow estimation. *Image Processing On Line*, 2013:137-150, 2013.

[59] Christian Plathow, Sebastian Ley, Julia Zaporozhan, Max Schöbinger, Ekkehard Gruenig, Michael Puderbach, Monika Eichinger, Hans-Peter Meinzer, Ivan Zuna, and Hans-Ulrich Kauczor. Assessment of reproducibility and stability of different breath-hold maneuvers by dynamic MRI: comparison between healthy adults and patients with pulmonary hypertension. *European radiology*, 16(1):173-179, 2006.

[60] Sven Plein, William HT Smith, John P Ridgway, Andrea Kassner, David J Beacock, Timothy N Bloomer, and Mohan U Sivananthan. Qualitative and quantitative analysis of regional left ventricular wall dynamics using real-time magnetic resonance imaging: comparison with conventional breath-hold gradient echo acquisition in volunteers and patients. *Journal of Magnetic Resonance Imaging: An Official Journal of the International Society for Magnetic Resonance in Medicine*, 14(1):23-30, 2001.

[61] Claudia Prieto, Philip G. Batchelor, D.L.G. Hill, Joseph V. Hajnal, Marcelo Guarini, and Pablo Irarrazaval. Reconstruction of undersampled dynamic images by modeling the motion of object elements. *Magnetic Resonance in Medicine*, 57(5):939-949, 2007.

[62] Maithra Raghu, Chiyuan Zhang, Jon Kleinberg, and Samy Bengio. Transfusion: Understanding transfer learning for medical imaging. In *Advances in neural information processing systems*, pages 3347-3357, 2019.

[63] Christopher M Rank, Thorsten Heußer, Maria TA Buzan, Andreas Wetscherek, Martin T Freitag, Julien Dinkel, and Marc Kachelrieß. 4D respiratory motion-compensated image reconstruction of free-breathing radial MR data with very high undersampling. *Magnetic Resonance in Medicine*, 77(3):1170-1183, 2017.

[64] Stefan Roth and Michael J Black. On the spatial statistics of optical flow. *International Journal of Computer Vision*, 74(1):33-50, 2007.

[65] Javier Royuela-del Val, Lucilio Cordero-Grande, Federico Simmross-Wattenberg, Marcos Martín-Fernández, and Carlos Alberola-López. Nonrigid groupwise registration for motion estimation and compensation in compressed sensing reconstruction of breath-hold cardiac cine MRI. *Magnetic Resonance in Medicine*, 75(4):1525-1536, 2016.

[66] Javier Royuela-del Val, Lucilio Cordero-Grande, Federico Simmross-Wattenberg, Marcos Martín-Fernández, and Carlos Alberola-López. Jacobian weighted temporal total variation for motion compensated compressed sensing reconstruction of dynamic MRI. *Magnetic Resonance in Medicine*, 77(3):1208-1215, 2017.

[67] Todd S Sachs, Craig H Meyer, Pablo Irarrazabal, Bob S Hu, Dwight G Nishimura, and Albert Macovski. The diminishing variance algorithm for real-time reduction of motion artifacts in
MRI. Magn Reson Med, 34(3):412-422, 1995.

[68] Andres Saucedo, Stamatios Lefkimmiatis, Novena Rangwala, and Kyunghyun Sung. Improved computational efficiency of locally low rank MRI reconstruction using iterative random patch adjustments. *IEEE transactions on medical imaging*, pages 1209–1220, 2017.

[69] Matthias Schloegl, Martin Holler, Andreas Schwarzl, Kristian Bredies, and Rudolf Stollberger. Infimal convolution of total generalized variation functionals for dynamic MRI. *Magnetic Resonance in Medicine*, pages 142–155, 2017.

[70] Randy M Setser, Stefan E Fischer, and Christine H Lorenz. Quantification of left ventricular function with magnetic resonance images acquired in real time. *J Magn Reson Imaging*, 12(3):430–438, 2000.

[71] Xiaohui Shen and Ying Wu. Sparsity model for robust optical flow estimation at motion discontinuities. In *2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, pages 2456–2463. IEEE, 2010.

[72] A. Sotiras, C. Davatzikos, and N. Paragios. Deformable medical image registration: A survey. *IEEE Transactions on Medical Imaging*, 32(7):1153–1190, 2013.

[73] Deqing Sun, Stefan Roth, and Michael J Black. A quantitative analysis of current practices in optical flow estimation and the principles behind them. *International Journal of Computer Vision*, 106(2):115–137, 2014.

[74] Deqing Sun, Stefan Roth, JP Lewis, and Michael J Black. Learning optical flow. In *European Conference on Computer Vision*, pages 83–97. Springer, 2008.

[75] Deqing Sun, Xiaodong Yang, Ming-Yu Liu, and Jan Kautz. PWC-Net: CNNs for optical flow using pyramid, warping, and cost volume. In *IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pages 8934–8943, 2018.

[76] Radu Timofte and Luc Van Gool. Sparse flow: Sparse matching for small to large displacement optical flow. In *IEEE Winter Conference on Applications of Computer Vision (WACV)*, pages 1100–1106, 2015.

[77] J Trzasko, A Manduca, and E Borisch. Local versus global low-rank promotion in dynamic MRI series reconstruction. In *Proc. Int. Symp. Magn. Reson. Med*, page 4371, 2011.

[78] Paul Tseng. Convergence of a block coordinate descent method for nondifferentiable minimization. *Journal of optimization theory and applications*, 109(3):475–494, 2001.

[79] Muhammad Usman, David Atkinson, Freddy Odille, Christoph Kolbitsch, Ghislain Vaillant, Tobias Schaeffter, Philip G. Batchelor, and Claudia Prieto. Motion corrected compressed sensing for free-breathing dynamic cardiac MRI. *Magnetic Resonance in Medicine*, 70(2):504–516, 2013.

[80] Julia V Velikina and Alexey A Samsonov. Reconstruction of dynamic image series from undersampled MRI data using data-driven model consistency condition (MOCCO). *Magnetic Resonance in Medicine*, pages 1279–1290, 2015.

[81] Martin von Siebenthal, Gabor Szekely, Urs Gamper, Peter Boesiger, Antony Lomax, and Ph Cattin. 4D MR imaging of respiratory organ motion and its variability. *Physics in Medicine & Biology*, 52(6):1547, 2007.

[82] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli. Image quality assessment: from error visibility to structural similarity. *IEEE Transactions on Image Processing (TIP)*, pages 600–612, 2004.

[83] Daniel S Weller, Luonan Wang, John P Mugler III, and Craig H Meyer. Motion-compensated reconstruction of magnetic resonance images from undersampled data. *Magnetic resonance imaging*, 55:36–45, 2019.

[84] Manuel Werlberger, Thomas Pock, and Horst Bischof. Motion estimation with non-local total variation regularization. In *2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, pages 2464–2471. IEEE, 2010.

[85] Lukas Wissmann, Claudio Santelli, William P Segars, and Sebastian Kozerke. MRXCAT: Realistic numerical phantoms for cardiovascular magnetic resonance. *Journal of Cardiovascular Magnetic Resonance*, page 63, 2014.
[86] Jonas Wulff and Michael J Black. Efficient sparse-to-dense optical flow estimation using a learned basis and layers. In *Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition*, pages 120–130, 2015.

[87] Li Xu, Jiaya Jia, and Yasuyuki Matsushita. Motion detail preserving optical flow estimation. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 34(9):1744–1757, 2011.

[88] Maxim Zaitsev, Julian Maclaren, and Michael Herbst. Motion artifacts in MRI: a complex problem with many partial solutions. *Journal of Magnetic Resonance Imaging*, pages 887–901, 2015.

[89] Deheng Zhang, Jinxu Tao, Zhongfu Ye, Bensheng Qiu, and Jinzhang Xu. Deformation corrected blind compressed sensing (DC-BCS) a novel framework for dynamic MRI reconstruction. In *Proceedings of the 2018 2nd International Conference on Algorithms, Computing and Systems*, pages 160–164, 2018.

[90] Tao Zhang, John M Pauly, and Ives R Levesque. Accelerating parameter mapping with a locally low rank constraint. *Magnetic Resonance in Medicine*, pages 655–661, 2015.

[91] Bo Zhao, Justin P Haldar, Anthony G Christodoulou, and Zhi-Pei Liang. Image reconstruction from highly undersampled (k, t)-space data with joint partial separability and sparsity constraints. *IEEE Transactions on Medical Imaging*, pages 1809–1820, 2012.

[92] N. Zhao, D. O'Connor, A. Basarab, D. Ruan, and K. Sheng. Motion compensated dynamic MRI reconstruction with local affine optical flow estimation. *IEEE Transactions on Biomedical Engineering*, pages 1–1, 2019.