Complex extension of optical flow and its practical evaluation for undersampled dynamic MRI

Abstract.

Reconstructing high-quality images from undersampled dynamic MRI data is a challenging task and important for the success of this imaging modality. To remedy the naturally occurring artifacts due to measurement undersampling, one can incorporate a motion model into the reconstruction so that information can propagate across time frames. Current models for MRI imaging are using the optical flow equation. However, they are based on real-valued images. Here, we generalise the optical flow equation to complex-valued images and demonstrate, based on two real cardiac MRI datasets, that the new model is capable of improving image quality.

1. Introduction

Magnetic resonance imaging (MRI) is a powerful medical imaging modality that allows for noninvasive visualisations of the human body. In particular, for dynamic MRI, there is a variety of applications, e.g.  cardiac imaging [25, 2], dynamic contrast-enhanced magnetic resonance imaging [13], [9] and flow imaging [24]. However, the imaging times in MRI are inherently slow, leading to motion artifacts that degrade image quality. This issue can be addressed from a software perspective where one aims for acquisition speedups by decreasing the number of measurements and using mathematical modelling to improve image quality. In this article, we follow this approach and solve the ill-posed inverse problem $A_{t}\rho_{t}=y_{t}$ (the index $t$ denotes the $t$-th motion state of the considered subject; the MRI forward operator $A_{t}$ includes subsampling of measurements and hence may depend on $t$) of undersampled dynamic MRI by considering the variational problem

$$\min_{\rho,v}\sum_{t}\|A_{t}\rho_{t}-y_{t}\|^{2}+\mathcal{R}(\rho,v),$$

where $\mathcal{R}$ is a regularisation term that incorporates prior knowledge about the images into the reconstruction. More precisely, we aim to jointly reconstruct the temporal image sequence $\rho$ and the velocity field $v$ associated with the temporal evolution. To this end, we model the underlying motion using the optical flow equation [6]. This was already done for MRI and real-valued images [1]. Here, we generalise this idea by explicitly taking into account that images are complex-valued, i.e.  we formulate an optical flow equation for complex-valued images. This is important since omitting the phase may lead to reconstruction artifacts. This problem can be seen in fig. 1. The figure shows reconstructions from fully sampled MRI measurements. The images on the right show a nonnegative real-valued reconstruction where some structures, marked by red boxes, of the complex-valued reconstructions on the left are missing. The top row shows the magnitude of the reconstructions, while the bottom row depicts complex values: The hue represents the normalised phase and the brightness represents the magnitude. In particular, one can see that structures where the phase is different are missing in the real-valued images.

1.1. Related work

Dynamic inverse problems in general have been studied extensively. An overview of the topic can be found in [14, 23]. The particular topic of undersampled (dynamic) MRI has also been widely studied and many different approaches have been proposed to enhance reconstruction quality. Many works leveraged compressed sensing making use of the fact that MRI images can be sparsely represented in suitable domains [17, 16, 10]. In dynamic imaging, the temporal relation between motion states should be taken into account to improve the results. One approach is to decompose the (temporal) sequence of images into a low-rank part (background between time frames) and a sparse part (modelling dynamics between frames) [12, 18] . This way, the image sequence is represented suitably for the undersampled MRI problem. Another approach is taken in [26], where the motion between different motion states is estimated based on the frame-wise undersampled reconstruction of each state. Subsequently, this motion information is used to get a final reconstruction incorporating all motion states into a single reconstruction. Joint reconstruction and motion estimation approaches are also promising. This was investigated in [1, 29] where the underlying motion was modelled using (locally affine) optical flow. Besides classical image reconstruction, learning-based approaches have been considered in dynamic MRI. Those aim at learning spatiotemporal dependencies within the image sequence, and/or mimic the iterative nature of classical reconstruction algorithms [21, 20, 28].

Besides MRI, dynamic inverse problems arise in other fields of imaging. For instance, in dynamic PET imaging the continuity equation was used as the motion model to describe temporal changes between time points [22]. In tomography, an approach based on a template and its temporal deformation was taken [7].

When it comes to working with complex-valued images, two approaches are conceivable: First, one can parameterise the images using real and imaginary part. This is the standard approach that is basically always followed. However, in MRI imaging often the magnitude of the image is of most interest and one would like to pose prior assumptions on the magnitude. This is difficult to incorporate using real and imaginary part because calculating the magnitude is a nonlinear operation, hence the overall regularisation would involve a nonlinear operator. Therefore, one can work with a magnitude and phase parameterisation of the complex-valued images. This makes regularising the magnitude easier but comes at the cost that now the MRI forward operator is nonlinear. This was considered in [11, 27].

2. Methods

2.1. MRI modelling

Solving the inverse problem in MRI imaging in its simplest form leads to inverting the discrete Fourier transform. Given enough measurements (i.e.  enough coverage of Fourier space/k-space) this problem is well-posed by the stable invertibility of the Fourier transform. In practice, however, multiple receiver coils (here are the measurements detected) and undersampling lead to an ill-posed inverse problem. The MRI forward operator maps a complex-valued image $\rho_{t}$ (which in our case is a time slice of a temporal sequence of images) to some complex space $Y_{t}$, where $Y_{t}$ depends on the undersampling factor and the number of receiver coils. For instance, in a continuous setting one could choose $\rho\in L^{p}((0,T);BV(\Omega))$, hence $\rho_{t}\in BV(\Omega)$, or similar spaces with more regularity [6]. In a discrete setting (which will be considered in this article) we take $\rho\in X=\mathbb{C}^{N_{t}\times N_{x}\times N_{y}}$, where $N_{t}$, $N_{x}$ and $N_{y}$ denote the number of dimensions in the temporal, $x$ and $y$ direction, respectively. Therefore, it is $\rho_{t}\in X_{t}=\mathbb{C}^{N_{x}\times N_{y}}$ for each time slice. We next introduce the forward operator describing the measurements. Denoting by $N_{c}$ the number of receiver coils, we have $Y_{t}=\mathbb{C}^{N_{c}\times N_{k,t}}$, $t=1,\ldots,N_{t}$, as the measurement space. Here, $N_{k,t}$ denotes the number of measured k-space points at time $t$ after the undersampling has been applied. With the Fourier operator $F$, the coil sensitivity operators $C_{t,i}$ and the undersampling operators $S_{t}$ (see section 4.4 for details) we can find for the $i$-th component (i.e.  coil) of the MRI forward operator $A_{t}$ at time point $t$

(1) $$A_{t}\colon X_{t}\to Y_{t},\quad(A_{t}\rho_{t})_{i}=S_{t}FC_{t,i}\rho_{t}.$$

Hence, we model MRI measurements as

$$y_{t}=A_{t}\rho^{\dagger}_{t}+\varepsilon_{t}$$

where $\rho^{\dagger}_{t}\in X_{t}$ is the ground truth image that we want to reconstruct and $\varepsilon_{t}$ is some additive noise.

Figure 1. Comparison between real-valued and complex-valued reconstructions in MRI. Left: Complex-valued reconstruction. Right: Real-valued reconstruction. Top row: Magnitude representation. The bottom row shows the same images but as a complex-valued image: The hue represents the normalised phase and the brightness represents the magnitude. The red boxes point to some missing structures in the real-valued image.

2.2. Modelling motion for complex-valued images

In this work, we use the optical flow equation as a motion model to regularise the dynamic inverse problem. The optical flow equation is often used in computer vision and has proven valuable for 2D images [1, 5, 6]. Its defining property is constant image intensity along a trajectory with $\frac{\mathrm{d}}{\mathrm{d}t}\rho=v$ [14]. Thus, the underlying PDE reads

$$\partial_{t}\rho+\langle v,\nabla\rho\rangle_{\mathbb{R}}=0$$

and during the reconstruction, both the image $\rho$ as well as the velocity field $v$ have to be reconstructed. Note, that since velocity and density are multiplied by each other, using the optical flow equation makes the problem nonconvex. To be able to apply the optical flow equation in MRI imaging (where complex-valued images are reconstructed), we propose to generalise this equation by allowing both the image $\rho$ and the velocity $v$ to be complex-valued. The inner product is therefore replaced by the standard inner product for complex vector spaces. This leads us to the model

(2) $$0=\partial_{t}(\rho^{1}+i\rho^{2})+\langle v^{1}+iv^{2},\nabla(\rho^{1}+i\rho^{2})\rangle_{\mathbb{C}}$$

$$\Longleftrightarrow\begin{cases}0=\partial_{t}\rho^{1}+\langle v^{1},\nabla\rho^{1}\rangle_{\mathbb{R}}+\langle v^{2},\nabla\rho^{2}\rangle_{\mathbb{R}}\\ 0=\partial_{t}\rho^{2}+\langle v^{2},\nabla\rho^{1}\rangle_{\mathbb{R}}-\langle v^{1},\nabla\rho^{2}\rangle_{\mathbb{R}}.\end{cases}$$

for $\rho=\rho^{1}+i\rho^{2}\in X$ and $v=v^{1}+iv^{2}\in V^{2}$ (note, that the velocity is a 2D vector). Here, $V$ is a suitable space for the velocities such that the derivatives needed for the priors are well-defined. In the discrete setting, we may choose $V=X$. Associated with the above two equations we introduce the optical flow operator

(3) $$M(\rho,v)=\begin{pmatrix}\partial_{t}\rho^{1}+\langle v^{1},\nabla\rho^{1}\rangle_{\mathbb{R}}+\langle v^{2},\nabla\rho^{2}\rangle_{\mathbb{R}}\\ \partial_{t}\rho^{2}+\langle v^{2},\nabla\rho^{1}\rangle_{\mathbb{R}}-\langle v^{1},\nabla\rho^{2}\rangle_{\mathbb{R}}\end{pmatrix}.$$

For our discrete space $X=\mathbb{C}^{N_{t}\times N_{x}\times N_{y}}$ the temporal derivative $\partial_{t}\colon X\to X$ and the inner product $\langle\cdot,\cdot\rangle_{\mathbb{R}}\colon X\times X\to X$ are to be understood pointwise.

We use the same discrete versions of the differential operators as in [6]. This means we use forward differences in time and central differences in space with Neumann boundary conditions for the differential operators involved in the optical flow equation. For the differential operators of other parts of the optimization function we use forward differences with Neumann boundary conditions (see [6] for more details).

Remark (Choices on optical flow model).

There are other possible generalisations of the optical flow equation to complex-valued images. For instance, one could introduce separate velocities for real and imaginary part, hence decoupling both components. Also, a single velocity for both the real and imaginary part would be an option. Finally, only one of the two parts could be regularised by the optical flow equation and the other one evolves freely. However, these variants do not explicitly model interchange between real and imaginary part. To be able to model such an interchange, we used the above approach.

Remark (Alternative choice for motion model).

An alternative motion model would be to use the continuity equation in density and momentum variables instead of density and velocity variables (see, for example, [4]). This way, the motion model does not incorporate a nonconvexity into the optimisation problem. Unfortunately, finding meaningful regularisers for the momentum variables is more challenging than for the velocities. For instance, an object moving along the $x$ axis with constant velocity will have a zero gradient of the corresponding velocities. The momentum, however, is nonconstant as it also changes when the mass of the underlying object changes. Therefore, we decided to use the optical flow equation.

3. Optimisation and Implementation

Based on [1, 6] we minimise

(4) $$F(\rho,v)=\sum_{t}||A_{t}\rho_{t}-y_{t}||^{2}+\alpha_{1}\mathcal{R}_{1}(\rho)+\alpha_{2}\mathcal{R}_{2}(v)+\alpha_{3} \mathcal{R}_{3}(M(\rho,v))$$

with respect to $\rho=\rho^{1}+i\rho^{2},\ v=(v_{x}^{1},v_{y}^{1})+i(v_{x}^{2},v_{y}^{2})$. Here, $y_{t}$ is the measured MRI signal at time point $t$, and $A_{t}$ and $M$ are operators related to the MRI measurements and the optical flow equation, respectively (see (1) and (3)). The regularisers for the variables $\rho$ and $v$ are chosen to promote smoothness (the “level of smoothness” can be controlled via the $\varepsilon_{i}$ parameters)

$$\mathcal{R}_{1}(\rho)=\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{1})+\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{2})$$

$$\mathcal{R}_{2}(v)=\mathcal{H}_{\varepsilon_{2}}(\nabla v^{1}_{x})+\mathcal{H}_{\varepsilon_{2}}(\nabla v^{1}_{y})+\mathcal{H}_{\varepsilon_{2}}(\nabla v^{2}_{x})+\mathcal{H}_{\varepsilon_{2}}(\nabla v^{2}_{y}),$$

where $\nabla$ denotes the spatial gradient and $\mathcal{H}_{\varepsilon}$ is the Huber loss with parameter $\varepsilon$. Defining for $x\in\mathbb{R}^{2}$

$$h_{\varepsilon}(x)=\begin{cases}\frac{\|x\|^{2}_{2}}{2\varepsilon}\ \ &\text{if }\|x\|_{2}\leq\varepsilon,\\ \|x\|_{2}-\frac{\varepsilon}{2}\ \ &\text{else},\end{cases}$$

we have $\mathcal{H}_{\varepsilon}(x)=\sum_{i=1}^{N}h_{\varepsilon}(x_{i})$ for $x\in(\mathbb{R}^{2})^{N}$. Finally, we use a soft constraint

$$\mathcal{R}_{3}(M(\rho,v))=\mathcal{H}_{\varepsilon_{3}}(M(\rho,v))$$

since we cannot expect the optical flow equation to be satisfied perfectly. Compared to the classical $\mathrm{L}^{1}$ loss the Huber loss is differentiable, and the parameter $\varepsilon$ controls the extent of smoothness within the Huber loss, i.e. it controls whether the Huber loss is more similar to the norm $\|\cdot\|_{2}$ or the squared norm $\|\cdot\|_{2}^{2}$. Tuning this parameter allows us to adapt to the considered problem and find a better model.

The considered optimisation problem (4) is nonconvex due to the optical flow equation. To deal with the nonconvexity we follow the same strategy as in [1, 6] and apply a block coordinate descent (see algorithm 2). This means we iteratively optimise the minimisation function only for $\rho$ or $v$ while keeping the other variable fixed. Each resulting subproblem is solved using the FISTA algorithm [3] (applied to $f+g$ with $f$ smooth and $g=0$, see algorithm 1) with early stopping (iterations where stopped if the relative $\mathrm{L}^{2}$ error was below $\delta=10^{-5}$). Furthermore, motivated by the coarse-to-fine approach in image registration, we smoothed the density $\rho$ and the velocity $v$ after each sub-optimisation using a Gaussian smoothing with standard deviation $\sigma/{i_{\mathrm{outer}}}$. Here, $\sigma$ is again a tunable parameter, and $i_{\mathrm{outer}}\in\{1,\ldots,n_{\mathrm{outer}}\}$ denotes the current outer iteration, i.e. the smoothing decays with growing iterations. For the maximal number of iterations we used $n_{\mathrm{outer}}=200$, $n_{\rho}=1400$ and $n_{v}=3200$. The optimisation algorithm was implemented in Python using SIRF [19] for reading the MR data and applying the associated measurement operators, and CIL [15] to implement the actual optimisation algorithm. The code will be made available upon publication and is attached for the review process.

Algorithm 1 FISTA with early stopping for minimizing smooth function $f$

1:Inputs:

2:      smooth function $f$ with Lipschitz constant $L$, initialisation $x_{\mathrm{init}}$ Algorithm parameters $n$, $\delta$

3:Initialize:

4:      $x^{0}\leftarrow x_{\mathrm{init}}$, $\hat{x}^{1}\leftarrow x_{\mathrm{init}}$, $t^{1}\leftarrow 1$

5:for $j=1$ to $n$ do

6:     $x^{j}\leftarrow$ $\hat{x}^{j}-\frac{1}{L}\nabla f(\hat{x}^{j})$

7:     $t^{j+1}\leftarrow\frac{1+\sqrt{1+4(t^{j})^{2}}}{2}$

8:     $\hat{x}^{j+1}\leftarrow x^{j}+(\frac{t^{j}-1}{t^{j+1}})(x^{j}-x^{j-1})$

9:     if $\|x^{j}-x^{j-1}\|_{2}/\|x^{j-1}\|_{2}<\delta$ then

10:         break      

11:return $x^{j}$

Algorithm 2 Optimization of $F(\rho,v)$

1:Inputs:

2:      Measurement $y$ Model parameters $\alpha_{1},\alpha_{2},\alpha_{3},\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}$ Algorithm parameters $\sigma,n_{\mathrm{outer}},n_{\rho},n_{v}$, $\delta$

3:Initialize:

4:      $v^{1}\leftarrow 0$, $\hat{v}^{1}\leftarrow 0$, $\rho^{1}\leftarrow 0$, $\hat{\rho}^{1}\leftarrow 0$

5:for $i=1$ to $n_{\mathrm{outer}}$ do

6:     $\rho^{i+1}\leftarrow\text{FISTA}(\rho\mapsto F(\rho,\hat{v}^{i}),\hat{\rho}^{i})$ with $n_{\rho}$ iterations

7:     $\hat{\rho}^{i+1}\leftarrow\mathrm{smooth}_{\sigma/i}(\rho^{i+1})$

8:     $v^{i+1}\leftarrow\text{FISTA}(v\mapsto F(\hat{\rho}^{i+1},v),\hat{v}^{i})$ with $n_{v}$ iterations

9:     $\hat{v}^{i+1}\leftarrow\mathrm{smooth}_{\sigma/i}(v^{i+1})$

10:     if $\frac{1}{2}\|v^{i+1}-v^{i}\|_{2}/\|v^{i}\|_{2}+\frac{1}{2}\|\rho^{i+1}-\rho^{i}\|_{2}/\|\rho^{i}\|_{2}<\delta$ then

11:         break      

12:return $\rho^{i+1},v^{i+1}$

3.1. Models for comparison

We compare our proposed method (4) (we call this model OF  and label corresponding reconstructions OF, too) to three slightly different variants:

• Frame-wise model (FW): For each time point a separate reconstruction is computed and no correlation between time points is assumed. This corresponds to optimising (4) for the density $\rho$ only and setting $\alpha_{2}$ and $\alpha_{3}$ to zero. Thus, we minimise

  $$F(\rho)=\sum_{t}||A_{t}\rho_{t}-y_{t}||^{2}+\alpha_{1}\mathcal{R}_{1}(\rho)$$

  with respect to $\rho=\rho^{1}+i\rho^{2}$ applying algorithm 1 with $n_{\rho}$ iterations, where

  $$\mathcal{R}_{1}(\rho)=\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{1})+\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{2})$$

  as above.

• Time derivative model (DT): The term containing the optical flow equation is replaced by a term containing time derivatives only. This means we fix the velocities to zero and optimise (4) for the density $\rho$ only. Thus, we minimise

  $$F(\rho)=\sum_{t}||A_{t}\rho_{t}-y_{t}||^{2}+\alpha_{1}\mathcal{R}_{1}(\rho)+\alpha_{3}\mathcal{R}_{3}(M(\rho,0))$$

  with respect to $\rho^{1}+i\rho^{2}$ applying algorithm 1 with $n_{\rho}$ iterations, where

  $$\mathcal{R}_{1}(\rho)=\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{1})+\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{2}).$$

  and

  $$\mathcal{R}_{3}(M(\rho,0))=\mathcal{H}_{\varepsilon_{3}}([\partial_{t}\rho^{1},\partial_{t}\rho^{2}]).$$

  as above promote spatial and temporal regularity, respectively.

• Optical flow using ground truth velocities (Cheat-OF): With this model we want to investigate how our method would perform if we had access to the velocities underlying the motion, i.e. this model gives us the best possible reconstruction within our framework. To this end, we use our model (4) with fully sampled data to determine ground truth velocities $v_{\mathrm{GT}}$. Then, we optimise (4) for the density $\rho$ only while fixing the velocity to be $v_{\mathrm{GT}}$. Thus, we minimise

  $$F(\rho)=\sum_{t}||A_{t}\rho_{t}-y_{t}||^{2}+\alpha_{1}\mathcal{R}_{1}(\rho)+\alpha_{3}\mathcal{R}_{3}(M(\rho,v_{\mathrm{GT}}))$$

  with respect to $\rho=\rho^{1}+i\rho^{2}$ applying algorithm 1 with $n_{\rho}$ iterations, where

  $$\mathcal{R}_{1}(\rho)=\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{1})+\mathcal{H}_{\varepsilon_{1}}(\nabla\rho^{2})$$

  and

  $$\mathcal{R}_{3}(M(\rho,v_{\mathrm{GT}}))=\mathcal{H}_{\varepsilon_{3}}(M(\rho,v_{\mathrm{GT}})$$

  as above.

3.2. Metrics

We consider two metrics to compare the quality of two images: PSNR (which is also used as the metric for the hyperparameter optimisation) and SSIM (which is the secondary alternative metric to compare image quality). To compute them, we use the respective functions provided by scikit-image where the ”data_range“ parameter is obtained from the ground truth image. Also, we compute the metrics for each time point and compute the final metric as the average over all time points.

3.3. Hyperparameter optimisation

Hyperparameter optimisation was carried out using scikit-optimize to find good sets of parameters for each model. We chose PSNR (computed on the dynamic mask, see section 4.2) as the metric to be minimised by the Bayesian optimisation. For the Bayesian optimisation we used the Optimizer class with a Gaussian process as the base estimator and expected improvement as the acquisition function. For the model in (4) and its optimisation we used the following hyperparameters. The regularisation parameters $\alpha_{1},\alpha_{2},\alpha_{3}$, the smoothness parameters $\varepsilon_{1},\varepsilon_{2},\varepsilon_{3}$ and the algorithm parameter $\sigma$ which gives the initial amount of smoothing of the density and velocity after each outer optimisation step. The number of tested parameters depends on the model. The more model parameters exist, the more combinations of parameters were considered within the Bayesian optimisation. More details can be found in the code.

4. Data

We work with two datasets taken from the OCMR dataset [8] containing cardiac cine series (long-axis view). The data consists of multi-coil k-space data acquired with Cartesian sampling on a Siemens MAGNETOM Avanto scanner (1.5T) (see [8] for more details). We use fully sampled data and apply the subsampling manually. A fully sampled reconstruction of both datasets can be found in fig. 2. Based on subject 1, we additionally consider a simulated dataset that resembles a real-world MRI image but satisfies the optical flow equation perfectly. This way, we can assess the model capacities and our optimisation procedure in an idealised setting.

4.1. Ground truth images

To assess the reconstruction quality, we need a suitable “ground truth” image to compare the results to. One possibility would be to use the fully sampled measurements to get reconstructions for each time frame separately using the FW  model. In doing so, however, one gets time trajectories that are slightly irregular and do not resemble realistic motion (see fig. 2 second column on the right) due to, e.g., in-frame motion. Therefore, we use a “dynamic ground truth” that results from minimising (4) with fully sampled data and high dynamic regularisation parameters. This results in reconstructions that are more regular in time (see fig. 2 first column on the right). This approach also has the advantage that we have “ground truth” velocities available that describe the dynamics within the ground truth image sequence. These can be used to assess the quality of our reconstructed velocities.

In the discussion section, we will briefly comment on the results being similar when using the frame-wise “ground-truth”.

4.2. Mask focusing on dynamic parts in images

We are primarily interested in the dynamic parts of the image and want to improve image quality in these regions particularly. Since large parts of the images are basically constant in time, computing loss metrics between ground truth images and reconstructions over the whole image would consequently result in focusing on the wrong parts. Therefore, we computed masks focusing on the dynamic parts of the images (see fig. 2 right column), and computed the image metrics on these masks only.

Figure 2. Each row shows a reconstruction for a different data set. The first column shows a dynamic ground truth (obtained as a reconstruction using (4) and fully-sampled data). The second column shows the same for a frame-wise ground truth (obtained from reconstructing each time frame separately (FW model) using fully-sampled data). The last column shows an overlay of the reconstructions with the dynamic masks. Within each column, the left side shows one time slice of the reconstructions, the right side shows time-space images along the red lines.

4.3. Simulated data

We test our proposed model on simulated data. To this end, we used a dataset from [8] (subject 1 in fig. 2) and solved (4) for fully sampled data. Then, we smoothed the obtained velocities in space and time and used the resulting velocities to simulate a sequence of images starting with the image at the first time point from solving (4). The thus obtained images and velocities are therefore guaranteed to satisfy the optical flow equation and the velocities are smooth enough to be reconstructable. Based on these simulated ground truth images we generated measurements by applying the MRI forward operator $A_{t}$. Finally, Gaussian noise (zero mean) was added where the standard deviation was chosen such that unregularised reconstructions yielded qualitatively similar images to unregularised reconstructions based on the original measured data.

4.4. Undersampling

We next describe the used undersampling. The data we are working with has been acquired with a standard Cartesian sampling scheme, and we consider four times undersampling in the phase encoding direction with full sampling of the rows in the readout direction. More precisely, we use the same $x=15$ percent central rows in each time frame and randomly select $25-x$ percent of the outer rows (same number of rows on each side of the central rows) in such a way that the selected outer rows between consecutive time frames do not overlap.

5. Experiments and results

5.1. Simulated data

The results for the simulated data are summarised in table 1 and can visually be inspected in fig. 3. The top row in that figure shows the magnitude of the simulated images. In the second row, the same images are represented as complex-valued images where the hue represents the normalised phase and the brightness represents the magnitude. From this row of images, we see that the phase varies spatially, but its temporal changes are small.

The other rows depict the difference between the reconstructions using the respective model and the simulated images (magnitude images are compared). We can clearly see that DT  improves over the frame-wise approach FW, and that OF  improves over DT. The improved reconstruction qualities can also be seen in PSNR values in table 1 (first column), where OF  improves over FW  and DT  by 10dB and 4dB, respectively. Moreover, we see, that Cheat-OF , for which ground truth velocities are available in the reconstruction, is basically not improving over OF. Similar observations can be made by comparing the SSIM values in table 1.

This good reconstruction result of the OF  model is backed by the reconstructed velocities depicted in fig. 4 (bottom row). We see that the general form of the velocities has been recovered compared to the ground truth velocities (top row), hence a good reconstruction quality is expected.

Figure 3. Results for the simulated data. The top row shows the evolution (magnitude) of the simulated ground truth images for the time points $t_{0}$, $t_{2}$, $t_{4}$, $t_{6}$ (out of $8$ time points). The second row shows the same images but as a complex-valued image: The hue represents the normalised phase and the brightness represents the magnitude. The other rows depict the difference between the magnitudes of the reconstructions and the magnitude of the simulated images on the dynamic mask. Rows 3-5: FW, DT  and OF.

Figure 4. Reconstructed velocities of one time frame for the simulated data using our proposed complex-valued optical flow method. First row: Ground truth velocities in $x$ and $y$ direction. Second row: Reconstructed velocities in $x$ and $y$ direction

5.2. Subject 1

The summarised results of subject 1 can be found in table 1 and fig. 5. We see that our proposed method OF  performs better than FW  or DT. The comparison between OF  and the frame-wise approach yields PSNR improvements (a similar pattern can be observed in the SSIM values) of almost 3dB, showing that the complex optical flow equation helps the reconstruction by describing the underlying motion. The gains of OF  over DT  are not very pronounced. However, they improve to over 1dB when looking at finer structures in more detail as can be seen in fig. 5 in the third row. Further, the results in table 1 and fig. 5 also show that the method Cheat-OF  including ground truth velocities outperforms all the other methods by a lot (which is no surprise since this method basically has access to the fully sampled measurements via the ground truth velocities obtained from a fully sampled reconstruction using (4)). This improvement can be explained by comparing the reconstructed to the “ground truth” velocities in fig. 7. We see that the reconstruction quality is decent for the real part, whereas the reconstructed imaginary part differs significantly from the “ground truth” velocities’ imaginary part. Overall, the velocities are not as well reconstructed as for the simulated data, which explains the improved results using Cheat-OF.

5.3. Subject 2

The summarised results of subject 2 can be found in table 1 and fig. 6. We see that OF  improves PSNR over FW  and DT  by 2dB and 1.4dB, respectively. Therefore, the proposed optical flow model helps with reconstructing dynamic MRI images. Again, a similar conclusion can be drawn from the SSIM values. Looking at the reconstructed phase (fig. 6 rows 3 and 4), we see that all methods reconstruct the phase well on an overall level (row 3). However, the magnified area in row 4 shows differences in the reconstructed phases in the detailed structures. Here, OF  gives a better result over FW  and DT. This improved phase reconstruction is associated with a better reconstruction of details, as can be seen in the second row. Again, Cheat-OF, as expected, improves over all other methods.

		simulation	subject 1 GT(dyn)	subject 2 GT(dyn)	subject 1 GT(FW)	subject 2 GT(FW)
	FW	$34.1\pm 1.2$	$27.3\pm 0.7$	$25.6\pm 0.6$	$27.7\pm 0.8$	$26.4\pm 0.9$
	DT	$40.1\pm 0.9$	$29.6\pm 0.3$	$26.2\pm 0.3$	$29.8\pm 0.4$	$27.5\pm 0.5$
	OF	$44.1\pm 0.6$	$30.2\pm 0.5$	$27.6\pm 0.4$	$30.2\pm 0.5$	$28.0\pm 0.5$
PSNR	Cheat-OF	$44.2\pm 0.3$	$34.4\pm 0.7$	$32.0\pm 0.7$	$33.1\pm 0.4$	$28.7\pm 1.0$
	FW	$0.943\pm 0.010$	$0.828\pm 0.014$	$0.812\pm 0.015$	$0.847\pm 0.014$	$0.852\pm 0.018$
	DT	$0.978\pm 0.004$	$0.869\pm 0.009$	$0.819\pm 0.012$	$0.880\pm 0.009$	$0.869\pm 0.011$
	OF	$0.991\pm 0.001$	$0.892\pm 0.010$	$0.873\pm 0.012$	$0.888\pm 0.011$	$0.880\pm 0.011$
SSIM	Cheat-OF	$0.992\pm 0.001$	$0.959\pm 0.009$	$0.942\pm 0.009$	$0.934\pm 0.004$	$0.872\pm 0.018$

Table 1. Summary of the results. The columns show the metrics for the respective data sets (with dynamic ground truth GT(dyn) and frame-wise ground truth GT(FW)). The displayed values are $(\text{average})\pm(\text{standard deviation})$ taken over the considered time points. The metric for each time point is computed on a dynamic mask (see section 4.2, fig. 2).

Figure 5. Reconstructions for subject 1. Left to right columns: Ground truth image, FW, DT, OF and Cheat-OF. First row shows magnitude of reconstructions. Second row: Difference of reconstruction to ground truth on the mask. Third row: Magnified area and corresponding metrics on magnification.

Figure 6. Reconstructions for subject 2. Left to right columns: Ground truth image, FW, DT, OF and Cheat-OF. First and second row show magnitude of reconstructions with respective metrics. Third and fourth row show the phase of the reconstructions. Second and last row are the indicated magnified areas.

Figure 7. Comparison between ground truth velocities (top row) and reconstructed velocities (bottom row) of subject 1 for a single time frame.

6. Discussion

Overall, our experiments demonstrated that explicitly considering the complex-valued optical flow equation (2) is useful for modelling motion between complex-valued images. In particular, our simulated experiment showed large reconstruction improvements. This, of course, is to be expected as the underlying ground truth motion perfectly satisfies the optical flow equation and the corresponding velocities follow the posed prior assumptions. Still, this simplified example shows that the optimisation procedure is suitable for the problem, and meaningful results are obtained. Moreover, we have mentioned already in section 5.1 that Cheat-OF  is basically not improving over OF  in the simulation, suggesting that OF  is able to produce nearly optimal reconstructions in case the underlying images satisfy the considered prior assumptions.

Comparing the ground truth velocities of the simulation (fig. 4) and of the measured subject 1 (fig. 7) we see that they are of different quality. In particular, the simulated velocity appears to be much smoother in space, hence they are more in accordance with our prior modelling. On the other hand, the ground truth velocities of subject 1 are much rougher in space which is not reflected in the reconstructed velocity. This of course influences reconstruction quality and partly explains the differences in reconstruction quality compared to the simulation. In future work, it will be interesting to investigate different priors that better resemble the underlying velocities’ structure. For instance, one could use Haar wavelets, include regularisation in time, or follow a learning-based approach that can derive very specialised and problem-specific priors.

For our “ground truth” images we argued in section 4.2 that using “dynamic ground truth” images is more appropriate in our setting. However, we repeated the procedures for subjects 1 and 2 with “frame-wise ground truth” images (for the simulation both versions of ground truth images coincide) with similar results, where PSNR improvements are slightly less pronounced (2.5dB over FW  for subject 1 and 1.6dB over FW  for subject 2, see table 1 right two columns). This suggests that the overall results do not depend on the choice of “ground truth” images.

7. Conclusions

Dynamic MRI (like cardiac MRI) is an important clinical imaging modality for which measurement undersampling is vital to get reasonably fast acquisitions. To remedy the image reconstruction artifacts naturally appearing due to this undersampling, motion models such as the optical flow equation can be used to propagate information across time frames. Current models have followed this approach for real-valued images. In this work, we have generalised the procedure to complex-valued images by proposing a complex version of the optical flow equation that can account for phase changes between time frames. We tested the new model on two real cardiac datasets and saw reconstruction improvements of the complex optical flow model over the comparison methods.

We used central differences for spatial derivatives in the optical flow equation. In future work, one could analyse which discrete differentiation approach works best in the context of the complex optical flow equation. Related is the question of whether a CFL condition changes when transitioning from the real to the complex optical flow equation.

We also saw that including “ground truth” velocities into the model further improved reconstruction quality implying that the currently chosen regularisation for the velocities is not optimal. In future work, it will be interesting to further investigate how to best choose regularisers for the velocities to achieve better performance (such an optimal choice is not even clear for the standard real-valued optical flow motion model).

Overall, we observed that the complex optical flow model is well-suited for cardiac MRI data. Of course, other options are conceivable. For instance, the continuity equation is a typical way of describing motion (especially in three dimensions). Also, a learning-based approach could be vital to adapt the model to the underlying problem so that problem-specific features can be taken into account. It will be interesting to compare different approaches with each other.