Computationally effective 2D and 3D fast phase unwrapping algorithms and their applications to Doppler optical coherence tomography

: We propose a simplification for a robust and easy to implement fast phase unwrapping (FPU) algorithm that is used to solve the phase wrapping problem encountered in various fields of optical imaging and metrology. We show that the number of necessary computations using the algorithm can be reduced compared to its original version. FPU can be easily extended from two to three spatial dimensions. We demonstrate the applicability of the two- and three-dimensional FPU algorithm for Doppler optical coherence tomography (DOCT) in numerical simulations, and in the imaging of a flow phantom and blood flow in the human retina in vivo . We introduce an FPU applicability plot for use as a guide in the selection of the most suitable version of the algorithm depending on the phase noise in the acquired data. This plot uses the circular standard deviation of the wrapped phase distribution as a measure of noise and relates it to the root-mean-square error of the recovered, unwrapped phase.

Doppler OCT methods can provide information similar to OCTA techniques, i.e., about the 3D architecture of active vascular systems [13,[22][23][24]. However, they can also enable the quantification of parameters characterizing the blood flow in the vessels of various organs, such as the eye [21,23], or the brain [21,[25][26][27]. Doppler OCT can therefore provide information that is important for the research and clinical diagnosis of pathologies affecting specific tissues or the entire circulatory system. However, there are a few practical limitations associated with this technique. One of the most often discussed limitations is the range of axial flow velocity values which can be measured. If the flow is too slow, it cannot be distinguished from the noise attributed to the phase stability of various components of the OCT system [19,[28][29][30][31][32][33][34]. If the flow is fast enough to introduce a phase change among the complex signals that exceeds 2π, a phase wrapping occurs that renders the estimation of the flow velocity obscured [24,32]. The same effect is encountered in electron holography, magnetic resonance imaging (MRI), radar interferometry, and other interferometric techniques [35][36][37][38]. If the flow is too fast, i.e., much faster than the imaging speed of the OCT system, interference fringe wash-out occurs [32,[39][40][41], and the flow cannot be detected. However, the severity of this effect depends on the degree of averaging of the interference fringes. In practice, it is prominent in the spectral-domain OCT systems and often negligible in the swept-source OCT systems.
Since the phase wrapping problem is commonly encountered in various imaging and detection techniques, several data analyses methods have been developed to solve it [38]. The simplest algorithm adds or subtracts 2π whenever the phase difference between consecutive phase values exceed π. Calculations of phase differences directly make this algorithm highly sensitive to noise and unpractical for use in most real life applications. More sophisticated algorithms were proposed by Goldstein [35] and Itoh [42]. These methods search for phase discontinuities by calculating phase gradients either directly or with use of Laplace operators, or by comparing phase distributions at different wavelengths. Therefore, these algorithms require numerous, time consuming calculations. To improve the efficiency of the computations, a fast phase unwrapping (FPU) method was developed by Schofield et al. [43]. The calculation of the phase derivatives with Laplace operators was replaced by Fourier transformations. This method has also been used in MRI for quantitative magnetic susceptibility mapping (QSM) [44], MRI venography [45], MRI flow estimation methods [46], magnetic resonance elastography [47], electron holography [48], and digital holographic microscopy [49].
The main reason for the failure of the phase unwrapping methods is attributed to the low signal-to-noise ratio (SNR) of reconstructed phase images. A phase unwrapping method for Doppler OCT must meet additional requirements, e.g., it must work in low SNR cases, and must also be insensitive to spatial variations of the SNR, including variations in depth (sensitivity roll-off), across the vessels (interference fringe wash-out at increasing flow velocities), or across the imaged pathology (weak OCT signal due to developing pathology). This method should also work in the presence of speckle noise and preferably operate in real time. To the best of our knowledge, to-this-date, the only methods that have been used for phase unwrapping in OCT are synthetic wavelength phase unwrapping [50] and various phase gradient analyses methods [51][52][53].
In this study, we have adapted the fast phase unwrapping (FPU) algorithm for spectraldomain Doppler OCT imaging, and demonstrate its applicability in data obtained with the joint spectral and time domain OCT (STdOCT) method [12]. We present the FPU method in a unified notation for multidimensional signals and propose a methodology to improve the calculation speed. The algorithm in the presented form has three properties that make it attractive for the OCT data processing. First, it uses the basic properties of Fourier transformations and postulates no assumptions on the input signal. Secondly, the method can be easily applied to either 2D or 3D data analyses. Lastly, it can be applied directly to OCT complex data without the need to extract the phase of the signal before its input to the algorithm. We compare the phase unwrapping performances of the FPU algorithm following its applications in numerical simulations and in data obtained from Doppler OCT or STdOCT imaging of a flow phantom which used a milk solution as the fluid flowing through a glass capillary. We demonstrate the applicability of the 2D and 3D FPU for Doppler OCT imaging of the human retina in vivo. We finally introduce a method to estimate the performance of the algorithms directly from the wrapped phases.

Fast phase unwrapping
The common observation in phase unwrapping methods is that the phase has to be corrected when an abrupt change occurs in the analyzed data (e.g., wrapping over a 2π range). One of the strategies used to correct the phase wraps is to calculate the phase derivatives and then compute the integral of the result. This leads to the estimation of the unwrapped phase, as will be explained in detail in a subsequent section. Depending on the data dimensionality, one, two, or three partial spatial derivatives have to be calculated. In the approach presented by Schofield et al. [43], a two-dimensional (2D) phase derivative calculation was performed using 2D Laplace operators. The key observation in this popular algorithm was that the use of Fourier transformations simplifies and speeds up the computations. Therefore, the name fast phase unwrapping (FPU) was used. Herein, we outline the FPU algorithm for data with an arbitrary number of spatial dimensions before we introduce the optimizations in the next section. The dataflow in the algorithm is depicted in Fig. 1. The wrapped phase distribution, ϕ (r), is the input data. It is used to calculate the phase estimate ψ est (r). The phase correcting function Q(r) is computed as the difference between the phase estimate and the input phase scaled by 2π and rounded to the nearest integer numbers, i.e., it is an integer map of 2π phase wraps. The phase wrap map is multiplied by 2π and added to the input phase to yield the output unwrapped phase distribution ψ (r).
If the spatial phase distribution ϕ (r) in the acquired data has phase wraps, a correcting function Q(r) needs to be identified to obtain the unwrapped phase distribution ψ (r): where r denotes the spatial coordinates, and Q(r) is a function used to correct the wrapped phase using an integer number of 2π radians. The estimate of Q(r) can be calculated as, Q′(r) is a function of real numbers (can be fractional, positive, and negative). Accordingly, ψ est (r) is the unwrapped phase estimate which should to be identified using only wrapped phase information. This can be done by integrating the derivatives of the phase, and in practice, by applying the Laplace operator ∇ 2 and its inverse ∇ -2 : Schofield proposed to solve this equation for ∇ 2 ψ est (r) by defining a function P(r) such that where j is the imaginary unit. This function has a key property in that it has the same values for wrapped and unwrapped phase, i.e., it is insensitive to phase wraps. The Laplacian of the P(r) function can be expressed as, which can be used to find ∇ 2 ψ est (r): In the above formula, the Laplacian of the unwrapped phase distribution is expressed with the use of function P(r) which can be obtained from the known (measured) wrapped phase distribution. In the original approach, Schofield et al. [43] used Eq. 4 and Euler's formula to derive the Laplacian of the unwrapped phase, Application of the inverse Laplace operator to both sides of this equation leads to an expression for the unwrapped phase estimate, The Laplace operators of an arbitrary function g(r) can be calculated using Fourier transformations (FT) [54], The above version of Schofield's algorithm (Eq. (10)) requires eight Fourier transformations, and constitutes a computational challenge if the phase distribution is multidimensional. A straightforward improvement in the numerical efficiency arises from the simple rearrangement of Eq. (8), as proposed by Jeught et al. [37]:

Optimized FPU algorithm
We propose the modification of the FPU algorithm based on the application of Laplace operators directly on the function P(r). This enables the calculation of the unwrapped phase estimator ψ est (r) from Eq. (6) with complex signals (Eq. (4)) as the input data. The Laplace operators can be computed again using Fourier transformations, The clear advantage of our approach is a reduction of the number of required Fourier transformations to only four (4FT).

Computational efficiency of FPU algorithms
The proposed calculation scheme can be used for phase unwrapping in any phase-sensitive imaging or detection technique. Equations. 10, 12, and 14, can be applied to solve multidimensional phase unwrapping problems. Additionally, the 3D phase distribution can be unwrapped either by the use of the 3D version of the algorithm, or the use of the 2D version applied to each of the 2D data slices. This property makes the approach attractive for Doppler OCT, whereby 3D data is acquired as a sequence of 2D images. Its advantage in Doppler optical coherence tomography is attributed to the readily available complex function P(r), as it will be shown in the next sections. The 2D version requires the computation of twodimensional Fourier transformations (FT), whereas in the 3D option, the three-dimensional FT is calculated. In this work, we have implemented two-dimensional (2D) and threedimensional (3D) versions of the FPU method. To make the implementation as clear as possible (see Code 1 in Ref. [55]) we decided to use MATLAB but without any optimizations. Therefore, all the Fourier transforms are implemented as complex Fast Fourier Transforms (FFTs). Although full optimization of the algorithms is beyond the scope of this article, we can make one remark on this topic. The previously reported versions of the FPU algorithm [30,32] calculated Laplacians on real-value input data. This allowed the calculation of Fourier transformations that led to Hermitian functions. This was followed by the multiplication with a real-valued, symmetrical coordinate function. The latter operation does not alter the phase of the transform, and the subsequent inverse Fourier transformation transforms the Hermitian function to a real-valued function. Therefore, each of the Laplacian calculations can be performed with the use of two real Fourier transformations that in principle require half of the operations needed by the complex Fourier transformation [56]. In the previous versions of the algorithm eight or six real Fourier transformations are required. In contrast, in our approach only one of the Laplacians is applied to real-valued function, while the second one is applied to complex-valued function and is calculated with the use of complex Fourier transformation. In Table 1 we compare our proposed 4FT FPU method with previously published approaches.
If we assume that real the FT requires half the calculation time compared to the complex FT of the same size, one can immediately notice that our approach requires exactly the same number of calculations as the version with six real FTs. It has to be noted, however, that in practical applications, the code that implements the real FT requires additional steps to obtain the results, and is in general less efficient than the complex FT. This is true in principle for most data sizes and available FFT libraries [57]. Therefore, it is in general better to use complex FTs instead of real FTs. Accordingly, we believe that our approach has potential in outperforming the previous implementations.
Experiments performed on the same simulated data sets with all three implementations of the FPU algorithms (using Code 1 run on MATLAB R2016b) showed that the obtained results are practically identical. The number of voxels with different values after phase unwrapping was in the range of 0.01% of the total number of voxels in tested data sets. All the results presented in the following sections were obtained with the 4FT version of the algorithm. In the case of the 3D version of the algorithm, the computational times for the same dataset (256×256×256 voxels) are 1.07±0.03 s, 1.77±0.02 s, and 2.30±0.03 s, while for the 2D version of the algorithm (applied to all the B-scans) the computation times are 0.64±0.06 s, 0.93±0.15 s, and 1.20±0.26 s, for 4FT, 6FT, and 8FT, respectively. This shows a linear dependence of the calculation time on the number of Fourier transforms for the 2D and 3D versions of the algorithms.

Doppler OCT
In the Fourier-domain OCT, an image line (A-scan) is obtained via the Fourier transformation of the interference signal collected in the wavelength space and linearized in the wavenumber space. 2D or 3D datasets are obtained by the lateral scanning of the object. Most frequently, a raster scan is performed in which the beam of light is scanned along the x-and y-axes of the Cartesian coordinates system. Fourier transformed FdOCT data can be represented at any given time t and for any spatial position r, using a complex amplitude A, where ϕ (r, t) is the phase of the interference pattern, and |A(r, t)| 2 is proportional to the light intensity backscattered at position r in the object at time t. To perform Doppler OCT analyses, a series of A-scans needs to be collected from the same spatial location in the object. In practice, the beam position changes are not greater than half of the beam's spot diameter [26,29,30,58]. The phase differences ϕ D between the pairs of A-scans are proportional to the axial component v z of the velocity vector v, where k z is the wavenumber, and n is the index of refraction of the moving medium. The phase wrapping occurs when the time lapse Δt between the acquisition of the A-scan pairs involved in the Doppler OCT analysis exceeds the limit where v z max is the maximum measurable axial velocity at Δt max . The measurements of the velocities exceeding the v z max value require the implementation of phase unwrapping algorithms. Given the representation of the FdOCT signal (Eq. (15)) and the representation of the phase differences in Doppler OCT (Eq. (16)), the function P(r) used to estimate the unwrapped phase in the FPU algorithm (Eq. (14)) can be defined as, The complex amplitudes A are obtained directly from the Fourier transformation of the interference spectra acquired in the FdOCT imaging. Our modified FPU algorithm is therefore naturally suited for the unwrapping of the phase differences in Doppler OCT.

Numerical simulations
To perform the tests of the FPU method and generate reference data for Doppler OCT imaging, a numerical model of the FdOCT signal corresponding to the flow of a scattering medium through a cylindrical tube was created. A laminar flow through a vessel with a circular lumen has a parabolic velocity distribution, changing from maximum in the center of the pipe to zero at the capillary walls according to where |r| is the radial distance from the center of the tube, R is the radius of the tube lumen, and v max is the maximum axial velocity in the center of the tube.
In the simulation, the phase differences (ϕ D ) were computed from the spatial velocity distribution with a given v max according to Eq. (17). The noise was simulated as a complex amplitude that changed randomly across the spatial and temporal coordinates according to A n (r, t) = A nRe (r, t) + jA nIm (r, t) with a normal (N) probability distribution (p) of the real (A nRe ) and imaginary (A nIm ) parts, where p(A nRe ) = N (µ, σ), and p(A nIm ) = N (µ, σ). The mean of the normal distribution is µ = 0 and the standard deviation σ sets the noise level. Pairs of complex OCT amplitudes at each location r separated in time Δt were simulated as sums of the complex signal A s and complex noise A n according to These signals where used to compute P (r, Δt) according to Eq. (18) and applied in the fast phase unwrapping algorithm.

Measures of error
In our experiments we have analyzed the performance of the FPU algorithms depending on the noise level in the input data. To provide a measure of the noise, we use a circular standard deviation of the phase distribution σ circ [59] according to, where the sum spans the spatial coordinates r in the area S where σ circ is calculated, and p(ϕ) is a probability density function which we have approximated with the histogram of the phase distribution. The key property of σ circ is its independence with respect to the phase wraps. Therefore, it can be used as a measure of phase variability of the wrapped input data. Care must be taken when choosing the area S to ensure that the variability of the phase is not dominated by phase gradients. In other words, Eq. (22) must be calculated within a range in which the velocity does not vary. Since the value of σ circ goes to infinity when the phase distribution approaches a uniform distribution, this parameter is sensitive to low-SNR conditions, where none of the methods will correctly unwrap the phase.
We have used the root-mean-square error (RMSE) as an error metric of the FPU algorithm outcomes, which measures the deviation of the unwrapped phase ψ (r) from the reference phase φ D ref distribution for a series of simulated data at increasing noise levels, as shown in Fig. 2.
In the simulations, the reference phase ϕ D ref (r) is identical to the arbitrarily given, noisefree and phase-wrap-free ϕ D (r). However, in the data obtained from the Doppler OCT imaging experiments (ϕ exp ), the true phase distribution that is undisturbed by wraps introduced during the OCT detection, is not known. In other words, RMSE(ψ exp ) cannot be calculated because of the lack of a reliable reference phase ϕ D ref (r). Therefore, the scatterplot of RMSE (ψ exp ) vs. σ circ (ϕ exp ) cannot be generated. However, the circular standard deviation σ circ (ϕ exp ) can still be calculated from the experimental data and referred to the simulated RMSE scatterplot to predict the outcome of the FPU algorithms for a particular noise level given by the same value of σ circ (ϕ). Such analysis can provide insights into the possibility of error-free phase unwrapping to enable the selection of the FPU algorithms that should be used for the output data that are least affected by errors.

Spectral-domain OCT experimental setup
Experiments were performed with a spectral domain OCT setup developed in our laboratory [60]. The light was emitted by a superluminescent diode (λ c = 860 nm, Δλ = 135 nm, Broadlighter T860, Superlum) providing an axial imaging resolution in tissue of ~3 µm. The imaging system was designed to provide a lateral resolution 17 µm in the experiments with the flow phantom and 7.3 µm in eye imaging. The repetition times of the CMOS camera used in the spectrometer were adjusted to obtain phase wrapping in the acquired data. The flow phantom contained 0.5 % fat milk, which was pumped through a glass capillary (600 µm lumen) with a syringe pump (neMESYS 290N, CETONI GmbH). The in vivo imaging of the human eye was performed in a healthy volunteer (I.G). Informed consent was obtained prior to the imaging in compliance with the Declaration of Helsinki [61]. The details of the imaging protocols are listed in Table 2. For the Doppler OCT computations we have used joint spectral and time domain OCT [5] in the fast flow detection regime, i.e., the flow velocity was calculated among groups of adjacent A-scans. To calculate flow based on the imaging of the human eye we have used groups of four A-scans. In the imaging of the flow phantom we have used groups of seven Ascans.

Simulation of Doppler OCT in a cylindrical tube
We have simulated 3D Doppler OCT data of laminar flow in a cylindrical tube with the procedure described in Section 2.4. The time interval Δt max and flow velocity were chosen to generate one phase wrap. One hundred different values of standard deviations σ of the normal probability distribution N (µ=0, σ) were used to simulate increasing noise levels. One hundred datasets were generated at each noise level. The 2D and 3D 4FT FPU algorithms were executed in all the simulated datasets. The RMSE values of the unwrapped phase distributions ψ (r) were calculated and plotted against the circular standard deviation σ circ which was computed from the phase-wrapped input data ϕ (r), as shown in Fig. 2(a).
These metric pairs were calculated in the ROI located in the center of the simulated cylindrical tube. The size of the ROI was selected to include the largest possible flow area with no phase wraps. Herein, we have arbitrarily selected a rectangular window and manually placed it in the center of the tube. However, other shapes of ROIs can be considered or even automatically determined. Regardless of the method of ROI selection, it should to be kept in mind that it is used to calculate the local standard deviation of the phase distribution characteristic of the blood flow in the analyzed vessel. If the ROI area is too small, then the σ circ values will exhibit increased dispersion due to noise. If the ROI area is too large, the phase distribution will no longer be characteristic of the local flow. Instead, it may attain all the values within the (-π, +π) range and rapidly increase. In the limit case of pure phase noise, when the phase distribution in the ROI includes with equal probability values form (-π, +π) range the value of σ circ will be infinite.
For the phase-wrap-free input data (the simulated phase distribution with increasing noise but with no phase wraps), the dependence between the RMSE and σ circ is linear (blue line in Fig. 2(a)). Their values increase proportionally at increasing noise levels. This idealized case serves as a reference for the analysis of more realistic scenarios in which the existing phase wraps have to be removed in the presence of noise. For low-noise levels (small σ circ values), the 2D and 3D versions of the 4FT FPU method yield correct unwrapped phase distributions. Example images obtained with σ circ = 0.69-1.04 are shown in Fig. 2(b). The plots of RMSE (ψ) vs. σ circ (ϕ) follow the reference line, i.e., the RMSE (ψ) is only affected by noise. At increasing noise levels, phase unwrapping errors begin to emerge in the case of the 2D 4F FPU algorithm. Patches of incorrect phase values appear within and outside the images of the tube. Their numbers and total areas increase at increasing σ circ values. As a consequence, the RMSE (ψ) vs. σ circ (ϕ) plot deviates considerably from the reference line. The 3D 4FT FPU method produces phase unwrapping errors at considerably higher noise levels. The RMSE (ψ) vs. σ circ (ϕ) plot deviates from the reference line at higher σ circ values. The 3D version of the algorithm is therefore more immune to noise compared to the 2D variant. It enables reliable phase unwrapping at higher noise levels. The standard deviation of the RMSE values, which is a measure of the repeatability of the outcomes of the algorithm, is considerably smaller in the case of the 3D 4FT FPU (dashed lines in Fig. 2(a)). We postulate that this is caused by the additional number of data from the third dimension which allows more precise derivatives calculations in the presence of noise.
It has to be noted that we did not threshold the Doppler OCT images (based on intensity) as it is commonly done to reject the phase noise from the areas with no OCT signal, neither in the simulated nor in the experimental data. This was because our intention was to demonstrate the influence of the FPU methods on all the components of the Doppler OCT image, including the detected flow, and the stationary parts of the object and noise (where there was no signal from the object). Similar simulations were performed to explore the performance of the algorithms in the case of multiple wraps in the simulated Doppler OCT data. The results are presented in Fig. 3(a-b), with an example of unwrapped phase image presented in Fig. 3(c). It can be observed that an increasing number of wrapping errors appear following unwrapping at lower σ circ values compared to the case of single wraps, as shown in Fig. 2(a). Once again, the 3D 4FT FPU algorithm performs better at lower σ circ values compared to 2D 4FT FPU. i.e., RMSE (ψ) is only affected by the noise level. The deviation of the 2D and 3D FPU RMSE plots from this reference line indicates increasing probability of phase unwrapping errors in ψ (r). As the number of phase wraps increases, the sensitivities of the algorithms to noise also increase. (c) An example of phase-wrapped input data φ (r) for σ circ . = 0.48 rad. 2D: Phase unwrapped using 2D 4FT FPU, and 3D: phase unwrapped using 3D 4FT FPU. Black rectangle denote the region-of-interest used to calculate σ circ .

Doppler OCT imaging of a flow phantom
The OCT experiments with a flow phantom were performed with the imaging parameters listed in Table 2 with Δt max equal to 10 μs. The value of v zmax was selected to provide one phase wrap in the Doppler OCT data. In Fig. 4(b), we present OCT and Doppler OCT images of the glass tube with 0.5 % fat milk as the circulating fluid. The tube was tilted along the Zaxis to ensure a nonzero axial component of the flow velocity vector. As a result, the image appears at increasing depths in subsequent B-scans. Given that the imaging sensitivity of the spectral-domain OCT system exhibits an inherent decay as a function of depth, the signalto-noise ratio decreases along the Z-axis. We have exploited this property in the analyses of the performances of the 2D and 3D 4FT FPU algorithms. Phase unwrapping was performed in the images of the glass tube cross-sections at increasing depth positions, i.e., at decreasing SNR values. Consequently, the circular standard deviation σ circ , which we calculated in the center of the tube, increases in subsequent images, as illustrated by the yellow rectangles in Fig. 4(b). The results are similar to the simulation outcomes. Phase unwrapping errors become evident as σ circ increases, first in the 2D version of the algorithm and in the 3D variant, at larger σ circ . Thus, the experiment confirms higher immunity to noise in the case of the 3D 4FT FPU algorithm. The plot of RMSE (ψ exp ) vs. σ circ (ϕ exp ) cannot be constructed using the experimental data because the phase distribution that is undisturbed by the wraps that were introduced by the OCT detection is not known. Therefore, RMSE (ψ exp ) cannot be reliably calculated. However, we can relate the σ circ (ϕ exp ) values obtained in the experiments with the use of the flow phantom to the circular standard deviations, and to the RMSE values estimated based on simulations. Fig. 4(a) shows increasing σ circ (ϕ exp ) values in subsequent B-scans. The ROI position was approximated by a line that goes form the center of the capillary on the first B-scan in the data set to the last center on the last one. Because a single dataset was used, the dispersion of the data points is high. Therefore, we fitted the linear function σ′ circ (ϕ exp ) to the data points and related in this way the mean values of the circular standard deviation to the simulation results.
Both algorithms yield phase unwrapping errors for lower values of σ′ circ (ϕ exp ) compared to simulations. This may be caused by the simplified noise model used in the simulations. This model does not take into account the few effects present in the experimental data. For example, the reductions of the SNR value owing to the interference signal washout caused by flow, or owing to depth, result in σ circ variations across each B-scan. Nevertheless, the results show that the calculations of σ′ circ (ϕ exp ) in the wrapped Doppler OCT image allow us to predict whether a given version of the algorithm can efficiently unwrap the phase.
To validate whether the 4FT FPU method is able to unwrap multiple phase wraps from experimental data, we have performed OCT experiments, see Fig. 5(a). The imaging parameters listed in Table 2 with Δt max being equal to 20 μs. Herein, the value of v zmax was set to give maximal phase value of approximately 11 rad to generate two phase wraps in the Doppler OCT data as visible in Fig. 5(b). The σ circ calculated from the center of the capillary was equal to 0.45. The results of phase unwrapping are presented in Fig. 5(c,d), whereby the 2D 4FT FPU algorithm results in an increased number of erroneously unwrapped voxels compared to the 3D 4FT FPU algorithm. This result is in agreement with our predictions based on the analysis of the RMSE (ψ exp ) vs. σ circ (ϕ exp ) plots presented in Fig. 3(a, b).

Doppler OCT imaging of the human eye in vivo
The imaging of the human volunteer eye in vivo was performed in regions within the vicinity of the optic nerve head where a) phase wrapping is most likely to occur owing to the presence of large vessels, and b) fast blood flows and increased inclinations relative to the scanning beam of light occur. Two example datasets were acquired at the same location of the eye, as shown in Fig. 6. In Fig. 6(e-h), the eye was intentionally slightly misaligned as compared to the data in Fig. 6(a-d) to decrease the overall imaging sensitivity. The brightness of the structural image is weaker, especially in the left part of the B-scans. The two datasets were acquired at approximately the same location in the eye. The SNR values calculated from the mean signal amplitude and noise variance for the data within the yellow rectangles in Fig. 6(a) and Fig. 6(e) are equal to 22 dB and 18 dB, respectively. Two large vessels are visible, as indicated by arrows A, C (first vessel) and B, D (second vessel), both with phase wrapping in the Doppler OCT image. A smaller vessel was also revealed as indicated by arrows E and F. In the second vessel, the noise was low (σ circ = 0.54 and 0.64), and the phase was unwrapped with no apparent errors associated with the applications of the 2D and 3D 4FT FPU algorithms. The analysis of the first vessel is more interesting. In Fig. 6(a) the signal strength was high in this area. Therefore, the circular standard deviation of the phase was low (σ circ = 0.88). Both variants of the algorithm unwrapped the datasets with no apparent errors. However, in Fig. 6(e), the noise was high in this area (σ circ = 1.84). The 2D 4FT FPU has not only failed but also generated considerable artifacts. The 3D version has only partially removed the phase wraps, thus leaving patches of unwrapped phase in the areas where it has failed. The flow in the small vessel F was correctly unwrapped only by the 3D 4FT FPU algorithm in Fig. 6(d), but was not unwrapped in vessel E by the 2D 4FT FPU algorithm, as observed in Fig. 6(c). Finally, the noise areas in the unwrapped images are not affected by the 3D algorithm, but they show clusters of correlated phase values in the results of the 2D 4FT FPU method. These findings are consistent with the simulations and experiments in the flow phantom.

Discussion and conclusions
The fast phase unwrapping algorithm is well suited for applications in Doppler OCT. Its main advantage is the direct access to the complex signal in the cases of Fourier domain OCT techniques. The modification of the FPU algorithm (implemented in the form of the 4FT FPU method) led to the reduction of the number of Fourier transformations from eight -as reported in the original algorithm -to four. The decreased computational burden accelerated the data processing required to generate the phase unwrapped Doppler OCT images. The availability of volumetric data enabled the extension of the algorithm from two to three dimensions. We have demonstrated that the inclusion of additional information on the phase distribution improved the performance of the 4FT FPU algorithm when it was applied to data with high noise levels. The tests performed on simulated Doppler OCT data obtained from a flow phantom and from the imaging of the human eye in vivo indicated that in lownoise data, 2D and 3D 4FT FPU algorithms provided reliable phase unwrapping outcomes. However, at increased noise levels, phase unwrapping errors occurred earlier in the 2D compared to the 3D versions of the algorithm. Additionally, the 3D 4FT FPU method was more immune to noise. However, it required more data processing. Correspondingly, increased computational computer powers were needed for its implementation. The choice between the 2D and 3D versions of the algorithm should thus be dictated by the noise levels of the Doppler OCT data. The measurement of the noise level in the data with phase-wraps required a metric insensitive to these phase ambiguities. We have used a circular standard deviation of phase distribution (σ circ ) as a measure of noise that was insensitive to phasewraps. We have related this metric to the RMSE which was calculated with the use of the unwrapped simulated data. The resultant plots served as reference plots and indicated which σ circ values of the 2D and 3D algorithms were suitable for the phase-sensitive imaging or for the measurements. For σ circ values larger than the noise level of the 2D algorithm, the 3D 4FT FPU method should be used. For noise levels higher than those associated with the 3D algorithm, the 4FT FPU algorithm cannot provide reliable phase unwrapping. The 4FT FPU algorithm applicability plot generated with the simulated data can thus be used to predict the outcomes of the 4FT FPU algorithm based on the data obtained from the experiments.
The analyzed algorithms unwrapped the phase distribution recorded in the data regardless of the origin of the phase wraps. They only required the application of Laplacian operators and adequately high SNR values. As a result, they have been shown to be insensitive to the origins of wraps, and to their spatial and temporal distributions. Phase wraps may appear in the vessels temporarily or locally in the 3D datasets. The temporal occurrence of phase wraps was connected to natural or induced changes in the blood circulation. Accordingly, natural changes can be caused, for example, by the pulsatile blood flow. Changes may also be induced as a part of diagnostic or therapeutic procedures, or as a result of onset of pathology. The local phase wraps may occur owing to vascular tortuosity. Steeper parts of the vessels may exhibit phase wraps, while parts of the same vessels with shallower inclinations relative to the scanning beam of light may be free from phase wraps. Phase can be unwrapped by the 3D phase unwrapping algorithm as long as the phase distribution is correlated in neighboring B-scans. In other words, phase can be unwrapped if the phase distribution does not change abruptly in successive B-scans due to sparse scanning and/or temporal changes in the blood flow when they slower than the acquisition speeds of the B-scans. Typically, scanning in 3D Doppler OCT is performed to ensure B-scan correlations. Additionally, during the typical acquisition of the 3D datasets, ~2-5 cardiac cycles were recorded in more than 300 B-scans. This ensured adequate temporal sampling for the changes in the blood flow along the vessels.
The Doppler OCT images presented in our study were not intensity thresholded to demonstrate the effects of the phase unwrapping in all the areas of interest, including the vessels, stationary tissue, and areas with no OCT signal (noise only). The images can be intensity thresholded to suppress the phase noise. However, thresholding is typically performed subjectively, and care needs to be taken to avoid loss of relevant information. More importantly, intensity thresholding and the resulting noise removal from selected areas of the images prior to the application of the phase unwrapping algorithms should be avoided. Firstly, the thresholding is subjective and may cause loss of information relevant to successful phase unwrapping. Secondly, phase values equal to zero are arbitrarily introduced to the dataset. Both may result in incorrect phase unwrapping and generation of artefacts.
The current overwhelming popularity of the OCT angiography slowed down the advancement of the Doppler OCT technique. However, the possibility of the quantitative analysis of blood flow should not be overlooked in the development of biomedical diagnostic techniques. Limitations often associated with Doppler OCT can be easily bypassed or eliminated by the application of appropriate experimental and other data analyses methods. In this study, we have introduced a phase unwrapping method which eliminates the phase ambiguity, thus extending the measurable velocity range in the Doppler OCT technique. In addition, we have provided an applicability plot that provides guidance for the use of our method. The 4FT FPU algorithm takes advantage of the specific OCT data acquisition, but can be applied to any phase-sensitive imaging or other measurement technique.