Background phase correction for quantitative cardiovascular MRI
09911062 ยท 2018-03-06
Assignee
Inventors
Cpc classification
G01R33/56518
PHYSICS
G01R33/56383
PHYSICS
G06V10/42
PHYSICS
G06T11/005
PHYSICS
International classification
G01R33/565
PHYSICS
Abstract
Systems and methods of correcting eddy current-induced background phase (EC-BP) in magnetic resonance imaging (PC-MRI) data. The method includes acquiring a slice of interest (SOI) at a first table position using a magnetic resonance imaging (MRI) scanner, the slice of interest having a predetermined imaging orientation and being acquired having predetermined gradient waveforms; acquiring at least one additional slice at a second table position using the MRI scanner, the at least one additional slice having a same imaging orientation as the slice of interest and being acquired using the same gradient waveforms as the slice of interest; determining time-averaged phase maps from the slice of interest and the at least one additional slice; determining a correction map from the time-averaged phase maps; and correcting a background phase (BP) of the slice of interest using the correction map.
Claims
1. A method of correcting eddy currents-induced background phase (EC-BP) in magnetic resonance imaging data, comprising: acquiring phase-contrast data from a slice of interest (SOI) using a magnetic resonance imaging (MRI) scanner, with a user-specified table position and SOI position and orientation; the SOI position and orientation, with respect to the magnet coordinate system, being controlled by varying a center frequency of an applied radio-frequency (RF) excitation pulse in combination with a variation of a slice select gradient amplitude and direction; acquiring phase-contrast data from at least one additional slice, the at least one additional slice having the same position and orientation, with respect to the magnet coordinate system, as the SOI and being acquired using the slice select gradient amplitude and direction, the at least one additional slice imaging a physically different location of the subject, achieved by either changing the table position with respect to the position used for the SOI acquisition or by physically moving the subject within the scanner; performing image reconstruction from k-space to image domain for the SOI and the at least one additional slice; subtracting a phase of a reconstructed additional slice from a phase of a reconstructed SOI to extract phase difference maps for the SOI and the at least one additional slice; determining a correction map by: processing the phase difference maps from both the SOI and the at least one additional slice; and estimating EC-BP in accordance with a first relationship: .sup.k1 represents the k coefficients of a 2D polynomial, wherein B and z are obtained by vertical concatenation of N temporally-averaged slices, such that B
[A.sub.1.sup.T|A.sub.2.sup.T| . . . |A.sub.N.sup.T].sup.T and z
[y.sub.1.sup.T|y.sub.2.sup.T| . . . |y.sub.N.sup.T].sup.T, and Bxz.sub..sup.2
(Bxz).sup.1(Bxz), wherein each column of matrix A
.sup.Mk defines the spatial distribution of one of the polynomial terms, wherein vector y
.sup.M1 is the measured temporally-averaged phase map, wherein
.sup.NMNM is a covariance matrix constructed by computing inner products of mean-subtracted temporal profiles across all pixel pairs in z; and correcting the BP of the SOI using the correction map.
2. The method of claim 1, wherein volumetric field of view is used instead of 2D slices.
3. The method of claim 1, wherein multiple encoding directions are used and correction is applied to each encoding direction.
4. The method of claim 1, wherein the background arises from sources other than eddy currents.
5. The method of claim 1, further enforcing an l.sub.1-norm penalty on polynomial coefficients to transform the first relationship as follows:
6. The method of claim 5, wherein x.sub.1 is replaced with Rx.sub.p with p0 and R is a linear or nonlinear transform.
7. The method of claim 1, further comprising applying a Cholesky decomposition, =LL, {tilde over (B)}L.sup.1B, {tilde over (z)}
L.sup.1z, and enforcing an l.sub.1-norm penalty on polynomial coefficients yielding:
8. The method of claim 7, wherein x.sub.1 is replaced with Rx.sub.p with p0 and R is a linear or nonlinear transform.
9. The method of claim 1, wherein a correction map, {circumflex over ()}.sub.mSAP, is then obtained by {circumflex over ()}.sub.mSAP=A.sub.i {circumflex over (x)}.sub.mSAP, wherein A corresponds to the SOI.
10. The method of claim 1, wherein only one slice (SOI) is collected and processed to find the correction map.
11. The method of claim 1, wherein only additional slices are used to find the correction map for SOI.
12. The method of claim 1, wherein the data from individual frames, before time-averaging, is used to construct B, z, {tilde over (B)}, and {tilde over (z)}.
13. The method of claim 1, wherein polynomials of different orders are used.
14. The method is claim 1, wherein functions other than polynomial functions are used for fitting the background phase.
15. The method of claim 1, further comprising constructing in a sparse format by setting entries with small values equal to zero.
16. The method of claim 1, further approximating by setting all off-diagonal entries to zero.
17. The method of claim 1, replacing with the identity matrix and only considering the pixels in z which have small temporal standard deviation.
18. The method of claim 1, further comprising removing rows from B and z that belong to a region with nearly zero spin-density.
19. The method of claim 1, wherein a solution to one of a first relationship, a second relationship or a third relationship, respectively below is used as prior information for removal of concomitant gradients and gradient non-linearity: .sup.k1 represents the k coefficients of a 2D polynomial, wherein B and z are obtained by vertical concatenation of N temporally-averaged slices, such that B
PI [A.sub.1.sup.T|A.sub.2.sup.T| . . . |A.sub.N.sup.T]T and z
[y.sub.1.sup.T|y.sub.2.sup.T| . . . |y.sub.N.sup.T].sup.T, and Bxz.sub..sup.2
(Bxz).sup.1(Bxz), wherein each column of matrix A
.sup.Mk defines the spatial distribution of one of the polynomial tennis, wherein vector y
.sup.M1 is the measured temporally-averaged phase map, wherein
.sup.NMNM is a covariance matrix constructed by computing inner products of mean-subtracted temporal profiles across all pixel pairs in z, wherein enforcing an l.sub.1-norm penalty on polynomial coefficients transforms the first relationship to the second relationship, and wherein a Cholesky decomposition, =LL, {tilde over (B)}
L.sup.1B, {tilde over (z)}
L.sup.1z, is applied to the second relationship to yield the third relationship.
Description
BRIEF DESCRIPTION OF THE DRAWINGS
(1) The components in the drawings are not necessarily to scale relative to each other. Like reference numerals designate corresponding parts throughout the several views.
(2)
(3)
(4)
(5)
(6)
(7)
DETAILED DESCRIPTION
(8) Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. Methods and materials similar or equivalent to those described herein can be used in the practice or testing of the present disclosure. While implementations will be described for remotely accessing applications, it will become evident to those skilled in the art that the implementations are not limited thereto, but are applicable for remotely accessing any type of data or service via a remote device.
INTRODUCTION
(9) Disclosed herein is an EC-BP correction scheme called multi-slice acquisition and processing (mSAP) where additional slices from the subject's own body are used in lieu of a static phantom. In mSAP, in addition to the slice of interest (SOI), data is also collected from at least one extra slice in the same imaging orientation using the same gradient waveforms but with different table position. As such, all slices collected under these conditions should have the same background phase; therefore, by leveraging information from additional slices with abundant static tissue, mSAP circumvents the uncertainties associated with available methods that utilize data from a single slice.
(10) The robustness of the disclosed method is strengthened by imposing l.sub.1 regularization to polynomial regression, effectively eliminating the need to precisely prescribe the polynomial order in advance.
(11) The method of the present disclosure provides a novel solution to EC-BP where information from additional anatomical slices is leveraged to perform accurate and robust correction of BP in the slice of interest. For cardiac PC-MRI, there is a lack of static tissue around the heart and great vessels; therefore, fitting performed on the SOI alone is prone to errors. If selected appropriately, HS can provide complementary information that is missing from the SOI. When processed jointly, SOI and HS can yield more accurate results than possible by fitting the SOI alone.
(12) The conventional method to correct for EC-BP is to divide the image into static and dynamic regions and perform polynomial fitting on the pixels in the static regions. The fitted polynomial surface is then subtracted from the phase map. The performance of this approach relies on user selected thresholds on magnitude and temporal variation of phase to segment the image. In contrast, the method(s) of the present disclosure use a generalized least squares method that employs the covariance matrix (or its approximation) to de-correlate and normalize the errors across the pixels. This approach not only deemphasizes regions where phase has high temporal variance but also take into account the independence of information, without using any tuning parameters.
(13) The polynomial fitting methods are sensitive to not only the extent of the static tissue in the close vicinity of the vessel but also to the selection of the polynomial order. Although data from multiple slices adds robustness to higher order fitting, one is still left with the important decision of selecting the polynomial order. For mSAP, a high order polynomial fitting is used with Lasso regression which enforces an l.sub.1-norm penalty on the polynomial coefficients. This constraint prompts unnecessary coefficients to be zero and thus effectively reduces the number of variables upon which the fitting is dependent. Although the l.sub.1-norm penalty has been extensively applied for compressive recovery of MRI from undersampled data, this has not been used as a safeguard against over-fitting for BP correction.
(14) The techniques herein correct for EC-BP by jointly processing data from anatomically different slices collected at the same position relative to isocenter using identical gradient waveforms. The BP-induced errors are reduced to less than 0.4%, resulting in an error of less than 5% in flow volume.
(15) Methodology Development
(16) The following presents the algorithmic details of mSAP for EC-BP correction.
(17) Joint Processing of Multi-Slice Data
(18) In PC-MRI, velocity information is encoded into the phase of the complex-valued image. To cancel the contributions from other sources, e.g., due to sensitivity maps of receive coils or B.sub.0 inhomogeneity, k-space at each cardiac phase is sampled twice, with two different values of VENC. Typically, the first step is to reconstruct complex image seriesone for each velocity encodingfrom the multi-channel k-space data using parallel imaging methods, e.g., SENSE or GRAPPA. A phase subtraction between the two differently encoded image series followed by scaling with VENC generates a spatiotemporal velocity map.
(19) Using the knowledge of the gradient waveforms and deviation from linear behavior, corrections can be applied to remove for the effects of concomitant gradients and gradient non-linearity. Since the velocity encoding gradient varies between the two encodings, the eddy currents generated by those two waveforms are different and thus do not completely cancel upon phase subtraction. Therefore, even after applying the correction for concomitant gradients, the velocity map is still contaminated with EC-BP. The EC-BP distribution has complicated dependence on gradient waveforms and is difficult to estimate even when the gradient waveforms are known precisely. In addition, because the EC-BP may not have long-term reproducibility, precomputing correction maps for EC-BP would be ineffective.
(20) The most commonly employed method for EC-BP is based on the polynomial fitting of the phase map, which can be written as:
(21)
(22) where vector x.sup.k1 represents the k coefficients of a 2D polynomial, each column of matrix A
.sup.Mk defines the spatial distribution of one of the polynomial terms, vector y
.sup.M1 is the measured temporally-averaged phase map, and P is a row-pruning operator (P
.sup.NM with N<M) based on user defined thresholds on the temporal variance of phase map and the magnitude of the complex image for which y is the phase. Alternatively, P is a diagonal weighting operator (P
.sup.MM) that deemphasizes pixels with high temporal variance. The resulting correction map, {circumflex over ()}, is then obtained by {circumflex over ()}=A{circumflex over (x)}.
(23) For multi-slice processing, x is estimated based on jointly fitting data from N slices, as follows:
(24)
(25) where B and z are obtained by vertical concatenation of N temporally-averaged slices, i.e., B[A.sub.1.sup.T|A.sub.2.sup.T| . . . |A.sub.N.sup.T]T and z
[y.sub.1.sup.T|y.sub.2.sup.T| . . . |y.sub.N.sup.T]T, and Bxz.sub..sup.2
(Bxz).sup.1(Bxz). The covariance matrix,
.sup.NMNM, is constructed by computing inner products of mean-subtracted temporal profiles across all pixel pairs in z. By using Cholesky decomposition, =LL,
L.sup.1B, {tilde over (z)}
L.sup.1z, and by enforcing l.sub.1-norm penalty on polynomial coefficients, Eq. 2 becomes:
(26)
(27) Eq. 3 may be solved using Fast Iterative Shrinkage-Thresholding Algorithm (FISTA), and the value of using Morozov discrepancy principle. The resulting correction map, {circumflex over ()}.sub.mSAP, is then obtained by {circumflex over ()}.sub.mSAP=A.sub.i{circumflex over (x)}.sub.mSAP, where A.sub.i corresponds to the SOI.
(28) The processing described above, enables performance of mSAP on multi-slice PC-MRI data using, e.g., Matlab (MathWorks, Natick, Mass.). mSAP will provide tuning-free correction of EC-BP.
(29) In the above, the Cholesky decomposition for large can be computationally expensive. Should that be the case, can be alternatively constructed in a sparse format by setting the entries with small values equal to zero. This sparse approximation of would facilitate faster Cholesky decomposition. Also, rows from B and z that belong to the region with nearly zero spin-density can be removed. The elimination of these entries will further speed up the processing without having significant impact on {circumflex over (x)}.sub.mSAP.
(30) Morozov discrepancy principle for selecting is based on the assumption that variance of each element of {tilde over (B)}x{tilde over (z)} is one. In the case of model mismatch, this assumption may not be valid. Should that be the case, L-curve criterion may be used to determine . Due the small size of the problem in Eq. 3, employing L-curve criterion is a viable alternative.
(31) It is not uncommon to encounter outliers in the phase map, especially at the tissue air boundary.
(32) Implementation of Acquisition Protocol
(33) To facilitate routine use, mSAP may be implemented on a clinical scanner. The implementation includes updating the existing protocol for PC-MRI acquisition. In the updated protocol, the user will specify the number and relative displacements of HS with respect to SOI. Upon the completion of SOI localization, the scanner will automatically move the table and perform low-resolution, free-breathing scans, one for each HS. Based on the contents of HS, e.g., the extent and location of static tissue, the user will specify which HS would be collected along with SOI. Each HS will be scanned using exactly the same scan parameters and gradient waveforms as the SOI; e.g., if the SOI is acquired as a segmented k-space acquisition during breath-hold, then each HS will be acquired in identical fashion. Upon the completion of multi-slice acquisition, the joint processing of the data as outlined in Task 1 will be performed on the scanner. The acquisition protocol for the mSAP method of the present disclosure will add less than 5 seconds to the overall data processing time.
(34) Thus, the above demonstrates of performance of multi-slice acquisition and processing on the scanner. By streamlining the workflow, the overhead of performing multi-slice acquisition will be minimized.
(35) Preliminary Data (Phantom Studies)
(36) This study was performed to determine whether the information from HS can be used to correct the BP in SOI if all the acquisitions are performed using the same gradient waveforms. A CardioFlow 5000 MR flow pump was used to generate pulsatile flow. The pump was connected to a flexible pipe that was bent into a u-shape such that two sections of the pipe were aligned in parallel inside the magnet. The imaging plane orientation was selected such that in-flow and return-flow were measured simultaneously. To mimic in vivo static tissue, a total of 9 water bottles of various sizes were placed around the pipe. Ten parallel slices were imaged, each at the isocenter of the magnet by changing the table position. The distance between adjacent slices was fixed at 2 cm. The phase maps from those ten slices are shown in
(37) One of the ten slices in
(38) As shown in Table 1, background correction in SOI* is ineffective, as evident by 16.9% discrepancy in stroke volume compare to SOI. To test mSAP, we jointly processed SOI* and HS (Eq. 3). The process was repeated individually for all HS. In this preliminary work, the covariance matrix, , was approximated by a diagonal matrix with entries indicating the temporal variance at each pixel. As evident from Table 1, all HS slices, when included individually, improved the flow quantification in SOI*, which affirms the key assumption made for mSAP: all slices collected under the same gradient waveforms but with different table positions have the same (or very similar) background phase. By employing information from a single HS, the error in the computation of stroke volume was reduced to 2.62% or less. Inclusion of more than one slice did not change the results significantly (data not shown), implying that the information from one HS may be adequate to suppress EC-BP.
(39) TABLE-US-00001 TABLE 1 Stroke volume (mL) in the two legs of the pipe (+ve for in-flow and ve for return-flow). Data used for Stroke Volume Correction fitting (% error) method SOI 34.3 34.4 Ground truth SOI* 33.1 (3.50%) 40.2 (16.9%) Traditional SOI* + HS1 34.2 (0.29%) 33.9 (1.45%) mSAP SOI* + HS2 34.3 (0.00%) 34.6 (0.58%) mSAP SOI* + HS3 34.4 (0.29%) 34.6 (0.58%) mSAP SOI* + HS4 34.5 (0.58%) 34.7 (0.87%) mSAP SOI* + HS5 34.2 (0.29%) 34.0 (1.16%) mSAP SOI* + HS6 34.5 (0.58%) 34.8 (1.16%) mSAP SOI* + HS7 34.7 (1.17%) 33.9 (1.45%) mSAP SOI* + HS8 34.9 (1.75%) 33.5 (2.62%) mSAP SOI* + HS9 34.6 (0.87%) 33.9 (1.45%) mSAP SOI indicates when the entire slice was used for fitting (ground truth), SOI* indicates when the static signal proximal to the pipes was removed to mimic cardiac imaging, SOI* + HSn indicates when the data from SOI* and HSn was jointly processed, where 1 n 9 identifies the specific HS used in fitting.
(40) Preliminary Data (Volunteer Imaging)
(41) A human volunteer was imaged on a 3T (Siemens Skyra) system, equipped with 18-channel coil array. Other imaging parameters were: FOV: 380260 mm.sup.2, matrix size: 192132, TE/TR: 2.94/5.0 ms, acceleration: 2, reconstruction: TGRAPPA, slice thickness: 6 mm, VENC: 100 or 150 cm/s depending on the view. The acquisition was performed using two different views to measure the flow (Qp) in MPA and the flow (Qs) in AA. With reference to
(42) TABLE-US-00002 TABLE 2 BP correction based on different combinations of slices; the resulting Qp/Qs are reported. Data used for Correction fitting Qp/Qs method uncorrected 1.112 None MPA: SOI + HS1 1.054 mSAP AA: SOI + HS1 MPA: SOI + HS1 1.046 mSAP AA: SOI + HS2 MPA: SOI + HS2 1.061 mSAP AA: SOI + HS1 MPA: SOI + HS2 1.053 mSAP AA: SOI + HS2
(43) Outcome
(44) mSAP adequately corrects the background phase, with the acquisition of only a single helper slice, and the inclusion of additional HS to result in only marginal additional improvement. The BP errors are below 0.4% of the dynamic range and the resulting miscalculation in flow volume to be below 5%.
(45) In the above, the gradient waveforms used to acquire SOI and HS may differ due to heart rate variations. The transitory effects of incrementing the phase encoding gradient with each R-wave, however, should be negligible. If this is not the case, it will be possible to record the triggering waveform from the SOI scan and use it to collect HS data.
(46) Although the covariance matrix used in mSAP is built from the spatiotemporal phase series, the fitting processing itself is applied to the time-averaged phase map. There may be small temporal fluctuations in BP across the cardiac cycle. If BP varies temporallyas observed by the disparity in fitting accuracy across different framesthe individual frames may be re-fitted by using the fitting obtained from time-averaged data as a prior, thereby allowing the BP of individual frames to vary slightly from the time-averaged value.
(47) It should be understood that the various techniques described herein may be implemented in connection with hardware or software or, where appropriate, with a combination of both. Thus, the methods and apparatus of the presently disclosed subject matter, or certain aspects or portions thereof, may take the form of program code (i.e., instructions) embodied in tangible media, such as floppy diskettes, CD-ROMs, hard drives, or any other machine-readable storage medium wherein, when the program code is loaded into and executed by a machine, such as a computer, the machine becomes an apparatus for practicing the presently disclosed subject matter. In the case of program code execution on programmable computers, the computing device generally includes a processor, a storage medium readable by the processor (including volatile and non-volatile memory and/or storage elements), at least one input device, and at least one output device.
(48) Although the subject matter has been described in language specific to structural features and/or methodological acts, it is to be understood that the subject matter defined in the appended claims is not necessarily limited to the specific features or acts described above. Rather, the specific features and acts described above are disclosed as example forms of implementing the claims.
BIBLIOGRAPHY
(49) 1. Weiss, M. B. et al. Myocardial blood flow in congestive and hypertrophic cardiomyopathy: relationship to peak wall stress and mean velocity of circumferential fiber shortening. Circulation 54, 484-494 (1976). 2. Clark, D. M., Plumb, V. J., Epstein, A. E. & Kay, G. N. Hemodynamic effects of an irregular sequence of ventricular cycle lengths during atrial fibrillation. J. Am. Coll. Cardiol. 30, 1039-1045 (1997). 3. Aurigemma, G. et al. Abnormal left ventricular intracavitary flow acceleration in patients undergoing aortic valve replacement for aortic stenosis. A marker for high postoperative morbidity and mortality. Circulation 86, 926-936 (1992). 4. Taylor, T. W. & Yamaguchi, T. Three-dimensional simulation of blood flow in an abdominal aortic aneurysmsteady and unsteady flow cases. J. Biomech. Eng. 116, 89-97 (1994). 5. Powell, A. J. & Geva, T. Blood flow measurement by magnetic resonance imaging in congenital heart disease. Pediatr. Cardiol. 21, 47-58 (2000). 6. Galderisi, M. et al. Recommendations of the European Association of Echocardiography How to use echo-Doppler in clinical trials: Different modalities for different purposes. Eur. J. Echocardiogr. 12, 339-353 (2011). 7. Moran, P. R. A flow velocity zeugmatographic interlace for NMR imaging in humans. Magn. Reson. Imaging 1, 197-203 (1982). 8. Pelc, N. J., Herfkens, R. J., Shimakawa, A. & Enzmann, D. R. Phase contrast cine magnetic resonance imaging. Magn. Reson. Q. 7, 229-254 (1991). 9. Gatehouse, P. D. et al. Flow measurement by cardiovascular magnetic resonance: a multi-centre multi-vendor study of background phase offset errors that can compromise the accuracy of derived regurgitant or shunt flow measurements. J. Cardiovasc. Magn. Reson. 12, 5 (2010). 10. Kilner, P. J., Gatehouse, P. D. & Firmin, D. N. Flow measurement by magnetic resonance: a unique asset worth optimising. J. Cardiovasc. Magn. Reson. 9, 723-728 (2007). 11. Lorenz, R. et al. Influence of eddy current, Maxwell and gradient field corrections on 3D flow visualization of 3D CINE PC-MRI data. Magn. Reson. Med. 72, 33-40 (2014). 12. Bernstein, M. A. et al. Concomitant gradient terms in phase contrast MR: Analysis and correction. Magn. Reson. Med. 39, 300-308 (1998). 13. Chernobelsky, A., Shubayev, O., Comeau, C. R. & Wolff, S. D. Baseline correction of phase contrast images improves quantification of blood flow in the great vessels. J. Cardiovasc. Magn. Reson. 9, 681-685 (2007). 14. Lankhaar, J. W. et al. Correction of phase offset errors in main pulmonary artery flow quantification. J. Magn. Reson. Imaging 22, 73-79 (2005). 15. Holland, B. J., Printz, B. F. & Lai, W. W. Baseline correction of phase-contrast images in congenital cardiovascular magnetic resonance. J. Cardiovasc. Magn. Reson. 12, 11 (2010). 16. Gatehouse, P. D. et al. A multi-center inter-manufacturer study of the temporal stability of phase-contrast velocity mapping background offset errors. J. Cardiovasc. Magn. Reson. 14, 72 (2012). 17. Noureldin, R. A. et al. The diagnosis of hypertrophic cardiomyopathy by cardiovascular magnetic resonance. J. Cardiovasc. Magn. Reson. 14, 17 (2012). 18. Fluckiger, J. U. et al. Left atrial flow velocity distribution and flow coherence using four-dimensional FLOW MRI: A pilot study investigating the impact of age and Pre- and Postintervention atrial fibrillation on atrial hemodynamics. J. Magn. Reson. Imaging 38, 580-587 (2013). 19. Mitchell, L., Jenkins, J. P., Watson, Y., Rowlands, D. J. & Isherwood, I. Diagnosis and assessment of mitral and aortic valve disease by cine-flow magnetic resonance imaging. Magn. Reson. Med. 12, 181-197 (1989). 20. Harloff, A. et al. In vivo assessment of wall shear stress in the atherosclerotic aorta using flow-sensitive 4D MRI. Magn. Reson. Med. 63, 1529-1536 (2010). 21. Moore, J. E., Xu, C., Glagov, S., Zarins, C. K. & Ku, D. N. Fluid wall shear stress measurements in a model of the human abdominal aorta: Oscillatory behavior and relationship to atherosclerosis. Atherosclerosis 110, 225-240 (1994). 22. Hope, T. A. et al. Comparison of flow patterns in ascending aortic aneurysms and volunteers using four-dimensional magnetic resonance velocity mapping. J. Magn. Reson. Imaging 26, 1471-1479 (2007). 23. Reeder, G. S., Currie, P. J., Hagler, D. J., Tajik, A. J. & Seward, J. B. Use of Doppler Techniques (Continuous-Wave, Pulsed-Wave, and Color Flow Imaging) in the Noninvasive Hemodynamic Assessment of Congenital Heart Disease. Mayo Clin. Proc. 61, 725-744 (1986). 24. Franois, C. J. et al. Renal arteries: isotropic, high-spatial-resolution, unenhanced MR angiography with three-dimensional radial phase contrast. Radiology 258, 254-260 (2011). 25. Stankovic, Z. et al. Normal and Altered Three-dimensional Portal Venous Hemodynamics in Patients with Liver Cirrhosis. Radiology 262, 862-873 (2012). 26. Hope, T. A. et al. Evaluation of intracranial stenoses and aneurysms with accelerated 4D flow. Magn. Reson. Imaging 28, 41-46 (2010). 27. Arenillas, J. F. et al. Progression and clinical recurrence of symptomatic middle cerebral artery stenosis: a long-term follow-up transcranial Doppler ultrasound study. Stroke. 32, 2898-2904 (2001). 28. Strandness, D. E., Schultz, R. D., Sumner, D. S. & Rushmer, R. F. Ultrasonic flow detection. Am. J. Surg. 113, 311-320 (1967). 29. Nordmeyer, S. et al. Four-dimensional velocity-encoded magnetic resonance imaging improves blood flow quantification in patients with complex accelerated flow. J. Magn. Reson. Imaging 37, 208-216 (2013). 30. Schubert, T., Bieri, O., Pansini, M., Stippich, C. & Santini, F. Peak velocity measurements in tortuous arteries with phase contrast magnetic resonance imaging: the effect of multidirectional velocity encoding. Invest. Radiol. 49, 189-194 (2014). 31. Gatehouse, P. D. et al. Applications of phase-contrast flow and velocity imaging in cardiovascular MRI. Eur. Radiol. 15, 2172-2184 (2005). 32. Frydrychowicz, A. et al. Three-dimensional analysis of segmental wall shear stress in the aorta by flow-sensitive four-dimensional-MRI. J. Magn. Reson. Imaging 30, 77-84 (2009). 33. Markl, M. et al. Estimation of global aortic pulse wave velocity by flow-sensitive 4D MRI. Magn. Reson. Med. 63, 1575-1582 (2010). 34. Dyverfeldt, P., Hope, M. D., Tseng, E. E. & Saloner, D. Magnetic resonance measurement of turbulent kinetic energy for the estimation of irreversible pressure loss in aortic stenosis. JACC Cardiovasc. Imaging 6, 64-71 (2013). 35. Tyszka, J. M., Laidlaw, D. H., Asa, J. W. & Silverman, J. M. Three-dimensional, time-resolved (4D) relative pressure mapping using magnetic resonance imaging. J. Magn. Reson. Imaging 12, 321-329 (2000). 36. Ebbers, T. et al. Higher order weighted least-squares phase offset correction for improved accuracy in phase-contrast MRI. in ISMRM 16th Annu. Meet. Exhib. 1367 (2008). 37. Eriksson, J. et al. Semi-automatic quantification of 4D left ventricular blood flow. J. Cardiovasc. Magn. Reson. 12, 9 (2010). 38. Hsiao, A. et al. Rapid pediatric cardiac assessment of flow and ventricular volume with compressed sensing parallel imaging volumetric cine phase-contrast MRI. Am. J. Roentgenol. 198, (2012). 39. Fair, M., Gatehouse, P. D., Greiser, A., Drivas, P. & Firmin, D. N. A novel approach to phase-contrast velocity offset correction by in vivo high-SNR acquisitions. J. Cardiovasc. Magn. Reson. 15, P56 (2013). 40. Lustig, M., Donoho, D. L., Santos, J. M. & Pauly, J. M. Compressed Sensing MRI. IEEE Signal Process. Mag. 25, 72-82 (2008). 41. Griswold, M. A. et al. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn. Reson. Med. 47, 1202-1210 (2002). 42. Pruessmann, K. P., Weiger, M., Scheidegger, M. B. & Boesiger, P. SENSE: Sensitivity encoding for fast MRI. Magn. Reson. Med. 42, 952-962 (1999). 43. Markl, M. et al. Generalized reconstruction of phase contrast MRI: Analysis and correction of the effect of gradient field distortions. Magn. Reson. Med. 50, 791-801 (2003). 44. Beck, A. & Teboulle, M. A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems. SIAM J. Imaging Sci. 2, 183-202 (2009). 45. Chen, W., Guo, Y., Wu, Z. & Nayak, K. MRI Constrained Reconstruction Without Tuning Parameters Using ADMM and Morozov's Discrepency Principle. in ISMRM 23rd Annu. Meet. Exhib. 3410 (2015). 46. Hansen, P. C. & O'Leary, D. P. The Use of the L-Curve in the Regularization of Discrete III-Posed Problems. SIAM J. Sci. Comput. 14, 1487-1503 (1993). 47. Tan, E. T., Brau, A. C. & Hardy, C. J. Nonlinear self-calibrated phase-contrast correction in quantitative cardiac imaging. J. Cardiovasc. Magn. Reson. 16, P349 (2014). 48. Iglewicz, B. & Hoaglin, D. How to detect and handle outliers. (ASQC Quality Press, 1993). 49. Breuer, F. A., Kellman, P., Griswold, M. A. & Jakob, P. M. Dynamic autocalibrated parallel imaging using temporal GRAPPA (TGRAPPA). Magn. Reson. Med. 53, 981-985 (2005). 50. Bernstein, M. A., Grgic, M., Brosnan, T. J. & Pelc, N. J. Reconstructions of phase contrast, phased array multicoil data. Magn. Reson. Med. 32, 330-334 (1994). 51. Dunn, O. J. Multiple Comparisons Using Rank Sums. Technometrics 6, pp. 241-252 (1964).