Deep variational network for rapid 4D flow MRI reconstruction

by   Valery Vishnevskiy, et al.

Phase-contrast magnetic resonance imaging (MRI) provides time-resolved quantification of blood flow dynamics that can aid clinical diagnosis. Long in vivo scan times due to repeated three-dimensional (3D) volume sampling over cardiac phases and breathing cycles necessitate accelerated imaging techniques that leverage data correlations. Standard compressed sensing reconstruction methods require tuning of hyperparameters and are computationally expensive, which diminishes the potential reduction of examination times. We propose an efficient model-based deep neural reconstruction network and evaluate its performance on clinical aortic flow data. The network is shown to reconstruct undersampled 4D flow MRI data in under a minute on standard consumer hardware. Remarkably, the relatively low amounts of tunable parameters allowed the network to be trained on images from 11 reference scans while generalizing well to retrospective and prospective undersampled data for various acceleration factors and anatomies.



page 2

page 4

page 5

page 6


Accurate quantification of blood flow wall shear stress using simulation-based imaging: a synthetic, comparative study

Simulation-based imaging (SBI) is a blood flow imaging technique that op...

End-to-End Variational Networks for Accelerated MRI Reconstruction

The slow acquisition speed of magnetic resonance imaging (MRI) has led t...

Memory-efficient Learning for High-Dimensional MRI Reconstruction

Deep learning (DL) based unrolled reconstructions have shown state-of-th...

Scale-Equivariant Unrolled Neural Networks for Data-Efficient Accelerated MRI Reconstruction

Unrolled neural networks have enabled state-of-the-art reconstruction pe...

Deep J-Sense: Accelerated MRI Reconstruction via Unrolled Alternating Optimization

Accelerated multi-coil magnetic resonance imaging reconstruction has see...

Compressive MRI quantification using convex spatiotemporal priors and deep auto-encoders

We propose a dictionary-matching-free pipeline for multi-parametric quan...

Free-breathing Cardiovascular MRI Using a Plug-and-Play Method with Learned Denoiser

Cardiac magnetic resonance imaging (CMR) is a noninvasive imaging modali...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.


Compressed Sensing 4D Flow Reconstruction

PC MRI encodes flow velocity at spatial location during cardiac phase () according to the following equation:


where is the velocity corresponding to a phase of ,

are the encoded velocity vector components. The four-point velocity encoding matrix is given as


Therefore, flow velocity can be calculated from the phase difference of reconstructed PC images .

Let be a discretized image on a grid corresponding to a cardiac phase and velocity encoding . Assuming Cartesian sampling on a regular

grid, the Fourier transform

, and coil sensitivity maps define the spatial encoding operator :


Considering a single velocity-encoded image sequence, let and be stacked column-vectors of signals and zero-filled k-space samples respectively, while defines the undersampling mask. Iterative image reconstruction methods seek for a maximum a posteriori (MAP) solution defined by the following optimization problem:


where the regularization term enforces prior assumptions about image regularities. Herein we consider the local low-rank (LLR) regularization 14 to leverage image correlations among cardiac phases:


where is the corresponding patch extraction operator, yielding patches, and is the nuclear norm. For LLR regularization the optimization problem (7) is convex and can be efficiently solved using operator splitting techniques such as the fast iterative shrinkage-thresholding algorithm (FISTA) 20.

FlowVN Training

We employ a =10 layer VN and perform 5 iterations of the ADAM algorithm (learning rate , , , batch size of 3) for training, during which we continually adjust with being the iteration number. On every layer, each 3D filter bank contains = filters of size = voxels. Activation functions are parametrized by = control knots with spacing =:


with gradients provided by the following formulas:


The acquired zero-filled k-space with undersampling mask was normalized by .

To enable backpropagation to be carried out with limited GPU memory, we employ spatio-temporal equivariance of the convolution and exploit the fact that k-space is fully sampled in the readout dimension

for Cartesian acquisitions. Therefore, to draw a training sample, we perform random cropping of width and in dimensions and respectively and simulate Fourier encoding in dimensions as illustrated in Fig. 2(a). The network was implemented using the Tensorflow framework 32. Fully-sampled and partial Fourier acquisition data from 11 healthy volunteers was used during training.

In Vivo Data Acquisition

As illustrated in Fig. 1, we used 11 subject for network training, and 7 healthy subjects and 1 patient for evaluation. All in-vivo work was performed upon written informed consent of the subjects and according to local ethics regulations.

Training datasets comprised 4D flow data measured in the aorta of 11 healthy subjects, 9 of them fully sampled and 2 acquired with partial Fourier 38 (factor 0.750.75).

For evaluation, data in the ascending aorta of 7 healthy subjects were acquired on a 3T Philips Ingenia system (Philips Healthcare, Best, the Netherlands) using a Cartesian 4-point referenced phase-contrast gradient-echo sequence with an encoding velocity , a spatial resolution of , , , cardiac phases and flip angle = 8. Exams for each of the 7 healthy subjects comprised a standard navigator-gated 2-fold accelerated parallel imaging 5 exam for reference, and a CS acquisition with an acceleration factor of =, using Cartesian pseudo-radial golden angle sampling pattern 39 and data driven-respiratory motion detection, as in 11. Only data in expiration were kept for reconstruction as shown in Fig. 1.

To evaluate reconstruction accuracy on pathological anatomy, 4D flow data was acquired in a single patient with dilation of the ascending aorta, combined aortic stenosis and regurgitation due to a bicuspid aortic valve on a 3T Philips Ingenia system (Philips Healthcare, Best, the Netherlands) using a navigator-gated 2-fold accelerated parallel imaging 5 scan.

A receiver coil with 28 channels was used for acquisition which were reduced to 5 channels using coil compression 40. Coil sensitivity maps were estimated with ESPIRiT 41. Concomitant field correction was applied to the signal phase according to Bernstein et al. 42 and eddy currents were corrected for with a third-order polynomial model fitted to stationary tissue 43, 44.


We compared the proposed FlowVN to the state-of-the-art compressed sensing LLR-regularized (8) reconstruction 14 and the variational network by Hammernik et al. 15 which we refer to as HamVN. The LLR implementation from the Berkeley advanced reconstruction toolbox (BART) 45 was used with patch size = and a maximum number of optimization iterations of 80. The optimal value of regularization parameter = was chosen via grid search to minimize the reconstructed flow field residual averaged over the manually segmented aorta on the retrospectively 12 undersampled acquisition. Since the original VN 15 was proposed for magnitude reconstruction of 2D and 3D data, we introduced the following modifications for the presented 4D flow evaluation: (i) 3D filters were grouped into 4 banks as in FlowVN (see Supplemental Algorithm 1), (ii) -norm of reconstruction was optimized (i.e. equation (3) with + ). We refer to this architecture as “HamVN”, the number of network layers, filters and control knots were the same as in FlowVN.

Retrospective Study

For simulated retrospective undersampling experiments, we used 2 PI data and simulate a pseudo-radial Golden angle sampling pattern 39 with acceleration factors of 6 to 22.

For each undersampling factor we evaluated the nRMSE of the image magnitude, the relative error (RelErr) of velocity magnitudes inside the aorta and the angular error (AngErr) of the estimated velocity vectors:


Additionally, we report the structural similarity index (SSIM) 46 with = on the reconstructed magnitude images.

Prospective Study

Using manual aorta segmentations we compute flow over cross sections of the aorta by integrating velocity components projected onto the cross section normal. The peak flow is then defined as the maximal flow over cardiac phases for a given cross section. Moreover, we calculate peak through-plane velocity defined as maximum velocity projection across cross sections of the aorta over cardiac phases.

To quantify agreement with the reference 2 PI reconstruction, we performed Bland-Altman analysis 47 of peak flow and peak through-plane velocities.

Data and code availability

The code for the network training and inference used in this study and network weights are avaliable on CodeOcean together with volunteer data: 48. The code for analysis is available on CodeOcean: 49.


  • 1 Markl, M., Frydrychowicz, A., Kozerke, S., Hope, M. & Wieben, O. 4D flow MRI. J. Mag. Res. Imag. 36, 1015–1036 (2012).
  • 2 Feinberg, D., Hale, J., Watts, J., Kaufman, L. & Mark, A. Halving MR imaging time by conjugation: demonstration at 3.5 kG. Radiol. 161, 527–531 (1986).
  • 3 Szarf, G. et al. Zero filled partial fourier phase contrast MR imaging: in vitro and in vivo assessment. J. Mag. Res. Imag. 23, 42–49 (2006).
  • 4 Walheim, J., Gotschy, A. & Kozerke, S. On the limitations of partial fourier acquisition in phase-contrast MRI of turbulent kinetic energy. Mag. Res. in Med. 81, 514–523 (2019).
  • 5 Pruessmann, K., Weiger, M., Scheidegger, M. & Boesiger, P. SENSE: sensitivity encoding for fast MRI. Mag. Res. in Med. 42, 952–962 (1999).
  • 6 Wiesinger, F., Boesiger, P. & Pruessmann, K. P. Electrodynamics and ultimate SNR in parallel MR imaging. Mag. Res. in Med. 52, 376–390 (2004).
  • 7 Lustig, M., Donoho, D. & Pauly, J. M. Sparse MRI: The application of compressed sensing for rapid MR imaging. Mag. Res. in Med. 58, 1182–1195 (2007).
  • 8 Kim, D. et al. Accelerated phase-contrast cine MRI using k-t SPARSE-SENSE. Mag. Res. in Med. 67, 1054–1064 (2012).
  • 9 Valvano, G. et al. Accelerating 4D flow MRI by exploiting low-rank matrix structure and hadamard sparsity. Mag. Res. in Med. 78, 1330–1341 (2017).
  • 10 Bollache, E. et al. k-t accelerated aortic 4D flow MRI in under two minutes: feasibility and impact of resolution, k-space sampling patterns, and respiratory navigator gating on hemodynamic measurements. Mag. Res. in Med. 79, 195–207 (2018).
  • 11 Walheim, J., Dillinger, H. & Kozerke, S. Multipoint 5D Flow Cardiovascular Magnetic Resonance - Accelerated Cardiac- and Respiratory-Motion Resolved Mapping of Mean and Turbulent Velocities. J. Cardiovasc. Magn. Res. (2019).
  • 12 Ma, L. E. et al. Aortic 4D flow MRI in 2 minutes using compressed sensing, respiratory controlled adaptive k-space reordering, and inline reconstruction. Mag. Res. in Med. (2019).
  • 13 Rich, A. et al. A bayesian approach for 4D flow imaging of aortic valve in a single breath-hold. Mag. Res. in Med. 81, 811–824 (2019).
  • 14 Zhang, T., Pauly, J. M. & Levesque, I. R. Accelerating parameter mapping with a locally low rank constraint. Mag. Res. in Med. 73, 655–661 (2015).
  • 15 Hammernik, K. et al. Learning a variational network for reconstruction of accelerated MRI data. Mag. Res. in Med. 79, 3055–3071 (2018).
  • 16 Mardani, M. et al. Deep generative adversarial neural networks for compressive sensing MRI. IEEE Trans. Med. Imag. 38, 167–179 (2018).
  • 17 Zhu, B., Liu, J. Z., Cauley, S. F., Rosen, B. R. & Rosen, M. S. Image reconstruction by domain-transform manifold learning. Nature 555, 487 (2018).
  • 18 Schlemper, J., Caballero, J., Hajnal, J. V., Price, A. N. & Rueckert, D.

    A deep cascade of convolutional neural networks for dynamic mr image reconstruction.

    IEEE Trans. Med. Imag. 37, 491–503 (2017).
  • 19 Maier, A. K. et al. Learning with known operators reduces maximum error bounds. Nat. Mach. Intel. 1, 373–380 (2019).
  • 20 Beck, A. & Teboulle, M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imag. Sc. 2, 183–202 (2009).
  • 21 Antun, V., Renna, F., Poon, C., Adcock, B. & Hansen, A. C. On instabilities of deep learning in image reconstruction-does AI come at a cost? arXiv preprint arXiv:1902.05300 (2019).
  • 22 Yang, G. et al. DAGAN: Deep de-aliasing generative adversarial networks for fast compressed sensing MRI reconstruction. IEEE Trans. Med. Imag. 37, 1310–1321 (2017).
  • 23 Quan, T., Nguyen-Duc, T. & Jeong, W. Compressed sensing MRI reconstruction using a generative adversarial network with a cyclic loss. IEEE Trans. Med. Imag. 37, 1488–1497 (2018).
  • 24 Narnhofer, D., Hammernik, K., Knoll, F. & Pock, T. Inverse GANs for accelerated MRI reconstruction. In Wav. and Spars., vol. 11138, 111381A (2019).
  • 25 Zhang, S., Block, K. & Frahm, J. Magnetic resonance imaging in real time: advances using radial FLASH. J. Mag. Res. Imag. 31, 101–109 (2010).
  • 26 Landweber, L. An iteration formula for Fredholm integral equations of the first kind. Am. J. of Math. 73, 615–624 (1951).
  • 27 Liu, Y. & Lew, M. S. Learning relaxed deep supervision for better edge detection. In CVPR, 231–240 (2016).
  • 28 Ravishankar, S. & Bresler, Y. MR image reconstruction from highly undersampled k-space data by dictionary learning. IEEE Trans. Med. Imag. 30, 1028–1041 (2010).
  • 29 Caballero, J., Price, A. N., Rueckert, D. & Hajnal, J. V. Dictionary learning and time sparsity for dynamic MR data reconstruction. IEEE Trans. Med. Imag. 33, 979–994 (2014).
  • 30 Lee, D., Yoo, J. & Ye, J. C. Deep artifact learning for compressed sensing and parallel MRI. arXiv preprint arXiv:1703.01120 (2017).
  • 31 LeCun, Y., Touresky, D., Hinton, G. & Sejnowski, T. A theoretical framework for back-propagation. In Proc. Connect. Mod., vol. 1, 21–28 (1988).
  • 32 Abadi, M. et al.

    Tensorflow: A system for large-scale machine learning.

    In Proc. USENIX, 265–283 (2016).
  • 33 Domke, J. Generic methods for optimization-based modeling. In Art. Int. and Stat., 318–326 (2012).
  • 34 Sun, J., Li, H., Xu, Z. et al. Deep ADMM-Net for compressive sensing MRI. In NIPS, 10–18 (2016).
  • 35 Jin, K. H., McCann, M. T., Froustey, E. & Unser, M. Deep convolutional neural network for inverse problems in imaging. IEEE Trans. on Imag. Proc. 26, 4509–4522 (2017).
  • 36 Vishnevskiy, V., Sanabria, S. J. & Goksel, O. Image reconstruction via variational network for real-time hand-held sound-speed imaging. In Intl. W. on Mach. Learn. for Med. Imag. Recon., 120–128 (2018).
  • 37 Vishnevskiy, V., Rau, R. & Goksel, O. Deep variational networks with exponential weighting for learning computed tomography. In MICCAI, 310–318 (2019).
  • 38 Cuppen, J. & van Est, A. Reducing MR imaging time by one-sided reconstruction. Mag. Res. Imag. 5, 526–527 (1987).
  • 39 Winkelmann, S., Schaeffter, T., Koehler, T., Eggers, H. & Doessel, O. An optimal radial profile order based on the Golden Ratio for time-resolved MRI. IEEE Trans. Med. Imag. 26, 68–76 (2006).
  • 40 Zhang, T., Pauly, J. M., Vasanawala, S. S. & Lustig, M. Coil compression for accelerated imaging with cartesian sampling. Mag. Res. in Med. 69, 571–582 (2013).
  • 41 Uecker, M. et al.

    ESPIRiT — an eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA.

    Mag. Res. in Med. 71, 990–1001 (2014).
  • 42 Bernstein, M. A. et al. Concomitant gradient terms in phase contrast MR: analysis and correction. Mag. Res. in Med. 39, 300–308 (1998).
  • 43 Busch, J., Giese, D. & Kozerke, S. Image-based background phase error correction in 4D flow MRI revisited. J. of MRI 46, 1516–1525 (2017).
  • 44 Walker, P. G. et al. Semiautomated method for noise reduction and background phase error correction in MR phase velocity data. J. Mag. Res. Imag. 3, 521–530 (1993).
  • 45 Tamir, J. I., Ong, F., Cheng, J. Y., Uecker, M. & Lustig, M. Generalized magnetic resonance image reconstruction using the Berkeley advanced reconstruction toolbox. In ISMRM W. on Dat. Sampl. and Imag. Recon. (2016).
  • 46 Wang, Z., Bovik, A. C., Sheikh, H. R. & Simoncelli, E. P. Image quality assessment: from error visibility to structural similarity. IEEE Trans. on Imag. Proc. 13, 600–612 (2004).
  • 47 Altman, D. G. & Bland, J. M. Measurement in medicine: the analysis of method comparison studies. J. R. Stat. Soc. 32, 307–317 (1983).
  • 48 Vishnevskiy, V., Walheim, J. & Kozerke, S. FlowVN: deep variational network for rapid 4D flow MRI reconstruction. CodeOcean URL (2020).
  • 49 Vishnevskiy, V., Walheim, J. & Kozerke, S. FlowVN: analysis. CodeOcean URL CodeOcean{}{}}{}{cmtt}. (2020).

1 Supplementary Material

Input:  — zero-filled k-space samples,  — undersampling mask
for to
Output: reconstructed image
Algorithm 1 Proposed variational reconstruction network model FlowVN.
=6   p-val. =8   p-val. =10   p-val. =12   p-val. =14   p-val. =16   p-val. =18   p-val. =20   p-val. =22   p-val. Image Magnitude nRMSE [%] LLR HamVN  FlowVN 1.90.3 2.30.2 1.70.2 0.0001 0.0017 2.00.3 2.50.2 1.90.2 0.0001 0.0076 2.20.3 2.80.3 2.20.2 0.0001 0.1418 2.50.3 3.00.3 2.40.2 0.0001 0.0256 2.80.3 3.30.3 2.60.3 0.0001 0.0001 3.30.4 3.60.3 2.90.3 0.0001 0.0001 3.60.4 3.80.4 3.00.3 0.0001 0.0001 3.80.4 4.00.4 3.20.3 0.0211 0.0001 4.20.5 4.10.4 3.40.3 0.9551 0.0001 Velocity Magnitude RelErr [%] LLR HamVN  FlowVN 9.41.9 12.82.5 8.72.4 0.0001 0.0225 10.52.3 13.52.8 10.02.8 0.0001 0.0790 11.82.2 14.93.0 11.42.8 0.0001 0.0385 13.22.4 16.42.9 12.42.4 0.0001 0.0003 15.12.5 18.03.5 14.02.9 0.0001 0.0001 17.42.3 20.23.2 15.53.3 0.0001 0.0001 19.02.7 22.42.9 16.92.9 0.0001 0.0001 19.93.4 23.14.1 17.53.4 0.0001 0.0001 21.92.8 25.33.6 19.23.3 0.0001 0.0001 Mean Flow AngErr [deg] LLR HamVN  FlowVN 7.01.8 9.62.6 6.91.9 0.0001 0.1956 7.92.0 10.52.9 8.02.3 0.0001 0.2228 8.92.2 11.53.1 9.22.5 0.0001 0.0204 10.22.7 12.73.4 10.32.8 0.0001 0.0225 11.32.9 13.63.8 11.43.2 0.0001 0.1166 12.23.1 14.34.0 12.33.4 0.0001 0.7785 13.33.4 15.34.2 13.03.6 0.0001 0.4799 13.73.7 15.54.4 13.63.9 0.0001 0.0073 14.83.6 16.14.5 14.33.8 0.0001 0.0005
Table 2: Comparison of reconstruction errors ( standard deviation) of 7 retrospectively undersampled acquisitions of healthy subjects for different acceleration factors . For each row the best performing method is highlighted in bold. P-values are obtained by comparing error metrics distributions of LLR and corresponding VN architectures using the two-sided Wilcoxon signed-rank test.
Figure 6: Retrospective ablation study. a, Comparison of VN architectures and modifications proposed in this paper, using retrospectively undersampled data from 7 healthy volunteers for different acceleration factors . HamVN-3D performs 3D filtering of dimensional data (assuming that data is fully sampled in the dimension). Modifications are accumulated as presented in the legend and, eventually, yield the FlowVN configuration. Undersampling-dependent modulation improves reconstruction quality only for higher acceleration factors. Note that the layer-wise exponential weighting of the reconstruction loss does not increase model complexity and, therefore, does not improve reconstruction accuracy for retrospective acceleration simulations. b, Velocity magnitude error maps for =16 shown at systolic peak flow.
Method Peak Velocity Error[%] p-value Peak Flow Error[%] p-value
HamVN-3D, single coil 19.3819.01 0.0001 21.3715.11 0.0001
HamVN-3D 6.4413.64 0.0001 3.73 12.12 0.0020
HamVN 3.6410.10 0.0001 1.909.69 0.3345
+linear activation functions 3.1810.35 0.0002 1.6510.70 0.3833
+unroll GD with momentum 2.959.89 0.0004 1.479.82 0.4749
+data activation functions 2.2910.12 0.0213 0.639.74 0.7957
+US-depending modulation 1.988.43 0.0645 0.4710.30 0.9272
+exp. weighting (FlowVN) 1.599.65 0.0875 0.059.79 0.9550
Table 3: Quantitative prospective comparison of VN architectures and modifications proposed in this paper using data from 7 healthy volunteers (). Peak velocities and peak flow are assessed over manually segmented aorta slices and reported relative to average values of the corresponding reference scan ( standard deviation). P-values are obtained by comparing error metrics distributions of reference PI data and the corresponding reconstruction method using the two-sided Wilcoxon signed-rank test. HamVN-3D performs 3D filtering of dimensional data (assuming that data is fully sampled in the dimension). “HamVN-3D, single coil” was trained to reconstruct each coil image separately and then to be combined using provided sensitivity maps. Modifications are accumulated as presented and, eventually, yield the FlowVN configuration.