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

Ewelina Pijewska; Iwona Gorczynska; Maciej Szkulmowski

doi:10.1364/BOE.10.001365

1. Introduction

Optical coherence tomography (OCT) is a micrometer-scale imaging technique [1,2] most commonly used in medical diagnostics to visualize structures and detect the functions of biological tissues and cells. OCT relies on the principle of white light interferometry. Analysis of the interference signals originating from backscattered light within the tissue and backreflected from the reference mirror of an interferometer provides information about the spatial distribution of scattering structures within the object. The broader the spectrum of the light is, the higher the axial imaging resolution is. In the case of the Fourier-domain variant of the OCT technique (FdOCT), the interference signal is recorded in the wavelength space. Interference spectra are acquired, digitized, and numerically processed (calibrated in the wavenumber domain, corrected for dispersion, spectrally shaped, and Fourier transformed) to generate complex OCT signals that contain information about the object. One of the advantages of the inherent access to the amplitude and phase of the interference signals is the possibility to devise methods for motion or flow detection and measurement. The motion/flow detection methods are typically used in biomedical OCT imaging for the visualization of the circulatory systems within the living tissues. Their most common application is the use of blood motion to generate contrast for the visualization of vessels using a family of techniques known as optical coherence angiography (OCA or OCTA) [3–5]. However, motion detection techniques can be also augmented with data analyses to provide quantitative information about at least one component of the flow velocity vector [6–11]. Most frequently, the flow measurement is performed with Doppler OCT techniques that either use the Doppler effect directly [12–16], or analyze the phase differences among complex OCT signals [17–21]. Quantitative motion detection and the Doppler OCT methods are often jointly referred to as OCT velocimetry.

Doppler OCT methods can provide information similar to OCTA techniques, i.e., about the 3D architecture of active vascular systems [13,22–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–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–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–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–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, i.e., 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–53].

In this study, we have adapted the fast phase unwrapping (FPU) algorithm for spectral-domain 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.

2. Methods

2.1 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.

Fig. 1 Graphical representation of the data flow in the fast phase unwrapping (FPU) algorithm. 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).

Download Full Size | PDF

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):

ψ (r) = φ (r) + 2 π \cdot Q (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) = {(2 π)}^{- 1} (ψ_{e s t} (r) - φ (r)) .

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 ∇² and its inverse ∇⁻²:

Q^{'} (r) = {(2 π)}^{- 1} \nabla^{- 2} (\nabla^{2} ψ_{e s t} (r) - \nabla^{2} φ (r)) .

Schofield proposed to solve this equation for ∇²ψ_est (r) by defining a function P(r) such that

P (r) = \exp (j ψ_{e s t} (r)) = \exp (j φ (r)),

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,

\nabla^{2} P (r) = \nabla (P (r) j \nabla ψ_{e s t} (r)) = - P (r) {(\nabla ψ_{e s t} (r))}^{2} + j P (r) \nabla^{2} ψ_{e s t} (r),

which can be used to find ∇²ψ_est (r):

\nabla^{2} ψ_{e s t} (r) = Im {P {(r)}^{- 1} \nabla^{2} P (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,

\begin{array}{l} \nabla^{2} ψ_{e s t} (r) = Im {\exp (- j φ (r)) \nabla^{2} \exp (j φ (r))} = \\ = \cos φ (r) \nabla^{2} \sin φ (r) - \sin φ (r) \nabla^{2} \cos φ (r) . \end{array}

Application of the inverse Laplace operator to both sides of this equation leads to an expression for the unwrapped phase estimate,

ψ_{e s t} (r) = \nabla^{- 2} [\cos φ (r) \nabla^{2} \sin φ (r)] - \nabla^{- 2} [\sin φ (r) \nabla^{2} \cos φ (r)] .

The Laplace operators of an arbitrary function g(r) can be calculated using Fourier transformations (FT) [54],

\begin{array}{l} \nabla^{2} g (r) = - {(2 π)}^{2} F^{- 1} {K^{2} F {g (r)}} \\ \nabla^{- 2} g (r) = - {(2 π)}^{- 2} F^{- 1} {K^{- 2} F {g (r)}}, \end{array}

where K is Fourier space conjugate of the vector r, and

F {}

and

F^{- 1} {}

are forward and inverse Fourier transformation operators. With this observation, Eq. (8) can be rewritten as,

\begin{matrix} ψ_{e s t} (r) = F^{- 1} {K^{- 2} F {\cos φ (r) \cdot F^{- 1} {K^{2} F {\sin φ (r)}}}} + \\ - F^{- 1} {K^{- 2} F {\sin φ (r) \cdot F^{- 1} {K^{2} F {\cos φ (r)}}}} . \end{matrix}

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]:

ψ_{e s t} (r) = \nabla^{- 2} [\cos φ (r) \nabla^{2} \sin φ (r) - \sin φ (r) \nabla^{2} \cos φ (r)],

which leads to a simplified version of Eq. (10) with six (instead of eight) Fourier transformations,

ψ_{e s t} (r) = F^{- 1} {K^{- 2} F {\cos φ (r) \cdot F^{- 1} {K^{2} F {\sin φ (r)}} - \sin φ (r) \cdot F^{- 1} {K^{2} F {\cos φ (r)}}}} .

If the periodic boundary conditions required for the existence of the Fourier transforms are not fulfilled, phase unwrapping errors arise, thus preventing direct use of the estimate ψ_est (r) as the final unwrapped phase distribution. To overcome this issue, Schofield proposed the use of the phase correcting function defined as a difference between unwrapped and wrapped phase distributions, scaled by 2π, and rounded to the nearest integer values,

Q (r) = round ({(2 π)}^{- 1} (ψ_{e s t} (r) - φ (r))) = round (Q^{'} (r)) .

Q(r) can be thought of as an integer map of 2π phase wraps. The phase correcting function was subsequently used in Eq. (1) to calculate the unwrapped phase.

2.2 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,

ψ_{e s t} (r) = F^{- 1} {K^{- 2} F {Im {P {(r)}^{- 1} \cdot F^{- 1} {K^{2} F {P (r)}}}}} .

The clear advantage of our approach is a reduction of the number of required Fourier transformations to only four (4FT).

2.3 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 two-dimensional Fourier transformations (FT), whereas in the 3D option, the three-dimensional FT is calculated. In this work, we have implemented two-dimensional (2D) and three-dimensional (3D) versions of the FPU method. To make the implementation as clear as possible (see Code 1 [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.

Table 1. Comparison of fast phase unwrapping (FPU) algorithms

View Table | View all tables in this article

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 [55] 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 data set (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.

2.4 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 data sets 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,

A (r, t) = | A (r, t) | \exp (j φ (r, t)),

where φ (r, t) is the phase of the interference pattern, and |A(r, t)|² 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,

φ_{D} (r, Δ t) = Im {\ln (A (r, t) A^{*} (r, t + Δ t))} = 2 k_{z} n v_{z} (r) Δ t,

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

Δ t_{\max} = \frac{π}{2 k_{z} n v_{z \max}},

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,

P (r, Δ t) = \frac{A (r, t) A^{*} (r, t + Δ t)}{| A (r, t) A^{*} (r, t + Δ t) |} = \exp (j φ_{D} (r, Δ t)) .

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.

2.5 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

v_{z} (r) = {\begin{matrix} v_{\max} - (v_{\max} / R^{2}) {| r |}^{2} & for | r | < R & (inside the capillary) \\ 0 & for | r | \geq R & (outside the capillary) \end{matrix}

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

A (r, t) = A_{s} (r, t) + A_{n} (r, t) = | A (r, t) | \exp [i φ (r, t)] + A_{n} (r, t),

\begin{array}{l} A (r, t + Δ t) = A_{s} (r, t + Δ t) + A_{n} (r, t + Δ t) = \\ = | A (r, t) | \exp [i φ (r, t) + i φ_{D} (r, Δ t)] + A_{n} (r, t + Δ t) . \end{array}

These signals where used to compute P (r, Δt) according to Eq. (18) and applied in the fast phase unwrapping algorithm.

2.6 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,

σ_{c i r c} (φ (r, t)) = \sqrt{2 \ln (1 / | \sum_{r \in S} p (φ (r, t)) \exp (i φ (r, t)) |)},

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} (r),

R M S E (ψ) = \sqrt{M^{- 1} \sum_{r} {(ψ (r) - φ_{D r e f} (r))}^{2}},

where

M

is the number of pixels in the analyzed region-of-interest (ROI) in the data set. In the absence of phase unwrapping errors, RMSE depends only on the noise level in the simulated data, and increases linearly at increasing circular standard deviations of the phase σ_circ (φ). Errors in phase unwrapping, i.e., in ψ (r), increase the RMSE. This leads to deviations from its expected linear dependence on σ_circ (φ).

To analyze the influence of noise on the outcomes of the FPU algorithms, we have generated a diagram in which we have plotted the RMSE values of the unwrapped phase ψ (r) as a function of the σ_circ values computed from the corresponding wrapped phase φ (r) distribution for a series of simulated data at increasing noise levels, as shown in Fig. 2.

Fig. 2 Result from the applications of 2D and 3D fast phase unwrapping in a set of simulated 3D Doppler OCT data at increasing noise levels. (a) A plot of phase unwrapping error metrics (RMSE (ψ)) vs. phase noise metrics (σ_circ (φ)). Dots represent mean RMSE (ψ) values computed from data sets comprising 100 data points each. Dashed lines indicate the standard deviations of the mean RMSE [2D FPU algorithm (green) and 3D FPU algorithm (black) outcomes]. The blue plot shows a linear dependence of the RMSE (ψ) on σ_circ (φ) in the absence of phase unwrapping errors in ψ (r), i.e. RMSE (ψ) is only affected by the noise level. The deviations of the 2D and 3D FPU RMSE plots from this reference line indicate an increasing probability of phase unwrapping errors in ψ (r). The inset presents an example of phase-wrapped input data φ (r). Yellow rectangle denote the region-of-interest used to calculate σ_circ. (b) Example cross-sectional images of the unwrapped phase ψ (r) extracted from four selected 3D data sets. The values of the circular standard deviation σ_circ (φ) in the input data are listed in the top left corners. The top row of images shows the results of the 2D FPU method, and the bottom row demonstrates the outcomes of the 3D FPU algorithm. The Doppler OCT images were intentionally not thresholded for noise elimination to demonstrate the influences of the 4FT FPU algorithms on the noise areas. The results of the algorithms for all the data points in the plots on the top are presented in Visualization 1.

Download Full Size | PDF

In the simulations, the reference phase φ_{D ref} (r) is identical to the arbitrarily given, noise-free 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.

2.7 Spectral-domain OCT experimental setup

Experiments were performed with a spectral-domain OCT setupdeveloped 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 of 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.

Table 2. Imaging protocols applied for the experiments with the use of the flow phantom and human eye.

View Table | View all tables in this article

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 A-scans.

3. Results

3.1 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 data sets were generated at each noise level. The 2D and 3D 4FT FPU algorithms were executed in all the simulated data sets. 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.

Fig. 3 Phase unwrapping using the 4FT FPU algorithm for multiple phase wraps in four sets of data at different maximum values of unwrapped phases ψ_max. (a, b) Plots of phase unwrapping error metric (RMSE (ψ)) as a function of phase noise metric (σ_circ (φ)) for 3D Doppler OCT data simulated at increasing noise levels. (a) 2D 4FT FPU and (b) 3D 4FT FPU. Dots represent the mean RMSE (ψ) values computed from data sets comprising 100 data points each. Dashed lines indicate the standard deviations of the mean RMSE. Colors represent different maximum values of unwrapped phases ψ_max. The blue solid line shows a linear dependence of the RMSE (ψ) on σ_circ (φ) in the absence of phase unwrapping errors in ψ (r), 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.

Download Full Size | PDF

3.2 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 Z–axis 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 signal-to-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.

Fig. 4 Phase unwrapping in the flow phantom. (a) Plot of the circular standard deviation values (σ_circ) in subsequent B-scans. The values were calculated based on the regions-of-interest indicated by the yellow rectangles in the wrapped phase images in (b); black line: linear fit to the data points. (b) A₀: cross-sectional structural OCT images of the glass tube with milk as the circulating fluid at four locations in the 3D data set; φ: wrapped phase; 2D: phase unwrapped with the 2D 4FT FPU method, arrows B and C indicate the phase unwrapping artefacts introduced in the areas of noise. 3D: outcomes of the 3D 4FT FPU method, arrow A: error-free phase unwrapping of the flow, arrow D: erroneous phase unwrapping (phase has remained wrapped). The colors in the images represent Doppler OCT phase values (red: blood flow direction against the incoming beam of light, and blue: blood flow along the light propagation direction). Increasing color intensities indicate increasing axial flow velocity values as indicated by the color bar. The thresholding of the Doppler OCT images to remove the noise was intentionally avoided to demonstrate how the 4FT FPU algorithms influenced the noise areas. Yellow rectangle denote the region-of-interest used to calculate σ_circ. A movie showing all the frames from the data set is presented in Visualization 2.

Download Full Size | PDF

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. Figure 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 data set 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).

Fig. 5 Comparison of 2D and 3D 4FT FPU methods in the imaging of the flow phantom with the use of multiple phase wraps. (a) Cross-sectional structural OCT image of the glass tube with milk as the circulating fluid. (b) Doppler OCT image with multiple phase wraps within the capillary lumen. Black rectangle denote the region-of-interest used to calculate σ_circ. (c) Outcome of the 2D 4FT FPU method. (d) Results of the 3D 4FT FPU method.

Download Full Size | PDF

3.3 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 data sets 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 data sets 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.

Fig. 6 Comparison of the 2D and 3D 4FT FPU methods in the imaging of the human retina in vivo in the vicinity of the head of the optic nerve. (e-h) The eye was intentionally misaligned to decrease the imaging sensitivity as compared to (a-d). (a, e) Structural OCT images of the eye fundus. Arrows A, B, C, and D, indicate large and visible vessels. (b, f) Doppler OCT image with phase wraps within the vessels. (c, g) Results following the use of the 2D 4FT FPU method. (d, h) Results following the use of the 3D 4FT FPU method. Arrows E and F point to small vessels which are not apparent in the structural image but are detected with Doppler OCT. Colors in the images represent the Doppler OCT phase values (red: blood flow direction against the incoming beam of light, blue: flow direction along the light propagation direction). Increasing color intensity indicates increasing axial flow velocity values as indicated by the color bar. The thresholding of the Doppler OCT images to remove the noise was intentionally avoided to demonstrate how the 4FT FPU algorithms influenced the noise areas. Yellow rectangles denote the regions-of-interest used to calculate the signal-to-noise ratio and σ_circ. Movies showing all the frames from both data sets are presented in Visualization 3 (top) and Visualization 4 (bottom).

Download Full Size | PDF

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 data sets 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.

4. 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 low-noise 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 phase-wraps. 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 data sets. 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 data sets, ~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 data set. 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.

Funding

Fundacja na rzecz Nauki Polskiej (POIR.04.04.00-00-2070/16-00).

Acknowledgements

We dedicate this work to our mentor Professor Andrzej Kowalczyk, a creator of the OCT School at the Nicolaus Copernicus University in Torun, Poland, who retired last year from his academic duties.

The project “FreezEYE Tracker – ultrafast system for image stabilization in biomedical imaging” was conducted within the TEAM TECH Program of the Foundation for Polish Science, co-financed by the European Union under the European Regional Development Fund.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef] [PubMed]

2. A. F. Fercher, C. K. Hitzenberger, G. Kamp, and S. Y. Elzaiat, “Measurement of Intraocular Distances by Backscattering Spectral Interferometry,” Opt. Commun. 117(1-2), 43–48 (1995). [CrossRef]

3. C. L. Chen and R. K. Wang, “Optical coherence tomography based angiography [Invited],” Biomed. Opt. Express 8(2), 1056–1082 (2017). [CrossRef] [PubMed]

4. S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, “Optical coherence angiography,” Opt. Express 14(17), 7821–7840 (2006). [CrossRef] [PubMed]

5. R. K. Wang, S. L. Jacques, Z. Ma, S. Hurst, S. R. Hanson, and A. Gruber, “Three dimensional optical angiography,” Opt. Express 15(7), 4083–4097 (2007). [CrossRef] [PubMed]

6. J. Tokayer, Y. Jia, A. H. Dhalla, and D. Huang, “Blood flow velocity quantification using split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Biomed. Opt. Express 4(10), 1909–1924 (2013). [CrossRef] [PubMed]

7. N. Uribe-Patarroyo, M. Villiger, and B. E. Bouma, “Quantitative technique for robust and noise-tolerant speed measurements based on speckle decorrelation in optical coherence tomography,” Opt. Express 22(20), 24411–24429 (2014). [CrossRef] [PubMed]

8. V. J. Srinivasan, S. Sakadzić, I. Gorczynska, S. Ruvinskaya, W. Wu, J. G. Fujimoto, and D. A. Boas, “Quantitative cerebral blood flow with optical coherence tomography,” Opt. Express 18(3), 2477–2494 (2010). [CrossRef] [PubMed]

9. V. J. Srinivasan, H. Radhakrishnan, E. H. Lo, E. T. Mandeville, J. Y. Jiang, S. Barry, and A. E. Cable, “OCT methods for capillary velocimetry,” Biomed. Opt. Express 3(3), 612–629 (2012). [CrossRef] [PubMed]

10. W. J. Choi, Y. Li, W. Qin, and R. K. Wang, “Cerebral capillary velocimetry based on temporal OCT speckle contrast,” Biomed. Opt. Express 7(12), 4859–4873 (2016). [CrossRef] [PubMed]

11. F. Jaillon, S. Makita, and Y. Yasuno, “Variable velocity range imaging of the choroid with dual-beam optical coherence angiography,” Opt. Express 20(1), 385–396 (2012). [CrossRef] [PubMed]

12. M. Szkulmowski, A. Szkulmowska, T. Bajraszewski, A. Kowalczyk, and M. Wojtkowski, “Flow velocity estimation using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 16(9), 6008–6025 (2008). [CrossRef] [PubMed]

13. Y. K. Tao, A. M. Davis, and J. A. Izatt, “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Opt. Express 16(16), 12350–12361 (2008). [CrossRef] [PubMed]

14. J. Walther and E. Koch, “Relation of joint spectral and time domain optical coherence tomography (jSTdOCT) and phase-resolved Doppler OCT,” Opt. Express 22(19), 23129–23146 (2014). [CrossRef] [PubMed]

15. A. Bouwens, D. Szlag, M. Szkulmowski, T. Bolmont, M. Wojtkowski, and T. Lasser, “Quantitative lateral and axial flow imaging with optical coherence microscopy and tomography,” Opt. Express 21(15), 17711–17729 (2013). [CrossRef] [PubMed]

16. R. A. Leitgeb, R. M. Werkmeister, C. Blatter, and L. Schmetterer, “Doppler optical coherence tomography,” Prog. Retin. Eye Res. 41, 26–43 (2014). [CrossRef] [PubMed]

17. R. Leitgeb, L. Schmetterer, W. Drexler, A. Fercher, R. Zawadzki, and T. Bajraszewski, “Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,” Opt. Express 11(23), 3116–3121 (2003). [CrossRef] [PubMed]

18. B. White, M. Pierce, N. Nassif, B. Cense, B. Park, G. Tearney, B. Bouma, T. Chen, and J. de Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical coherence tomography,” Opt. Express 11(25), 3490–3497 (2003). [CrossRef] [PubMed]

19. B. Vakoc, S. Yun, J. de Boer, G. Tearney, and B. Bouma, “Phase-resolved optical frequency domain imaging,” Opt. Express 13(14), 5483–5493 (2005). [CrossRef] [PubMed]

20. A. Szkulmowska, M. Szkulmowski, A. Kowalczyk, and M. Wojtkowski, “Phase-resolved Doppler optical coherence tomography--limitations and improvements,” Opt. Lett. 33(13), 1425–1427 (2008). [CrossRef] [PubMed]

21. S. Tozburun, C. Blatter, M. Siddiqui, E. F. J. Meijer, and B. J. Vakoc, “Phase-stable Doppler OCT at 19 MHz using a stretched-pulse mode-locked laser,” Biomed. Opt. Express 9(3), 952–961 (2018). [CrossRef] [PubMed]

22. R. K. Wang and L. An, “Doppler optical micro-angiography for volumetric imaging of vascular perfusion in vivo,” Opt. Express 17(11), 8926–8940 (2009). [CrossRef] [PubMed]

23. A. Szkulmowska, M. Szkulmowski, D. Szlag, A. Kowalczyk, and M. Wojtkowski, “Three-dimensional quantitative imaging of retinal and choroidal blood flow velocity using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express 17(13), 10584–10598 (2009). [CrossRef] [PubMed]

24. I. Grulkowski, I. Gorczynska, M. Szkulmowski, D. Szlag, A. Szkulmowska, R. A. Leitgeb, A. Kowalczyk, and M. Wojtkowski, “Scanning protocols dedicated to smart velocity ranging in spectral OCT,” Opt. Express 17(26), 23736–23754 (2009). [CrossRef] [PubMed]

25. B. J. Vakoc, R. M. Lanning, J. A. Tyrrell, T. P. Padera, L. A. Bartlett, T. Stylianopoulos, L. L. Munn, G. J. Tearney, D. Fukumura, R. K. Jain, and B. E. Bouma, “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. 15(10), 1219–1223 (2009). [CrossRef] [PubMed]

26. V. J. Srinivasan, S. Sakadzić, I. Gorczynska, S. Ruvinskaya, W. Wu, J. G. Fujimoto, and D. A. Boas, “Quantitative cerebral blood flow with optical coherence tomography,” Opt. Express 18(3), 2477–2494 (2010). [CrossRef] [PubMed]

27. A. Bouwens, T. Bolmont, D. Szlag, C. Berclaz, and T. Lasser, “Quantitative cerebral blood flow imaging with extended-focus optical coherence microscopy,” Opt. Lett. 39(1), 37–40 (2014). [CrossRef] [PubMed]

28. S. Yazdanfar, C. Yang, M. Sarunic, and J. Izatt, “Frequency estimation precision in Doppler optical coherence tomography using the Cramer-Rao lower bound,” Opt. Express 13(2), 410–416 (2005). [CrossRef] [PubMed]

29. B. Park, M. C. Pierce, B. Cense, S. H. Yun, M. Mujat, G. Tearney, B. Bouma, and J. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 microm,” Opt. Express 13(11), 3931–3944 (2005). [CrossRef] [PubMed]

30. B. J. Vakoc, G. J. Tearney, and B. E. Bouma, “Statistical properties of phase-decorrelation in phase-resolved Doppler optical coherence tomography,” IEEE Trans. Med. Imaging 28(6), 814–821 (2009). [CrossRef] [PubMed]

31. E. Koch, J. Walther, and M. Cuevas, “Limits of Fourier domain Doppler-OCT at high velocities,” Sens. Actuators A Phys. 156(1), 8–13 (2009). [CrossRef]

32. H. C. Hendargo, R. P. McNabb, A. H. Dhalla, N. Shepherd, and J. A. Izatt, “Doppler velocity detection limitations in spectrometer-based versus swept-source optical coherence tomography,” Biomed. Opt. Express 2(8), 2175–2188 (2011). [CrossRef] [PubMed]

33. A. C. Chan, E. Y. Lam, and V. J. Srinivasan, “Comparison of Kasai autocorrelation and maximum likelihood estimators for Doppler optical coherence tomography,” IEEE Trans. Med. Imaging 32(6), 1033–1042 (2013). [CrossRef] [PubMed]

34. S. Makita, F. Jaillon, I. Jahan, and Y. Yasuno, “Noise statistics of phase-resolved optical coherence tomography imaging: single-and dual-beam-scan Doppler optical coherence tomography,” Opt. Express 22(4), 4830–4848 (2014). [CrossRef] [PubMed]

35. R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: Two-dimensional phase unwrapping,” Radio Sci. 23(4), 713–720 (1988). [CrossRef]

36. S. Chavez, Q. S. Xiang, and L. An, “Understanding phase maps in MRI: a new cutline phase unwrapping method,” IEEE Trans. Med. Imaging 21(8), 966–977 (2002). [CrossRef] [PubMed]

37. V. S. Jeught, J. Sijbers, and J. J. Dirckx, “Fast Fourier-Based Phase Unwrapping on the Graphics Processing Unit in Real-Time Imaging Applications,” J. Imaging 1(1), 31–44 (2015). [CrossRef]

38. D. C. Ghiglia and M. D. Pritt, Two-dimensional phase unwrapping: theory, algorithms, and software (Wiley, New York, 1998).

39. S. H. Yun, G. Tearney, J. de Boer, and B. Bouma, “Motion artifacts in optical coherence tomography with frequency-domain ranging,” Opt. Express 12(13), 2977–2998 (2004). [CrossRef] [PubMed]

40. J. Walther, G. Mueller, H. Morawietz, and E. Koch, “Signal power decrease due to fringe washout as an extension of the limited Doppler flow measurement range in spectral domain optical coherence tomography,” J. Biomed. Opt. 15(4), 041511 (2010). [CrossRef] [PubMed]

41. J. Walther and E. Koch, “Impact of a detector dead time in phase-resolved Doppler analysis using spectral domain optical coherence tomography,” J. Opt. Soc. Am. A 34(2), 241–251 (2017). [CrossRef] [PubMed]

42. K. Itoh, “Analysis of the phase unwrapping algorithm,” Appl. Opt. 21(14), 2470 (1982). [CrossRef] [PubMed]

43. M. A. Schofield and Y. Zhu, “Fast phase unwrapping algorithm for interferometric applications,” Opt. Lett. 28(14), 1194–1196 (2003). [CrossRef] [PubMed]

44. W. Li, A. V. Avram, B. Wu, X. Xiao, and C. Liu, “Integrated Laplacian-based phase unwrapping and background phase removal for quantitative susceptibility mapping,” NMR Biomed. 27(2), 219–227 (2014). [CrossRef] [PubMed]

45. H. Bagher-Ebadian, Q. Jiang, and J. R. Ewing, “A modified Fourier-based phase unwrapping algorithm with an application to MRI venography,” J. Magn. Reson. Imaging 27(3), 649–652 (2008). [CrossRef] [PubMed]

46. M. Loecher, E. Schrauben, K. M. Johnson, and O. Wieben, “Phase unwrapping in 4D MR flow with a 4D single-step laplacian algorithm,” J. Magn. Reson. Imaging 43(4), 833–842 (2016). [CrossRef] [PubMed]

47. E. Barnhill, P. Kennedy, C. L. Johnson, M. Mada, and N. Roberts, “Real-time 4D phase unwrapping applied to magnetic resonance elastography,” Magn. Reson. Med. 73(6), 2321–2331 (2015). [CrossRef] [PubMed]

48. T. Latychevskaia, P. Formanek, C. T. Koch, and A. Lubk, “Off-axis and inline electron holography: Experimental comparison,” Ultramicroscopy 110(5), 472–482 (2010). [CrossRef]

49. D. Parshall and M. K. Kim, “Digital holographic microscopy with dual-wavelength phase unwrapping,” Appl. Opt. 45(3), 451–459 (2006). [CrossRef] [PubMed]

50. H. C. Hendargo, M. Zhao, N. Shepherd, and J. A. Izatt, “Synthetic wavelength based phase unwrapping in spectral domain optical coherence tomography,” Opt. Express 17(7), 5039–5051 (2009). [CrossRef] [PubMed]

51. Y. Wang, A. A. Fawzi, O. Tan, X. Zhang, and D. Huang, “Flicker-induced changes in retinal blood flow assessed by Doppler optical coherence tomography,” Biomed. Opt. Express 2(7), 1852–1860 (2011). [CrossRef] [PubMed]

52. Y. Wang, D. Huang, Y. Su, and X. S. Yao, “Two-dimensional phase unwrapping in Doppler Fourier domain optical coherence tomography,” Opt. Express 24(23), 26129–26145 (2016). [CrossRef] [PubMed]

53. S. Xia, Y. Huang, S. Peng, Y. Wu, and X. Tan, “Robust phase unwrapping for phase images in Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt. 22(3), 36014 (2017). [CrossRef] [PubMed]

54. T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. Orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A 13(8), 1670–1682 (1996). [CrossRef]

55. E. Pijewska, I. Gorczynska, and M. Szkulmowski, “Matlab implementation of various forms of the Fast Phase Unwrapping algorithm for 2 and 3 dimensional data,” figshare (2019) [retrieved 11 Feb 2019], https://doi.org/10.6084/m9.figshare.7308413.

56. H. Sorensen, D. Jones, M. Heideman, and C. Burrus, “Real-valued fast Fourier transform algorithms,” IEEE Trans. Acoust. Speech Signal Process. 35(6), 849–863 (1987). [CrossRef]

57. M. Frigo and S. G. Johnson, “benchFFT” (2019), accessed on 15–01–2019, http://www.fftw.org/benchfft/.

58. N. Uribe-Patarroyo and B. E. Bouma, “Velocity gradients in spatially resolved laser Doppler flowmetry and dynamic light scattering with confocal and coherence gating,” Phys. Rev. E 94(2-1), 022604 (2016). [CrossRef] [PubMed]

59. N. I. Fisher, Statistical Analysis of Circular Data (Cambridge University Press, 1995).

60. D. Ruminski, B. L. Sikorski, D. Bukowska, M. Szkulmowski, K. Krawiec, G. Malukiewicz, L. Bieganowski, and M. Wojtkowski, “OCT angiography by absolute intensity difference applied to normal and diseased human retinas,” Biomed. Opt. Express 6(8), 2738–2754 (2015). [CrossRef] [PubMed]

61. World Medical Association, “World Medical Association Declaration of Helsinki: ethical principles for medical research involving human subjects,” JAMA 310(20), 2191–2194 (2013). [CrossRef] [PubMed]

Name	Description
Code 1	Matlab implementation of various forms of the Fast Phase Unwrapping algorithm for 2 and 3 dimensional data.
Visualization 1	2D and 3D fast phase unwrapping results in a set of simulated 3D Doppler OCT data with increasing noise level. Top: a plot of phase unwrapping error-metric (RMSE (psi)) vs. phase noise-metric (scirc (fi)). Dots represent mean RMSE (psi) values comput
Visualization 2	Phase unwrapping in the flow phantom. Top: plot of the circular standard deviation values (scirc) in subsequent B-scans. The values were calculated in the region of interest indicated by yellow rectangles in wrapped phase images. Blue line – linear f
Visualization 3	Comparison of 2D and 3D complex FPU method in imaging of the human retina in vivo with the eye intentionally misaligned to decrease the imaging sensitivity. Colors in the images represent DOCT phase values: red – blood flow direction against the inco
Visualization 4	Comparison of 2D and 3D complex FPU method in imaging of the human retina in vivo with the eye intentionally misaligned to decrease the imaging sensitivity. Colors in the images represent DOCT phase values: red – blood flow direction against the inco

Algorithm	Number of multidimensional Fourier Transforms (FT)	Formula** for the unwrapped phase distribution ψ_est with a given wrapped phase φ
Original FPU [43] – Eq. (10)	Eight real FTs	$ψ_{e s t} = \nabla^{- 2} (\sin φ \cdot \nabla^{2} \cos φ) - \nabla^{- 2} (\cos φ \cdot \nabla^{2} \sin φ)$
Re-arranged FPU [37] – Eq. (12)	Six real FTs	$ψ_{e s t} = \nabla^{- 2} (\sin φ \cdot \nabla^{2} \cos φ - \cos φ \cdot \nabla^{2} \sin φ)$
4FT FPU – this article Eq. (14)	Two real FTs + two complex FTs	$ψ_{e s t} = \nabla^{- 2} Im {\exp (- j φ) \nabla^{2} \exp (j φ)}$

	Flow phantom	Human eye
A-scans × B-sans	500 × 800	600 × 700
Transverse scan range (x, y)	1.6 × 0.9 mm	0.9 × 0.9 mm
Camera repetition time (Δt_max in Eq. (16)	10 μs (20 μs)	30 μs
v_zmax (Eq. (17)	21.5 mm/s (10.75 mm/s)	7.2 mm/s

Algorithm	Number of multidimensional Fourier Transforms (FT)	Formula** for the unwrapped phase distribution ψ_est with a given wrapped phase φ
Original FPU [43] – Eq. (10)	Eight real FTs	$ψ_{e s t} = \nabla^{- 2} (\sin φ \cdot \nabla^{2} \cos φ) - \nabla^{- 2} (\cos φ \cdot \nabla^{2} \sin φ)$
Re-arranged FPU [37] – Eq. (12)	Six real FTs	$ψ_{e s t} = \nabla^{- 2} (\sin φ \cdot \nabla^{2} \cos φ - \cos φ \cdot \nabla^{2} \sin φ)$
4FT FPU – this article Eq. (14)	Two real FTs + two complex FTs	$ψ_{e s t} = \nabla^{- 2} Im {\exp (- j φ) \nabla^{2} \exp (j φ)}$

	Flow phantom	Human eye
A-scans × B-sans	500 × 800	600 × 700
Transverse scan range (x, y)	1.6 × 0.9 mm	0.9 × 0.9 mm
Camera repetition time (Δt_max in Eq. (16)	10 μs (20 μs)	30 μs
v_zmax (Eq. (17)	21.5 mm/s (10.75 mm/s)	7.2 mm/s

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

Abstract

1. Introduction

2. Methods

2.1 Fast phase unwrapping

2.2 Optimized FPU algorithm

2.3 Computational efficiency of FPU algorithms

2.4 Doppler OCT

2.5 Numerical simulations

2.6 Measures of error

2.7 Spectral-domain OCT experimental setup

3. Results

3.1 Simulation of Doppler OCT in a cylindrical tube

3.2 Doppler OCT imaging of a flow phantom

3.3 Doppler OCT imaging of the human eye in vivo

4. Discussion and conclusions

Funding

Acknowledgements

Disclosures

References

Supplementary Material (5)

Cited By

Figures (6)

Tables (2)

Equations (23)

Biomedical Optics Express