A comparison of phase unwrapping methods in velocity‐encoded MRI for aortic flows

The phase of a MRI signal is used to encode the velocity of blood flow. Phase unwrapping artifacts may appear when aiming to improve the velocity‐to‐noise ratio (VNR) of the measured velocity field. This study aims to compare various unwrapping algorithms on ground‐truth synthetic data generated using computational fluid dynamics (CFD) simulations.


INTRODUCTION
Phase contrast MRI is an important technique for evaluating vascular hemodynamics. 1,2 This technique encodes velocity information in the data phase, so all values are restrained to the range [− , ), where the maximum encoded velocity V enc is mapped to . If V enc is set lower than the actual velocity, an image error known as "phase wrapping" occurs as the phase value is wrapped back into the range [− , ). Conversely, choosing a V enc that is too high decreases the velocity-to-noise ratio, thus lowering the quality of the data. 3 Since the range of measured velocities in MRI can be very large, the selection of the proper V enc may involve repeated scans. This makes it attractive to postprocess the acquired velocity data with phase unwrapping techniques to correct the wraps introduced by using a low V enc . This also allows the utilization of high-velocity-to-noise ratio data acquired with a lower V enc .
This can be performed in a number of different ways by assuming regularity in the spatial [4][5][6] or the temporal dimension [7][8][9] or both. 10,11 Loecher et al. 4 introduce a four-dimensional (4D) Laplacian method, and its performance is confirmed on both simulated flow and a volunteer dataset of the aorta via visual qualitative assessment from experts. The method presented in Bhalerao et al. 5 use 4D MRI data from a phantom and a single volunteer, with various V enc s and (simulated) noise levels. The method based on temporal continuity shown by Xiang 9 was tested on a single acquisition of two-dimensional (2D) slices of the aorta of a healthy volunteer and evaluated qualitatively. An improved temporal approach was shown in Untenberger et al., 7 and the evaluation used computational fluid dynamics (CFD) data as well as MRI acquisitions of a phantom and 10 healthy volunteers with 2 V enc s. A different approach is used in another paper by Loecher et al., 10 where the probability of a voxel being wrapped is computed using the gradients with neighboring voxels. This is shown to be effective on 4D acquisitions of the aorta with several V enc s evaluated visually. For magnetic resonance elastography, Barnhill et al. 11 compare three different algorithms as well as two commercial algorithms on volunteer acquisitions of various organs but not with regards to blood flow. All these methods have been shown to be effective, but there is often only a visual or qualitative evaluation of their effectiveness rather than a quantitative one. Generally only a single spatial and temporal resolution is used, which makes it difficult to establish the limitations of the algorithms. In the same manner, the sensitivity to noise is rarely investigated.
So far, there has been no thorough comparison of these different algorithms. Franco et al. 12 provide a comparison of dual-venc algorithms, but do not include single-venc approaches. For the vast majority of the algorithms, the code is also not publicly available and hence methods cannot be compared with the same datasets. Moreover, ground truth data in subjects or experimental work would involve for instance multi-venc acquisitions-in order to get high velocity-to-noise ratio perfectly unwrapped images-which require advanced scan protocols not yet available in clinical practice.
Therefore, the aim of this article is to compare four different phase-unwrapping methods using synthetic data generated using CFD simulations where a ground truth is available. We establish the advantages and drawbacks of the various methods and the cases where they are appropriate. Additionally, the code we used for each of the algorithms will be made publicly available under https://git.web.rug.nl/p305235/ Phase_Unwrapping_Comparison. Moreover, we confirm some of the findings on 2D PC-MRI datasets published previously.
The remainder of the paper is structured as follows. First, the phase wrapping problem is introduced and we present the different unwrapping methods. Then we describe the details of the different datasets and algorithms. The results of the unwrapping algorithms for the different datasets are then compared and discussed.

Problem statement
In PC-MRI, the velocity information is encoded into the phase of the transverse component of the magnetization such that the measured phase is for a selected encoding gradient M(G), where 0 is the background phase and is the gyromagnetic ratio. By taking two-phase measurements 1 , 2 with different encoding gradients, the velocity can be retrieved as where the velocity encoding gradient V enc is defined as . This parameter is manually set for the PC-MRI acquisition.
Since the phase is limited to the interval [− , ), the measured velocity is limited to the interval [−V enc , V enc ). The phase, being the angle of a vector, is a cyclical entity and will "wrap around" to the negative part of the interval if it exceeds the upper limit, and vice versa. As a result, if the true velocity is higher than the V enc , the measured velocity will be where k ∈ Z is determined by how far the true velocity u true exceeds V enc , and is unknown if only the measured velocity is known. Analogously, the estimated phase will be w = true − 2k ⋅ where the true phase may now exceed the interval [− , ). This artifact is known as velocity aliasing or phase wrapping and occurs in other MRI procedures and different fields as well. Resolving these wraps to reconstruct the true phase or velocity is called unwrapping.
Out of the available methods, we review three algorithms solving the 4D phase unwrapping problem: one relying on temporal smoothness, one relying on spatial smoothness for solving a Laplacian equation over the entire domain, and one using both temporal and spatial smoothness to iteratively unwrap the picture pixel-by-pixel by assigning each pixel a probability of being wrapped. To avoid confusion, we will use "3D + T" (three spatial dimensions plus time) for the 4D problem, and analogously "2D + T," "3D," and "2D."

2.2
Phase unwrapping methods

Temporal unwrapping
Temporal phase unwrapping has first been introduced in Reference 8. This method is based on the assumption that at each voxel the phase only varies slowly in time or that the temporal resolution is high enough. Inter-voxel dependency is neglected. Given a time series of N phase maps w (x, t ), = 1, … , N, the differential phase maps are computed as: According to the aforementioned assumption, these differential maps cannot contain any phase wraps of their own. Therefore any value |D (x)| greater than has to be the result of a phase wrap occurring in one of the phase maps. To regain the "correct" differential value, it is sufficient to wrap the differential phase maps into the range [− , ). Once the wrapping-free differential phase maps have been computed, the unwrapped phase maps can be calculated by integrating over the differential maps, starting at a reference time frame. This requires the existence of a timeframe that does not have any wrapped pixels, which will have to be manually selected. When trying to measure blood velocities, the cardiac cycle contains times with very low velocities, for example, in diastole, so such a reference timeframe should exist when the temporal resolution is high enough. 9 The unwrapped phase u with a reference timeframe w (x, t ref ) at time t ref can then be described as This algorithm is simple to implement and efficient in terms of computational effort. However, it requires a sufficiently high temporal resolution as well as a wrap-free reference frame, which has to be manually selected. The method can unwrap multiple wraps (which occur if the V enc is less than 50% of the maximal velocity), as long as they occur in separate timesteps, that is, a wrapped pixel in one timeframe being wrapped again in a later timeframe. This method is sensitive to noise and may falsely unwrap pixels due to noise. Movement of the vessel of interest also negatively impacts the results for this method as it breaks the continuity of the pixels between the timeframes.

3D + T gradient-based unwrapping
In gradient-based unwrapping methods 10 each pixel is assigned a probability that it is wrapped based on the gradients with the neighboring pixels, and all pixels with a probability higher than a certain threshold are unwrapped. This method evaluates the data pixel-by-pixel and often requires several consecutive iterations. The method assumes smoothness in both space and time. The probability that a pixel is wrapped is estimated by the gradients of the neighboring pixels to the evaluated pixel via where N(x) is the set of the spatial neighbors of the pixel x and N(t ) is the set of the temporal neighbours of the time position t . The constants w s , w t are weighting factors. It is recommended to set w s to 1 for the spatial dimension and w t to 2.5 for the temporal dimension. 10 Each step is detailed in Algorithm 1.
The algorithm uses a dual threshold approach, where the algorithm first completes a number of iterations with a low threshold in order to break up larger regions of Algorithm 1.
3D+T gradient-based unwrapping algorithm Given a wrapped phase image w , a set of thresholds r low and r high , and iteration numbers n low , n high : Repeat n low times: Repeat again n high times with r high instead of r low Return wrapped pixels, and then uses a higher threshold for further iterations to accurately unwrap the remaining pixels. The suggestion for the optimal thresholds are r low = 0.32 and r high = 0.75, for 50 iterations each. 10 This algorithm is very simple to implement and requires little computational effort per iteration. However, depending on the number of iterations, the time required may increase above that of other methods. This method can also only unwrap single wraps and therefore requires that V enc > 50% of the maximal velocity to unwrap successfully.

Laplacian unwrapping
Laplacian unwrapping is a popular single-step technique that can be used with 2D, three-dimensional (3D), 2D+T, and 3D+T MRI data. Starting from the relationship the method assumes the continuity and differentiability of and w and aims to first find the Laplacian of the true phase. Specifically, by first showing the relation with Δ and ∇ the spatio(-temporal) Laplacian and gradient operators, respectively, and combining it with Equation (7), the following equation is obtained: Simplifying the right-hand side leads to Therefore, the unwrapped phase is computed by solving the problem: whereΔ denotes the Laplacian operator but applying some set of boundary conditions. In the case of a spatiotemporal Laplacian operator for 2D+T or 3D+T datasets, it is necessary to introduce a scaling constant for the temporal dimension. In practice, this scaling is quite robust and was simply set to 1 s 2 ∕cm 2 , as it was in Reference 4. Note that, in discrete data, the Laplacian operation on the right-hand-side's terms is approximated by the sum of the second derivatives along each dimension.
The results differ depending on the chosen boundary conditions-Periodic, Dirichlet, or Neumann. In this work, periodic boundary conditions were used on all boundaries, as in the implementation provided by Reference 4.
The Laplacian method can be computationally expensive, depending on the method used to apply the inverse Laplacian and the size of the input data, but it can unwrap nested wraps and it is reported to handle noise well since it provides a simultaneously smoothing field.

2.2.4
Optimal multiple motion encoding Optimal multiple motion encoding (OMME) 13 is an unwrapping technique utilizing data acquired with multi- The results of this method are most robust to noise for V enc values with a proportion of 1 2 This method relies on having multiple acquisitions with different V enc values and is not directly comparable to single-V enc methods which only require one acquisition. However, it is currently the best-performing method that the authors are aware of, and its results are therefore included as a "gold-standard" reference for achievable unwrapping results.

Synthetic data
A computational fluid dynamics simulation in physiological regimes was generated as detailed in the supporting information. Let u true be the velocity field generated by this simulation. Then the MRI magnetization measurements were modeled in the following way: where 0 is the background phase, C corresponds to the magnitude, and is a complex Gaussian measurement noise with a mean of zero. PC-MRI requires an additional measurement to recover an estimate of u true since the background phase is unknown. This can by done by turning off the motion encoding, resulting in the magnetization Now the velocity can be estimated from the two magnetizations via the equation These synthetic measurements are interpolated into a tetrahedral box mesh placed around the aorta geometry.
We consider two different resolutions for the box mesh: 1.5 × 1.5 × 1.5 mm 3 and 2.5 × 2.5 × 2.5 mm 3 . The reference solution is undersampled in time to a resolution of 1 = 0.03 s or 2 = 0.06 s, leading to a total of 28, respectively, 14 time steps. Two different noise levels of 15 dB or 12 dB in the complex magnetization was applied, and three V enc s were set to 120%, 60%, and 30% of the maximum velocity, resulting in 192, 96, and 48 cm/s, respectively. Only the velocity in the foot-head direction was acquired.
Next to the 3D+T measurements, we also generated synthetic 2D+T "slice" measurements by selecting a single 2D slice perpendicular to the foot-head direction and interpolating the velocity into a rectangular hexahedral mesh with a spatial resolution of 2.5 × 2.5 mm 2 (Figure 1). The subsampling in time as well as the noise levels remained the same. The chosen V enc s were 155, 77, and 38 cm/s, corresponding to 120%, 60%, and 30% of the maximum velocity, respectively.
Fifty independent realizations of each noise level per V enc and spatial resolution were made. Also, 50 realizations of two additional noise levels (signal-to-noise ratio [SNR] 9 and SNR 6) were acquired on the 3D geometry for a spatial resolution of 2.5 × 2.5 mm 2 and a V enc of 60% of the maximal velocity. The results in Section 4 show the averages and SDs over these 50 realizations.

In vivo data
We used a database of 2D+T PC-MRI of 26 subjects originally reported in Reference 13. Measurements were acquired in a 1.5T Achieva scanner with a five-element cardiac coil. The following scan parameters were used: flip Reference solution of the blood flow from the CFD simulation and meshes used for simulating measurements. (A) CFD velocity solution at t = 0.24 s; (B) 3D aorta and slice mesh used for the measurements. angle = 15 • , repetition time = 5.5 ms, echo time = 3.7 ms, matrix size = (320, 232). The acquired data are slices perpendicular to the ascending aorta just above the valsalva sinus, in 2D plus time. The spatial resolution was 1 mm × 1 mm and the temporal resolution ranged from 35 to 48 ms depending on the patient's heart rate. Three different V enc s were used: 150, 75, and 50 cm/s. The MRI images were segmented using MATLAB so that nonzero entries are only present in the ascending and descending aorta.
To investigate the noise sensitivity, the experiments were repeated with a Gaussian noise with a zero mean and a SD of 0.15 added to the phase of the segmented measurements.

Implementation
All algorithms were implemented in both Python and Matlab. The codes are included as Data S1, including the input CFD data and scripts generating the results. For the temporal method, the reference timeframe was selected at diastole, where the flow velocity is low enough to ensure that the frame is aliasing-free. For the synthetic data, this frame was the first time frame. For the volunteer data, the last timeframe was used instead. The selection of the reference timeframe does not change the results if it is aliasing-free.
For the 3D+T gradient method, we used the threshold values suggested in Reference 10, that is, t low = 0.32 and t high = 0.75. Their work suggests using 50 iterations of either threshold, but we adjusted this value empirically to 10 iterations of either threshold.
The implementation of the Laplacian method is a translation of the code provided in Reference 4 from MAT-LAB to Python. This provides both an implementation in spatial as well as spatiotemporal dimensions.

Error assessment
We compare the unwrapping methods on the error norm where u ref is the reference solution.
To assess the statistical significance of the difference between the unwrapping results, a Friedman test with a post-hoc Wilcoxon Signed-Rank test was used, with p < 0.05 being considered statistically significant after a Bonferroni correction to account for multiple comparisons. The statistical analysis was performed in Python.

Computing times
The processing times of the studied algorithms, running on a notebook with 16 GB RAM and an Intel Core i5 processor, are shown in Table 1.

3D+T data Figures 2 and 5 show the errors and SD for each unwrap-
ping method for all considered spatial and temporal resolutions for the Aorta mesh, using a V enc of 60% of the maximal velocity. Also included are the errors for the wrapped data and unwrapped with OMME 14 using V enc1 = 120% and V enc2 = 60% or 30% of V max . In Figure 2A, it can be seen that all unwrapping methods lead to an improvement in the error over the original data. The Temporal method is very close to the baseline OMME solution, whereas the other methods perform worse. The 3D+T gradient method has a significantly higher error with a larger SD than the other methods. The pattern persists with the increase in noise for SNR 12.
Using the finer grid in Figure 2B, all methods except 3D+T gradient perform significantly better, with the spatiotemporal Laplacian now achieving similar errors to the temporal method for both noise levels. Additionally, the SD of the errors has decreased.

T A B L E 1
Processing times of the algorithms for one realization on each dataset.

F I G U R E 3
Examples of the unwrapped velocity on the Aorta mesh (CFD simulation), sliced at z = 5.7 cm, with signal-to-noise ratio 12, t = 60 ms, at t = 0.18, on the coarse mesh with spatial resolution h = 2.5 mm. (A) Wrapped; (B) optimal multiple motion encoding; (C) temporal; (D) spatial Laplacian; (E) spatiotemporal Laplacian; (F) 3D+T gradient.
For SNR 12, the differences between the temporal results and the OMME results were not statistically significant (p = 0.5 for h = 2.5 mm, and p = 0.11 for h = 1.5 mm), indicating that the temporal method performs at the same level as OMME with regards to both means and variance in these cases. All other methods were pairwise statistically significantly different (p < 0.05).
For the lower temporal resolution of t = 60 ms in Figure 2C, the temporal method has a significantly higher error and larger SD as compared to t = 30 ms. Also, the method's sensitivity to noise appears to have increased. The other methods are not impacted significantly by the change in temporal resolution. As a result, the spatiotemporal Laplacian now has similar errors to the temporal method and lower ones for SNR 12. This pattern persists in Figure 2D with the higher spatial resolution. This improves the error values for all methods except the temporal method, leading to the spatiotemporal Laplacian having a significantly lower error and standard deviation compared to the temporal method. For SNR 12, the spatial Laplacian outperforms the temporal method as well.
For the lower temporal resolution, the differences in the performance of the methods were in each case found to be statistically significant for each pairwise combination of methods.
The results for V enc = 30% of the maximal velocity are depicted in the lower part of Figure 2. It can be seen that the performance of all single-V enc methods worsens as compared to the higher V enc , while OMME still performs well. For the lower timestep ( t = 30 ms) in Figure 2E,F, the Temporal method outperforms the Laplacian and 3D+T Gradient methods, but remains significantly worse than OMME. For the higher timestep ( t = 60 ms, Figure 2G,H), the Temporal method and the Spatiotemporal method perform worse, and all the single-V enc methods lead to only a small decrease in the error. All differences were found to be statistically significant. Figure 3 shows the results of the unwrapping methods on the aorta mesh on an example 2D slice including the

F I G U R E 4
Mean error of each method by signal-to-noise ratio on the 3D+T aorta data (CFD simulation) with V enc = 96 m/s, a spatial resolution of h = 2.5 mm, and a temporal resolution of t = 30 ms.

F I G U R E 5
Average error and standard deviation by each method for different t on a hexagonal element slice mesh (CFD simulation) taken at z = 5.7cm with V enc = 77 cm/s (60% of the maximal velocity) in ( ascending and descending aorta. It can be seen that the Temporal method and both Laplacian methods unwrap perfectly or nearly perfect, whereas the 3D+T Gradient only leads to a very partial unwrap. In Figure 4, the mean errors of the different unwrapping methods are depicted for varying noise levels, for the case of the CFD simulation of the 3D aorta geometry, with a V enc of 60% of the maximal velocity, a spatial resolution of h = 2.5 mm and a temporal resolution of 30 ms. It can be seen that the Temporal method, while performing better than other methods at low noise levels, responds significantly worse to an increase in noise. The other methods all react similarly. It can also be seen that we can expect a reduction in the error from the wrapped data from the Laplacian methods and OMME even with a SNR of 6, given this temporal and spatial resolution.

2D+T data
The results on the 2D+T slices with hexagonal elements and a V enc of 60% are depicted Figure 5A,B. Here, all methods except for the 3D+T gradient method perfom very similarly for t = 30 ms. The spatiotemporal Laplacian has a much lower SD than the other methods, especially for SNR 12. For the coarser temporal resolution, the temporal method does significantly worse than the Laplacian methods, with an error twice as large as that of the spatiotemporal Laplacian and a much larger SD. The 3D+T gradient method does not lead to a significant reduction in the error. The only methods without statistically significant differences were OMME and temporal for t = 30 ms (p = 1 and p = 0.17, respectively) as well as OMME and spatiotemporal Laplacian for SNR15, t = 30 ms (p = 0.28),

F I G U R E 6
Examples of the unwrapped velocity on a slice taken at z = 5.7 cm, with signal-to-noise ratio 12, t = 30 ms, at t = 0.18. (A) Wrapped; (B) optimal multiple motion encoding; (C) temporal; (D) spatial Laplacian; (E) spatiotemporal Laplacian; (F) 3D+T gradient.
The results with a lower V enc of 30%, hence exhibiting double wraps, are shown in Figure 5C,D. All single-V enc methods perform significantly worse compared to the higher V enc , and both Laplacian methods no longer lead to a substantial decrease in the error. All methods had statistically significantly different results from each other with the exception of the spatiotemporal Laplacian and the 3D+T Gradient in the case of t = 60 ms as well as the 3D+T Gradient and the wrapped data in the case of SNR 12, t = 30 ms.
Examples for the unwrapping results on the 2D+T data can be seen in Figure 6, with a temporal resolution of t = 30 ms and a noise level of SNR 12. It can again be observed that the temporal method unwraps most wrapped voxels, but also generates a falsely wrapped voxel with a very high velocity. The spatial Laplacian leads to a good unwrap with only a few wrapped pixels remaining, while the spatiotemporal Laplacian unwraps perfectly. However, the 3D+T gradient method only leads to a partial unwrap, with many aliased voxels remaining.

In vivo data
The errors of the unwrapping methods on each of the volunteer MRI acquisitions for V enc = 75 cm/s are depicted in Figure 7.
It can be seen that the reduction in error is lesser than for the synthetic data. Still the error is reduced by more than 50% on average after unwrapping. The temporal method and spatial Laplacian method perform similarly on this dataset, while the spatiotemporal Laplacian has the lowest error and SD of all methods. The 3D+T gradient method, as with the synthetic data, only marginally reduces the error.
For the errors on the noised data in Figure 7 show that the errors for all methods increased slightly, but the comparative performance of the methods remained the same. In both cases, none of the methods reach the same performance as OMME.
For this type of data, the differences between temporal and spatial Laplacian were not statistically significant for both the original and noised data (p = 0.79 and p = 0.35, respectively). Additionally, for the noised data, there were no statistically significant differences between the input data and the results of the 3D+T Gradient method (p = 0.88), confirming that no effective unwrapping was done. All other methods were statistically significantly different from both the input data and each other.
For the lower V enc = 50 cm/s (results in Figure 7C,D), the spatiotemporal Laplacian achieves the best results. Nonetheless, the results are considerably worse than for the higher V enc , with a higher error and SD. Especially for the noised data, only a small reduction in error is achieved, and the range of the SD shows that the error may even increase over the original data in some realizations.

F I G U R E 7
Average error and standard deviation by each method on the 2D+T PC-MRI data from 26 different subjects with V enc = 75 cm/s on the top row, V enc = 50cm∕s below.

F I G U R E 9
Examples of the unwrapped velocity on MRI data acquisition with added noise from volunteer 7 at t = 4. (A) Wrapped; (B) optimal multiple motion encoding; (C) temporal; (D) spatial Laplacian; (E) spatiotemporal Laplacian; (F) 3D+T gradient.
It was found that the Temporal and the Spatiotemporal Laplacian method did not have statistically significant differences for both noise levels (p = 0.73 and p = 0.71, respectively). All other differences were statistically significant. Figure 8 shows example pictures of the unwrapping results from volunteer 7. The Temporal method unwraps very well except for several voxels on the edges of the ascending and descending aorta, where the movement of the aorta disrupts the temporal continuity of the voxels between one timestep and the next. The spatial Laplacian method fails to unwrap a small area on the edge of the ascending aorta, but unwraps everything else successfully. The spatiotemporal Laplacian methods removes all aliased voxels. In contrast, the 3D+T gradient method removes only a subset of aliased voxels, leading to a partial, fractured unwrapping.
The results of the data with additional noise are shown in Figure 9. It can be seen that all the methods except 3D+T Gradient deliver similar results to the original data except for a few additional falsely wrapped voxels, while the unwrapping results of 3D+T Gradient are significantly worse, with no effective unwrapping done.

DISCUSSION
The results on the different datasets are consistent with each other in regards to the relative performance of the different unwrapping methods. This shows that the performance of the unwrapping algorithms on the synthetic data used here is likely a good indicator of their performance on in-vivo data.
On each of the datasets, the temporal method and the spatiotemporal method generally outperform the spatial Laplacian and the 3D+T gradient method. The 3D+T gradient method fails to perform adequately, as the resulting images are still heavily aliased, and is additionally more computationally expensive than the other methods due to the number of iterations. The spatial Laplacian offers no advantage over the spatiotemporal Laplacian and performs worse than the temporal method in almost all cases as well, with the exception of the PC-MRI acquisitions, where there is no statistically significant difference between the spatial Laplacian and the temporal method. It is only slightly less computationally expensive than using the higher-dimensional Laplacian and significantly more expensive than the temporal method.
With a small timestep, the temporal method achieves equal or better results than the spatiotemporal Laplacian method, especially on a coarser mesh. It is also significantly cheaper to compute than the spatiotemporal Laplacian. However, it is more sensitive to noise, and has a much larger SD, especially on the 2D+T datasets. As a result, the spatiotemporal Laplacian outperforms the temporal method on datasets with a finer grid, a larger timestep, or more noise. It is overall the most reliable method and often generates completely unaliased images.
An interesting aspect is that the spatiotemporal Laplacian method does not drastically outperform the temporal method on the full 3D+T Aorta datasets, despite having the advantage of using an additional spatial dimension as compared to the PC-MRI slice datasets.
It is notable that all the methods seem to be limited by the V enc . With a V enc of less than 50% of the maximal velocity, the performance of the unwrapping methods (except for OMME) decreases severely, likely due to the presence of double wraps. It is therefore not recommendable to use a very low V enc to take advantage of the low noise, as the error will still be higher than with a successfully unwrapped higher-V enc acquisition.
A limitation of this study is that the algorithms were not tested on full 3D+T in-vivo MRI acquisitions, nevertheless the synthetic data presented here is realistic in terms of spatiotemporal resolution and noise.

CONCLUSION
In this work we compared four different algorithms to correct the aliased data and retrieve correct velocity values. Using two different synthetic datasets and in-vivo data, we have shown the advantages and disadvantages of the algorithms and have come to the conclusion that the spatiotemporal Laplacian described in Reference 4 outperforms the others in terms of performance and reliability with acceptable computational cost. The temporal method is computationally cheaper, but more sensitive to noise and low temporal resolution, and is therefore only preferable if the temporal resolution is high compared to the spatial resolution or computational capacities are highly limited. The other two algorithms (spatial Laplacian and 3D+T gradient) offer no advantages compared to these two. The presented single-V enc methods however do not reach the reliable performance of OMME, and are not effective for a V enc of lower than 50% of the maximal velocity.