1. Introduction

Many interferometric techniques have been developed for the purpose of quantitative measurement of physical properties in biological specimens based on retrieval of phase information from interference signals [1–5]. For example, phase contrast microscopy [1] can be used to image transparent biological samples by converting phase information to intensity differences qualitatively. Phase shifting interferometry and Fourier phase microscopy [2] can be used to obtain cellular images. Optical coherence tomography (OCT) [6] is a noninvasive technology that is currently used for in vivo high resolution sectional imaging of microstructure in biological tissues [7–11]. Based on OCT technique, nanometer scale phase microscopy [12–14] was developed to detect depth resolved information from a given sample. In addition to morphological structure imaging, the measurement of blood flow in biological tissues can be obtained by Doppler Fourier-domain (FD) OCT [15–23], which shows great medical importance for vascular circulation investigation in human body.

In Doppler FD-OCT [24–28], light reflected by moving blood incurs a Doppler frequency shift Δν, which is proportional to the flow velocity V, Δν = −2nVcosα/λ₀, where n is the refractive index of the medium, α is the angle between the OCT probe beam and blood flow, and λ₀ is the center wavelength of the light source. This frequency shift introduces a phase shift in the detected OCT spectral interference pattern. After fast Fourier transform, the result is a complex function F(z) characterized by amplitude A(z) and phase Φ(z). The phase difference (or the Doppler phase), φ(z) = Φ_{j + 1}(z)- Φ_j(z), between two sequential axial lines (A-lines) at each pixel is calculated to determine Doppler frequency shift (Δν = φ(z)/(2πT)) using the algorithm [28]:

Phase unwrapping is an essential step for data processing in phase related interferometric field. In the noise-free case, the unwrapping process can be achieved through searching jump in the wrapped phase signal, and adding an appropriate integer multiple of 2π to each data pixel, or accumulating the phase difference (phase gradient) from one point to the other [29–32]. However, real experiment data always contain noise. Defects in the interference signals, such as noise, shadow and phase discontinuity, are the main difficulties in the phase retrieval. To remove these invalid data points and recover precise phase results can be a complex and time-consuming process.

Wrapped phase signal can be unwrapped correctly after a proper noise filtering, but such a noise filter is not easy to be obtained. To reduce noise before unwrapping, several methods have been developed based on the path following algorithm [33–38], such as the multilook filter method [33], the Fourier transform-based method, and the local statistic filter which divides the phase data into small blocks to allow the signals of each block to be filtered separately [34, 35]. Based on windowed Fourier transform, algorithms are proposed to reduce noise through processing signals locally [37, 38]. However, a proper window size is critical for achieving good result. Moreover, the choice for the best path is not always obvious, especially when the noise statistics are unknown or in the case of object-dependent. The other type of phase unwrapping algorithm is based on a least-square approach [39–43]. These methods use a principle that minimizes the difference between the partial derivatives of wrapped phase differences and the partial derivatives of the solution. They fit the phase surface to provide the solution, and spread the effect of a singular point to the whole area of imaging data. This may introduce distortion in the regular regions without noise.

Despite many phase unwrapping algorithms have been developed for various applications, there is no agreement among the current phase unwrapping algorithms for different applications [44], because special procedures have to be designed in different algorithms to deal with many different issues due to the existence of disturbance in the measured data. In Doppler OCT, effective phase unwrapping is highly desirable for accurate volume flow quantification and automated data processing. In endovascular Doppler OCT studies, quality-guided phase unwrapping method is used to recover Doppler phase map [45, 46]. But incorrect phase unwrapping happens near the vessel wall due to the existence of high shear rate. Xu et al. present two noise filtering methods [47], the “sine/cosine average filter” and the “phase tracker method” (PTM). The filtered Doppler image can then be unwrapped with less error. Unfortunately [47], did not discuss how to correct phase discontinuities. In this paper, we present the theoretical analysis and experimental result to demonstrate a novel 2D phase unwrapping in Doppler FD-OCT. By calculating phase difference between two adjacent points in the axial direction, real phase information is recovered iteratively. Random phase noise is suppressed by means of window moving average in the sampled Doppler phase image using complex Doppler OCT data. Using the 2D phase unwrapping method, residual discontinuous phase points in the wrapped Doppler OCT phase map can be identified accurately.

2. Materials and Instruments

This research adhered to the tenets of the declaration of Helsinki in the treatment of human subjects. The raw data used below were taken from two healthy subjects (male) in 2010 when two of the authors (Wang and Huang) worked at the University of Southern California (USC). The research protocol was approved by the institutional review board of USC.

A commercial spectrometer-based FD-OCT system was used in this study (RTVue, Optovue Inc.). Our OCT system operated at a center wavelength of 840 nm with an axial resolution of 5 μm. The transverse resolution was about 20 μm as limited by optical diffraction of the eye. Time interval between two sequential axial scans was 36.7 μs. The maximum measurable Doppler frequency shift was 13.6 kHz at phase wrapping limit. This corresponded to a maximum measurable axial velocity component of 4.2 mm/s.

3. Analysis

3.1 Theory

Figure 1(b) shows a Doppler OCT phase image sampled from human retina (image size is 1.0 x 0.96 mm), and Fig. 1(a) is its structural image. For vessel V_r at right side of the image, it has negative Doppler phase signal (black signal) with wrapped phase (white signal) at its central area.

 

Fig. 1 (a) Retina structure image; (b) Doppler OCT phase image; (c) Phase difference Δφ of the adjacent two pixels for the data between points P₁ and P_n in image (b); (d) Plot of Doppler FD-OCT complex data in a complex plane.

Download Full Size | PDF

To describe the algorithm for phase unwrapping, we restore complex data f(z) = A(z)exp(iφ (z)) for each Doppler image pixel using its OCT amplitude A(z) and Doppler phase φ(z). At the position marked as dashed line P₁P_n across vessel V_r in Fig. 1(b), complex value of each pixel from P₁ (at the upper boundary of retina) to P_c (center of the vessel) are plotted in a complex plane in Fig. 1(d), where P₁ and P_c are shown as white circular symbols in phase area IV and II separately. It can be seen that from retinal surface (P₁) to the vessel center (P_c), complex vector of each Doppler OCT data point rotates continuously from vector OP₁ to OP_c clockwise with phase φ increasing gradually from vessel wall to the central area of blood flow.

For two adjacent points P_m and P_m-1 in Fig. 1(d), phase φ_m for point P_m can be expressed as φ_m = φ_m-1 + Δφ, in which φ_m-1 is the phase of point P_m-1, and Δφ is the phase difference between those two points. Thus, in one A-line of Doppler OCT phase map, setting the unwrapped Doppler phase ψ at upper tissue boundary P₁ (with Doppler phase φ₁) as ψ₁ = φ₁, restored phase result ψ_m of the m_th data point in the A-line can be calculated as:

3.2 Random phase noise

In practice, however, the presence of noise in the wrapped Doppler phase map can affect effective phase unwrapping. Noise may make the detection of phase variation among Doppler points erroneous, and lead to wrong unwrapped result. Considering three adjacent points A, B, and C in one Doppler OCT A-line, in which A and C are signal points, and B is the point with random phase noise. To calculate phase value at point C using Eq. (2) through the route A→B→C, the influence of noise B has to be considered. In Fig. 2(a), points A and C are shown in the phase area III of a complex plane. Phase difference Δφ_ca between points C and A is negative, which corresponding to the angle rotation from vector OA to OC clockwise. From point A or C, we draw lines across coordinate center O to point A’ or C’ in phase area I separately. Then the complex plane is divided into four regions, COA, AOC’, C’OA’ and A’OC. If noise B is within the region COA, as shown in Fig. 2(a), obviously, Doppler phase at point C can be determined through A→B→C without mistake.

 

Fig. 2 Effect of noise B on phase unwrapping; (a) B in the region COA; (b) B in the region AOC’; (C) B in the region A’OC; (d) B in the region C’OA’.

Download Full Size | PDF

If noise B falls into the region AOC’, as shown in Fig. 2(b), phase difference Δφ_ba between points B and A is positive, which corresponding to the angle rotation from vector OA to OB counter clockwise. Phase difference Δφ_cb between points C and B is negative, which corresponding to the vector rotation from OB to OC clockwise. Counter clockwise rotation from OA to OB will be compensated by the clockwise rotation from OB to OC, and then we can get the correct phase at point C as ψ_c = φ_a + (Δφ_cb + Δφ_ba). Similarly, if noise B is within the region A’OC, as shown in Fig. 2(c), Doppler phase at point C can still be determined correctly, because the angle rotation from vector OA to OB will be compensated by the reverse rotation from OB to OC. However, if noise B falls into the region A’OC’, as shown in Fig. 2(d), phase vectors rotate continuously in one direction, from OA to OB and from OB to OC counter clockwise. The accumulated rotation angle is larger than π, and the calculated phase difference Δφ_ca between points C and A is positive, which is incorrect.

From Fig. 2, it can be seen that random noise can affect phase unwrapping process if it happens in the opposite region A’OC’ of signal points A and C. To eliminate the influence of random noise, point B should be moved out of the region A’OC’. One simple way to do this is to average Doppler phase points in a complex plane. Random noises generally appear when OCT signals having low signal to noise ratio (SNR), or when the amplitudes of the OCT signals are less than those of normal data. The amplitude weighted average helps the data points with sufficiently large SNR to dominate in the filtering operation, and the effect of random noise can be suppressed. In this study, 3 by 3 window moving average is implemented over the x-z 2D Doppler OCT phase map. For a Doppler image point φ(x, m), its complex value f(x,m) can be averaged as:

3.3 Discontinuous phase point

Blood flow distribution inside a vessel is complex. In Doppler OCT, sampled phase image may contain discontinuous signals. The true phase changes by more than half a cycle (π radians) between two consecutive Doppler points with phase discontinuity, and the wrapped phase will be shifted one cycle. In this case, the calculation of phase difference Δφ_m.m-1 in Eq. (2) isn’t reliable, and phase aliasing will be created, such as unnatural gaps of unwrapped phase with the error of 2Kπ radians (K is an integer number). Thus, discontinuous phase points cannot be restored easily using Eq. (2) only in one dimension. To correct phase discontinuities, 2D phase unwrapping is carried out in the section 4.3 in this study.

4. Results

4.1 Testing window filtering

Wrapped Doppler phase data in Fig. 1(b) from point P₁ to P_n were processed using Eq. (2) firstly, and shown in Fig. 3. Black triangle symbols represent the original wrapped phase φ. The plus symbols show the unwrapped result ψ without using 3x3 window smoothing. After unwrapping, phase jumping in the sampled data is removed and Doppler signal shows a continuous profile. The minimum phase shift is −3.93 radian, which exceeds –π to + π interval. In Fig. 3, the white circular symbols represent the restored Doppler phase after window smoothing with Eq. (3). It can be seen that window filtering can smooth the recovered flow profile effectively without introducing errors in the phase distribution.

 

Fig. 3 Doppler OCT phase profile before and after unwrapping.

Download Full Size | PDF

4.2 Restoring 2D Doppler phase map

Through processing each A-line separately in a phase image after 2D window smoothing, phase unwrapping is done for a vessel in two dimensions. To determine boundary condition φ₁, we set up an amplitude threshold A₀ based on our OCT signal strength. For each A-line, after five-pixel smoothing of the OCT amplitude, retina surface (up-boundary) is searched from image top to the bottom. For the first pixel with the averaged amplitude $\overline{A}(z)>A_0$, its phase, φ(z), is chosen as the boundary value φ₁ = φ(z).

Figure 4 shows 2D Doppler phase maps sampled from human retina vessels, in which column a represents the original wrapped phase images 1, 2, and 3 (images 1 and 2 are from the same vessel sampled at different time). Images in column b are the unwrapped phase result ψ without using 3x3 window smoothing. Column c shows the wrapped Doppler phase after noise filtering, and images in column d are their recovered phase results. In Fig. 4, the grey scale is set from −5 to 5 radians (brightness from 0 to 255). Brightness of pixels with phase less than −5 radians is set as zero. For images b2 and b3, solid black lines can be seen in the recovered phase map, which indicate wrong unwrapping result induced by noise. After noise filtering, those lines disappear in images d2 and d3, and wrapped Doppler OCT phase is recovered successfully.

 

Fig. 4 Column a: Original Doppler phase φ of three vessels; Column b: Phase unwrapping result without noise filtering; Column c: Wrapped phase map after window smoothing; column d: Unwrapped phase result after noise filtering.

Download Full Size | PDF

In image b2, phase unwrapping fails in two A-lines. For the left failed A-line L_rn1, its sampled Doppler phase is plotted in Fig. 5(a) as the white circular symbols, and the black triangle symbols are the unwrapped result without using window smoothing. In Fig. 5(a), the horizontal axis represents the number of phase pixels started from retinal surface axially. The range from number 60 to 70 corresponds to the tissue area under the vessel with background phase signal. There should be no phase wrapping happened there, and the restored phase ψ should be the same as the original phase φ. However, Fig. 5(a) shows a −2π phase difference between ψ and φ in the area under the vessel, which is not correct.

 

Fig. 5 (a) Wrapped (circular symbols) and recovered (triangle symbols) phase for the A-line L_rn1; (b) Doppler phase difference before and after window smoothing; (c) Noise point P_r in a complex plane; (d) Recovered Doppler phase with (triangle symbols) and without (circular symbols) noise filtering.

Download Full Size | PDF

For each pixel in the A-line L_rn1, differences between the window smoothed phase ${\overline{\phi }}_m$ and the original phase φ_m are calculated and plotted in Fig. 5(b), in which the vertical axis represents the phase difference ${\overline{\phi }}_m-{\phi }_m$. It can be seen that Doppler phase at the pixel P_r has a big change (close to –π) after window smoothing. In Fig. 5(c), the adjacent three points P_r-1, P_r, and P_{r + 1} are shown as black dots in a complex plane. Obviously, P_r is in the opposite region of points P_r-1 and P_{r + 1}. As discussed in Fig. 2(d), point P_r can be identified as the noise point which induces −2π phase gap in the unwrapped phase result in the A-line L_rn1. In Fig. 5(c), the white triangle symbols represent those three points after window smoothing, and noise point P_r is moved to the new position ${\overline{P}}_r$. After noise filtering, the recovered phase for the A-line L_rn1 is plotted in Fig. 5(d) as the black triangle symbols, which shows a smooth phase distribution without −2π phase gap. As a comparison, the white circular symbols in Fig. 5(d) represent the failed unwrapping result without window smoothing. Therefore, complex data averaging can suppress phase noise effectively in Doppler phase unwrapping.

For vessel 2 in Fig. 4, its wrapped Doppler OCT phase image a2 is also unwrapped using Xu’s PTM algorithm for comparison. The recovered result is shown in Fig. 6(a). A 3x3 sliding window, which is the same size as that employed in Fig. 4, is used for least square fitting at each image pixel. As can be seen in Fig. 6(a), after PTM denoising, random noises under the vessel (in image b2 of Fig. 4) are suppressed successfully. At the position marked as the dashed line in Fig. 6(a), phase profile is plotted in Fig. 6(b) as black triangle symbols. White square symbols in Fig. 6(b) show the phase profile at the same position in image d2 of Fig. 4. It can be seen that phase unwrapping result using our complex data averaging is comparable with that of PTM algorithm. For vessel 2, the sampled Doppler OCT phase image consists of 500 A-lines (several blood vessels sampled in one B-scan including vessel 2). In the axial direction, the unwrapping depth is set as U = 100 to cover retinal vessels. The processing time for a data set of 500 x 100 pixels was less than 1 second using our complex data averaging algorithm (computer: CPU laptop, Dell Inspiron). However, to get image 6(a) using the same data set with PTM algorithm, it took about 53 seconds using the same CPU due to exhaustive searching of three fitting parameters in each sliding window.

 

Fig. 6 (a) Unwrapped Doppler OCT phase map using PTM algorithm; (b) Recovered Doppler phase profile with PTM method(triangle symbols) and complex averaging algorithm (square symbols).

Download Full Size | PDF

4.3 2D unwrapping for discontinuous phase point

To investigate discontinuous phase points, we choose a Doppler phase map sampled at vessel branching. Figures 7(a) and 7(b) show the original phase image and the unwrapped phase map after 3x3 window smoothing respectively. It can be seen that after denoising, there are still four failed A-lines found in the restored phase image 7(b). They are labeled as line L₁, L₂, L₃ and L₄ separately, in which lines L₂ and L₃ are adjacent to each other.

 

Fig. 7 (a) Sampled Doppler OCT phase image; (b) Restored phase image by means of 3x3 window smoothing.

Download Full Size | PDF

To recover Doppler phase correctly in those failed A-lines, phase unwrapping is carried out in the lateral image direction. The method is described in Fig. 8. For the adjacent two failed A-lines L₂ and L₃, Fig. 8(a) shows four Doppler OCT A-lines L_left, L₂, L₃ and L_right, in which L_left is the A-line at left side of L₂ and L_right is the A-line at right side of L₃. Doppler phase at A-lines L_left and L_right has been recovered successfully, as shown in Fig. 7(b). In Fig. 8(a), the numbers in the vertical direction represent Doppler image pixels started from retinal surface axially, from 1 to U.

 

Fig. 8 (a) Four A-lines distribute laterally; (b) Phase unwrapping result for the A-line L_right. Black dot: Phase result unwrapped axially; White circle: Phase result unwrapped laterally.

Download Full Size | PDF

For the A-lines L_left and L_right, their corrected phase values are ψ_L(m) and ψ_R(m) separately, in which m is the pixel number vertically. Phase ψ_L(m) and ψ_R(m) are used as references (boundary values) to figure out phase discontinuities in A-lines L₂ and L₃. At the j_th image depth (the j_th point axially), setting boundary condition ψ₁ = ψ_L(j) and using Eq. (2), phase unwrapping is performed in the lateral direction through the path L_left(j) → L₂(j) → L₃(j) → L_right(j). The laterally unwrapped results for A-lines L₂, L₃ and L_right are named as β₂(j), β₃(j) and β_R(j) separately. If β_R(j) = ψ_R(j), the lateral unwrapping process at the j_th image depth is successful. Then phase values β₂(j) and β₃(j) will be kept as the recovered results for pixels L₂(j) and L₃(j) respectively. If β_R(j)≠ψ_R(j), the lateral unwrapping process is failed, and the j_th pixels L₂(j) and L₃(j) will be identified as the suspicious points. They may be the singular points which introducing mistakes during phase unwrapping in the axial direction.

Figure 8(b) shows the restored Doppler phase for the A-line L_right, in which the black dots represent the corrected phase ψ_R(m) recovered axially, and the white circular symbols represent phase result β_R(m) got through lateral unwrapping. Comparing phase ψ_R(m) and β_R(m) in Fig. 8(b), it can be seen that they are different at three points, which labeled as j, j + 1 and j + 2. Therefore, for A-lines L₂ and L₃, pixels at these three positions are suspicious points, and marked as black dots in lines L₂ and L₃ in Fig. 8(a). These three points are adjacent to each other and form a suspicious region.

To identify if the suspicious region consists of singular points, vertical unwrapping test is performed. For the failed A-line L₃, its laterally unwrapped result β₃(m) is plotted in Fig. 9 aswhite circular symbols. For the suspicious region from j to j + 2, phase unwrapping is done axially from pixel (j-1) (the point at upside of the suspicious region) to pixel (j + 3) (the point at downside of the suspicious region). Phase β₃(j-1) is used as the boundary value. The unwrapped results for those points are α₃(j), α₃(j + 1), α₃(j + 2) and α₃(j + 3), shown as black dots in Fig. 9. If α₃(j + 3) = β₃(j + 3), it means there is no singular points among those suspicious points, and the suspicious region is induced by phase discontinuities in the lateral unwrapping direction. Values α₃(j), α₃(j + 1) and α₃(j + 2) can be kept as the corrected phase results for those three points. If α₃(j + 3)≠ β₃(j + 3), that means the suspicious region cannot be unwrapped successfully either in the vertical direction, or in the lateral direction. Therefore, the suspicious region is a failed region which consists of singular points. In Fig. 9, obviously, α₃(j + 3) is not equal to β₃(j + 3). Thus, the suspicious region from pixel j to j + 2 is a failed region, as labeled in Fig. 8(a).

 

Fig. 9 Lateral phase unwrapping result for A-line L₃.

Download Full Size | PDF

From Fig. 8(b), it can be seen that after lateral phase unwrapping, Doppler phase at the A-line L₃ has been restored successfully except the pixels in the failed region. To determine true phase values for those failed points, phase unwrapping is done in two different paths axially. The first path is from pixel j-1 to j + 2, and the result is shown as the black dots in Fig. 10. The restored values α₃(j), α₃(j + 1) and α₃(j + 2) are negative, and their summation is α₃(j) + α₃(j + 1) + α₃(j + 2) = −15.31 radians. The other path is from pixel j + 3 to j, and the result is shown as white circular symbols in Fig. 10. The unwrapped phase values at those three points are positive, and their summation is 3.54 radians. For the A-line L_right with corrected phase ψ_R(m) (shown as black dots in Fig. 8(b)), summation of phase values at those three depth is ψ_R(j) + ψ_R(j + 1) + ψ_R(j + 2) = −13.29 radians, which is close to the value −15.31 radians. Considering the continuity of phase distribution in the lateral direction, the unwrapping result through the path j-1 to j + 2 is correct. Thus, true phase information at A-line L₃ is got, which is shown as black dots in Fig. 10. In the restored phase at the A-line L₃, phase values at pixels j + 2 and j + 3 are −6.23 and −1.19 radians separately. The absolute value of their difference is 5.04 radians, which is larger than π. Therefore, the point L₃(j + 2) is the discontinuous phase point, which introducing mistakes on the A-line L₃ during phase recovering axially.

 

Fig. 10 Determination of true phase values for A-line L₃.

Download Full Size | PDF

Using the method described above, phase result for the other three failed A-lines, L₁, L₂ and L₄ are recovered. Restored phase for A-lines L₂ and L₁ are shown in Figs. 11(a) and 10 (b) separately. For the A-line L₄, Fig. 11(c) shows the laterally unwrapped result for the A-line next to it at right side, in which the white circular symbols are the laterally unwrapped result β, and the black dots represent the corrected phase ψ unwrapped axially. They are the same. There is no suspicious point found in Fig. 11(c). Thus, for the A-line L₄, the laterally unwrapped result β₄(m) is chosen directly as the recovered phase and shown in Fig. 11(d). From the restored phase of those three A-lines, it can be seen that absolute values of phase difference between points P_d1 and P_d2 in Fig. (a), P_d3 and P_d4 in Fig. (b), and P_d5 and P_d6 in Fig. (c) are all larger than π. Therefore, points P_d1, P_d3 and P_d5 are discontinuous points, and abrupt phase changes happen after them.

 

Fig. 11 (a) Corrected phase result for the A-line L₂; (b) Corrected phase result for the A-line L₁; (c) Searching failed points for the A-line L₄; (d) Corrected phase result for the A-line L₄.

Download Full Size | PDF

After correction of discontinuous points in the failed A-lines in Fig. 7(b), its Doppler phase map is shown in Fig. 12(a). Wrapped phase core inside the blood vessel disappears, and the grey scale is extended to the range from −7 to 7 radians. Figure 12(b) shows its three dimensional phase distribution. There is no wrapped phase found on the recovered flow profile. The execution of the proposed algorithm is described in a flow chart in Fig. 13.

 

Fig. 12 (a) Doppler OCT phase map after 2D unwrapping; (b) 3D phase distribution at the area marked as a dashed window in (a); Unit for the grey scale bar: radian.

Download Full Size | PDF

 

Fig. 13 Flow chart of the proposed algorithm.

Download Full Size | PDF

4.4 Phase unwrapping for multiple vessels

To demonstrate the feasibility of the proposed method, we investigated phase unwrapping for the Doppler phase map containing multiple blood vessels. A circular scanning pattern was employed to acquire Doppler OCT signal from human retina. The scanning circle was centered at the optic nerve head and had a radius of 1.8mm. Each circle transected all branch retinal arteries and veins. There were 4000 A-lines sampled in each circle. The phase difference for every four A-lines were calculated and averaged to get the Doppler flow information. Thus, each OCT frame consisted of 1000 vertical lines.

Figure 14(a) shows the retina structure image sampled using circular scan. The image size is 11.3 mm x 1.27 mm (horizontal x vertical). Figure 14(b) is the acquired Doppler OCT phase image. Doppler signal induced by moving blood within major branch retinal vessels around the optic nerve head is visible in this image.

 

Fig. 14 (a) Retina structure image; (b) Doppler OCT image; (c) Unwrapped phase with discontinuous points; (d) Recovered phase image after the correction of discontinuous points.

Download Full Size | PDF

In the Doppler phase image in Fig. 14(b), three vessels, V₁, V₂ and V₃ have wrapped phase signal. Detailed phase structure of them can be seen in the enlarged windows (magnification 3:1) in Fig. 14(b). Vessels V₁ and V₂ have positive phase signal with negative phase wrapping, and vessel V₃ has negative phase signal with positive wrapping. Because phase data in image (b) have been averaged for every four A-lines laterally during Doppler data processing, we choose 1 by 3 window (smoothing axially) to suppress random noise. Figure 14(c) shows unwrapped phase image after noise filtering. One failed line L_f can be identified as white solid line in the recovered image. After correction of discontinuous points for the A-line L_f, the restored Doppler phase map is shown in Fig. 14(d). No wrapped phase information can be found in vessels V₁, V₂ and V₃. Doppler phase map with multiple vessels containing wrapped phase signal is corrected successfully.

5. Discussion

There are many different methods developed to do phase unwrapping in interferometric imaging field. Most of them are divided into two basic classes. The first class is the path-following method. The second is the least-square method, which fits the phase surface to provide the solution. In Doppler OCT, least-square fitting has been employed in Xu’s study [47]. The method described in this study is based on phase gradient calculation through processing complex Doppler OCT data, and there is no need to formulate the unwrapping solution. The approach to search discontinuous phase points may be able to be classed into path-following method.

A major challenge for effective phase unwrapping is to reduce phase errors induced by noise. An effective denoising method should have the ability to suppress random phase noise, while preserving the abrupt phase changes induced by flow discontinuities for accurate phase unwrapping. In Doppler OCT, random noise usually predominates at the area with low OCT signal, for example in the shadow under a big blood vessel. This can be seen in images b2 and b3 in Fig. 4. After 2D window moving average using complex Doppler OCT data, random noise is suppressed successfully. In Doppler OCT study, one method for volume blood flow calculation is through integrating the Doppler phase signal over a blood vessel [15, 16]. In Fig. 4, images in column a show the original wrapped Doppler OCT phase φ. Integration over those images is meaningless because phase profile inside those images does not show the real flow profile. Phase image b1 is restored without doing window smoothing. It contains the real flow information for vessel 1, and is used as the reference. Integration of the phase data in the whole area of image b1 is −2668.5 radians. For image d1, the integration result is −2664.5 radians. Therefore, after 3 by 3 window moving average, the restored volume phase accuracy is −2664.5/(−2668.5) = 99.9%.

Note that in Doppler OCT, complex data averaging has been used to increase phase accuracy [28, 48]. Particularly in Szkulmowska’s paper [48], complex data averaging is used to avoid blood velocity underestimation during Doppler data processing, however only Doppler phase is used in the study. How to do phase unwrapping was not investigated in the paper, as can be seen in Figs. 4c and 4h in [48]. On the other hand, in our phase unwrapping studies, both Doppler OCT amplitude and phase are used. We find that complex Doppler data averaging can suppress random phase noise effectively without altering phase gradient, and the abrupt phase change induced by discontinuous points is not smeared out, showing that our phase smoothing method is not a simple low pass filter. After denoising, phase points at background tissue outside of a blood vessel are determined as the reliable points. We choose the pixel at up retinal surface as the starting point, and use the pixels under a blood vessel as the criterion to evaluate the reliability of unwrapping result in one A-line. This simplifies the process to determine boundary conditions.

After filtering of random phase noise, mistakes in the unwrapped phase data are introduced by discontinuous residues with abrupt phase change. In Doppler OCT, phase discontinuity can be induced by shear rates, flow deformation, and high flow rates. At discontinuous points, trend of the flow speed distribution is changed, and the absolute value of phase change is larger than π. It is impossible to correct them through phase derivative analysis only in the axial dimension. In this study, the 2D phase unwrapping is employed to investigate phase discontinuities.

In vascular circulation study with Doppler OCT, effective phase unwrapping is highly desirable. One advantage of this unwrapping method is its simplicity to program and implement. There is no complex algorithm employed. Doppler OCT phase can be recovered accurately without losing original phase information, even for discontinuous points in the deformed blood flow. In Xu’s phase unwrapping study [47], the tubing diameter is 1.6 mm with lateral OCT sampling step of 3 μm, corresponding to effective 533 A-lines in each OCT image. Assuming each vertical image pixel is also 3 μm (this number is not given in [47]), each effective data set would be 533 x 533 pixels in two dimensions, and the corresponding processing time is about 3 seconds [47]. As a comparison, the time for processing a similar data set with 500 x 500 pixels without discontinuous points (because in [47], discontinuous points were not considered) using our method is about 0.9 seconds (including pixel smoothing and unwrapping), less than 1/3 of that in [47]. In the comparison above, we used the CPU of a low-end laptop computer (Dell Inspiron purchased around 2010), while the processor used in [47] is a GeForce GTX750Ti graphics processing unit (GPU) from nVidia Inc.

In order to recover Doppler phase successfully using Eq. (2), physical displacement between two adjacent Doppler image points should be sufficiently small, so that their phase difference Δφ is within the range –π to + π. In the axial direction, this condition can be satisfied if the OCT resolution is sufficiently fine. For example, for a blood vessel (radius R = 50 μm) with a central maximum flow speed of 20 mm/s, at a Doppler angle of 70 degrees, assuming a parabolic flow profile [49], the maximum flow speed is 6840 μm/s at the OCT probe beam direction. The detected flow speed has a function V(r) = −2.74r² + 6840 μm/s, where r is the radius coordinate of the vessel. The first derivative of this speed function is dV/dr = −5.48r, with a maximum value dV/dr|_max = 274 μm/s per micrometers at the vessel wall (when r = R). In our FD-OCT system, the maximum measurable velocity is V_max = 4.2 mm/s at phase π limit, and the corresponding axial resolution requirement is 15.3 μm (V_max/(dV/dr|_max)). Therefore, in the axial direction, if the physical displacement between two adjacent pixels were less than 15.3 μm, this method could be safe to be used. Such a condition can be readily satisfied in our study, because the axial resolution of our OCT system is 5 μm, much less than 15.3 μm.

6. Conclusions

In conclusion, we have investigated a novel 2D method for phase unwrapping in Doppler FD-OCT based on phase gradient approach. Random phase noise is removed through 2D window moving average using complex Doppler OCT data. Phase gradient is preserved during complex data averaging. Discontinuous points in the wrapped Doppler OCT phase map, which induce errors in the restored phase, can be corrected accurately through the proposed 2D unwrapping process. Wrapped Doppler phase map with multiple retinal blood vessels is recovered successfully. The algorithm is computationally efficient and easy to be implemented. Its execution time varies depending on the number of discontinuous phase points. As one of the several options for phase unwrapping in Doppler OCT, we anticipate that the method and results presented in this paper will benefit the applications of various phase related OCT techniques. Future work will be carried out to improve the performance of the proposed approach through comparing it with other existing phase unwrapping methods, and extend its applications to other phase related OCT techniques.