Research Papers: Imaging

Robust phase unwrapping for phase images in Fourier domain Doppler optical coherence tomography

[+] Author Affiliations
Yong Huang, Shizhao Peng, Yanfeng Wu, Xiaodi Tan

Beijing Institute of Technology, School of Optoelectronics, Beijing, China

J. Biomed. Opt. 22(3), 036014 (Mar 29, 2017). doi:10.1117/1.JBO.22.3.036014
History: Received November 22, 2016; Accepted February 23, 2017
Text Size: A A A

Open Access Open Access

Abstract.  To solve the 2π phase ambiguity for phase-resolved Doppler images in Doppler optical coherence tomography, we present a modified network programming technique for the first time to the best of our knowledge. The proposed method assumes that error of the discrete derivatives between unwrapped phase image and wrapped phase image can be arbitrary values instead of integer-multiple of 2π, which makes the real-phase restoration accurate and robust against noise. We compared our proposed method with the network programming method. Parameters including root-mean-square-error and noise amplification degree were adopted for comparison. The experimental study on simulated images, phantom, and real-vessel OCT images were performed. The proposed method consistently achieves optimal results.

Figures in this Article

Accurate phase extraction can provide valuable information for various sensing and imaging techniques. Currently, phase images have been widely used in the following areas: interferometric synthetic aperture radar,1,2 magnetic resonance imaging,3,4 fringe projection profilometry,58 digital holography,912 microscopic imaging,1315 and phase-resolved Doppler optical coherence tomography (DOCT).1618

DOCT is a noninvasive high-resolution three-dimensional imaging modality and has been widely used for various flow studies.1923 However, phase-resolved Doppler images often suffer from wrapping issue that limits phase values in the range of (π,π]. Deviation from real phase profile will occur if measured phase range is over 2π. Thus, phase unwrapping is an essential and important process for accurate phase profile reconstruction.

To solve phase wrapping issue, various methods have been proposed. Generally, these methods include path-following method,2426 minimum-norm method,2729 multiwavelength method,3033 and other algorithms.3439 Usually preprocess and in-process noise reduction techniques are combined with these methods.4042

For path-following method, Goldstein et al.’s24 algorithm is a classic one, which walks along the path with balanced residues using branch cut algorithm. This method requires a starting point to build the unwrapping path for integration. Thus, different starting points may result in different unwrapping phases. Usually some preprocess noise reduction technique and quality mask are adopted to guide the unwrapping path to deal with noise.25,26 However, it is difficult to obtain a reasonable integral path under high noise.

Least-squares27 is the most common algorithm in minimum-norm method. To solve this minimization, fast Fourier transform and discrete cosine transform operations are used for unwrapping. Later to reduce the noise influence, weighted minimum-norm algorithm was proposed.28 In the weighted minimum-norm methods, network flow algorithm is the most attractive, which was proposed by Costantini29 for phase unwrapping in SAR interferometry. Network flow algorithm assumes that the error of the derivatives between the unwrapped phase and wrapped phase is an integer multiple of 2π. Thus, the phase wrapping problem can be regarded as the linear programming problem and solved by network flow approach. Although this method can prevent errors from spreading, the assumption is unreasonable due to the existence of noise.

Multiple-wavelength method generates a longer synthetic wavelength to increase the phase measurement range using two or multiple phase images.3033 This method relies upon pixel-to-pixel comparison and thus is very sensitive to noise. Other methods mainly include Chinese remainder theorem (CRT) algorithms, intelligent algorithms, and combined noise reduction technique for phase unwrapping. CRT can be used for measuring profilometry of the microphase34,35 because phase unwrapping can be treated as a congruence problem. Same as other methods, this method shows limitation when faced with large amounts of noise.

In intelligent algorithms, many heuristic algorithms, such as genetic algorithm,36 cellular-automata,37 and swarm intelligence,38 are introduced for unwrapping optimization. Path-independent phase unwrapping using total-variation denoising adopts noise reduction techniques during the phase unwrapping process.39

To solve phase wrapping problem in OCT phase-resolved Doppler image, we propose a robust algorithm based on network programming method. In the proposed method, we assume that the error between unwrapped phase image derivatives and wrapped phase image derivatives can be any value, not necessarily the integer multiple of 2π. In fact, this hypothesis is more reasonable for wrapped phase images with large noise in real situations. Then we transform the phase unwrapping problem into an optimization problem to be solved.

Acquisition of Unwrapped Phase Gradients

Due to the 2π ambiguity of phase-resolved Doppler image, the relationship between true phase image ϕ and wrapped phase image φ for every pixel can be expressed as Display Formula

φ(i,j)=W[ϕ(i,j)]=ϕ(i,j)+2πk(i,j)(i=1,,m;j=1,,n),(1)
where (i,j) denotes position of every pixel, k(i,j) is an arbitrary integer at pixel location (i,j), and m and n are image height and width, respectively. W is the wrapping operator, which denotes modulo operation on true phase image to get wrapped phase image. It should be stated that we use mod(ϕ+π,2π)π (mod is modulo operation) to realize the effect of wrapping operator W. After wrapping operation, phase is limited to the range (π,π]. From Eq. (1), one can see that if no noise is present, unwrapped phase image can be obtained by adding an integer multiple of 2π to wrapped phase image for every pixel.

Here, we use the gradient operation G, which denotes image derivative and can be performed in either x or y direction, to process the both sides of Eq. (1) Display Formula

G[φ(i,j)]=G[ϕ(i,j)+2πk(i,j)].(2)

Since gradient operation G is a linear operator, Eq. (2) becomes Display Formula

G[φ(i,j)]=G[ϕ(i,j)]+G[2πk(i,j)].(3)

It is obvious that G[2πk(i,j)] still is integer multiple of 2π, so Eq. (3) can be simplified as Display Formula

G[φ(i,j)]=G[ϕ(i,j)]+2πk1(i,j),(4)
where k1(i,j) is an integer.

Then, we use the unwrapping operator U (U is the same modulo operation as W mathematically. We define U as unwrapping operator to give the reader a better understanding and emphasis of the unwrapping process.) Display Formula

U{G[φ(i,j)]}=U{G[ϕ(i,j)]+2πk1(i,j)}.(5)

It can be seen from Eq. (1) that the meaning of unwrapping operator U is that an integer multiple of 2π is added to the wrapped object, thus Eq. (5) can be simplified as Display Formula

U{G[φ(i,j)]}={G[ϕ(i,j)]+2πk1(i,j)}+2πk2(i,j),(6)
where k2(i,j) is an integer.

With the simplified structure, Eq. (6) can be written as Display Formula

U{G[φ(i,j)]}=G[ϕ(i,j)]+2π[k1(i,j)+k2(i,j)].(7)

In Eq. (7), the value range of U{G[φ(i,j)]} is (π,π] due to the use of unwrapping operator U. In general, true phase image ϕ is continuous, thus the value range of G[ϕ(i,j)] is also (π,π]. So, to keep the same value range for both sides of Eq. (7), we can conclude that Display Formula

k1(i,j)+k2(i,j)=0,  (i,j).(8)

Based on the conclusion of Eqs. (8) and (7) it can be written as Display Formula

G[ϕ(i,j)]=U{G[φ(i,j)]}.(9)

From Eq. (9), we can know that the derivatives of wrapped phase image after using unwrapping operator U are equal to the derivatives of true phase image. The relationship of Eq. (9) can be extended in the directions of x and y, which can be expressed as φx(i,j)=Δϕx(i,j), φy(i,j)=Δϕy(i,j), where φx,φy denotes the derivative of wrapped phase image after using unwrapping operator U and Δϕx,Δϕy denotes the derivative of true phase image. Specifically, their derivatives are given by Display Formula

φx(i,j)=U[φ(i+1,j)φ(i,j)],(10)
Display Formula
φy(i,j)=U[φ(i,j+1)φ(i,j)],(11)
Display Formula
Δϕx(i,j)=ϕ(i+1,j)ϕ(i,j),(12)
Display Formula
Δϕy(i,j)=ϕ(i,j+1)ϕ(i,j).(13)

Figure 1 shows the special relationship between true phase image and wrapped phase image directly. Figure 1(a) is the true phase image and Figs. 1(a-1) and 1(a-2) are the corresponding gradients in x and y directions. Figure 1(b) is wrapped phase image and Figs. 1(b-1) and 1(b-2) are the corresponding gradients in x and y directions. From Figs. 1(a) and 1(b), we can clearly see that after using unwrapping operator, the wrapped phase image gradients are equal to the true phase image gradients, as Fig. 1(a-1)is equal to Figs. 1(b-3) and 1(a-2) is equal to Fig. 1(b-4). Figure 1(c) shows gradient along the y direction of row 20 in Figs. 1(a-2) and 1(b-2). We can see that unwrapping operator can restore big gradient values caused by the phase jumps.

Graphic Jump Location
Fig. 1
F1 :

Gradient analysis of simulated phase image with no noise. (a) True phase image, (a-1) gradient of true phase image in the x direction, (a-2) gradient of true phase image in the y direction, (b) wrapped phase image, (b-1) gradient of wrapped phase image in the x direction, (b-2) gradient of wrapped phase image in the y direction, (b-3) gradient in the x direction after unwrapping operator, (b-4) gradient in the y direction after unwrapping operator, (c) gradient along the y direction of row 20 in (a-2) and (b-2). Red line denotes true phase image gradient and black line denotes wrapped phase image gradient with phase jumps. Black arrow and red arrow denote the phase jump and unwrapping operation process.

However, when the wrapped phase image is corrupted by noise, it is impossible to guarantee Eq. (1). Therefore, error term e(i,j) caused by noise is introduced into Eq. (9). So Eq. (9) becomes Display Formula

G[ϕ(i,j)]=U{G[φ(i,j)]}+e(i,j).(14)

So the error amounts between φx,y and Δϕx,y can be expressed as Display Formula

ex(i,j)=Δϕx(i,j)φx(i,j),ey(i,j)=Δϕy(i,j)φy(i,j).(15)

Due to the uncertain property of noise, error amounts ex,ey can be any values, instead of integer multiple of 2π. Dropping the constraint for error amounts, unwrapping process becomes natural compared to network programming algorithm.

Based on Eq. (15), constraints for adjacent pixels can be given by Display Formula

ex(i,j+1)ex(i,j)ey(i+1,j)+ey(i,j)=[Δϕx(i,j+1)Δϕx(i,j)Δϕy(i+1,j)+Δϕy(i,j)][φx(i,j+1)φx(i,j)φy(i+1,j)+φy(i,j)].(16)

From Eqs. (12) and (13), we can calculate that both Δϕx(i,j+1)Δϕx(i,j) and Δϕy(i+1,j)Δϕy(i,j) are equal to ϕ(i+1,j+1)ϕ(i+1,j)ϕ(i,j+1)+ϕ(i,j). Thus Δϕx(i,j+1)Δϕx(i,j)Δϕy(i+1,j)+Δϕy(i,j) equals 0 for every pixel. Hence, Eq. (16) can be simplified as Display Formula

ex(i,j+1)ex(i,j)ey(i+1,j)+ey(i,j)=R(i,j),(17)
where R(i,j) is [φx(i,j+1)φx(i,j)φy(i+1,j)+φy(i,j)].

As the error values ex,ey is limited, thus process of solving ex,ey can be treated as the following optimization problem: Display Formula

mini,jcx(i,j)|ex(i,j)|+i,jcy(i,j)|ey(i,j)|s.t  ex(i,j+1)ex(i,j)ey(i+1,j)+ey(i,j)=R(i,j),(18)
where cx,cy are weighted coefficients to balance the objective function,29 and R(i,j) is [φx(i,j+1)φx(i,j)φy(i+1,j)+φy(i,j)]. To solve this problem, we can turn Eq. (18) into linear programming problem and use certain optimizing tool to obtain ex,ey values. Equation (18) uses the error terms ex(i,j),ey(i,j) to share the values R(i,j) caused by noise. Thus, the proposed unwrapping method can tolerate the effect of phase noise during the unwrapping process.

Unwrapped Phase Restoration

In Sec. 2.1, every pixel can be affected by error amounts ex,ey, so estimation of the gradient of true phase image Δϕx,Δϕy can be obtained by Eq. (15). Eventually, estimated true phase image, which we refer as unwrapped phase image in later description, can be restored using either one of the following two equations: Display Formula

ϕ(i,j)=φ(i,j)+jΔϕ1(1,j)+iΔϕ2(i,j),(19)
Display Formula
ϕ(i,j)=φ(i,j)+iΔϕ2(i,1)+jΔϕ1(i,j).(20)

Equations (19) and (20) are equivalent. In this paper, we use Eq. (20). Theoretically, unwrapped phase image should be the same as wrapped phase image except for wrapping regions. However, when wrapped phase image contains large noise, the estimated gradients of some pixels are not identical to the true phase image gradients, thus gradient differences caused by noise are accumulated for latter pixels in the integration path. As a result, unwrapped phase image will have a constant shift relative to the true phase image. As the final step for the phase unwrapping, we applied median filtering to remove the relative constant shift.

Evaluation Metrics for Phase Unwrapping

Following four parameters are used to compare the unwrapping results of different methods objectively.

Root-mean-square-error

It can be used to describe the closeness of unwrapped phase image I to the true phase image ϕ. Display Formula

RMSE1=i=1mj=1n(Ii,jϕi,j)2m×n,(21)
where i and j denote pixel position, m and n are the height and width of the image, respectively.

Noise amplification degree

It is used to evaluate the effect of unwrapping algorithm on noise. Noise amplification degree (NAD) is defined as Display Formula

NAD=i=1mj=1n|Ii,jϕi,j|i=1mj=1n|ϕi,jLi,j|,(22)
where Li,j is the phase values of the noisy true phase image with noise added in the position i,j. i=1mj=1n|Ii,jϕi,j| describes the closeness of the unwrapped phase image I to the true phase image ϕ. i=1mj=1n|Ii,jLi,j| describes noise intensity added to the true phase image, thus NAD reflects the ability of noise amplification. A good unwrapping algorithm should not amplify noise. Thus, the value of NAD should be less than or equal to 1.

Residual wrapped map

In situations where it is hard to obtain the true phase image, root-mean-square-error (RMSE) and NAD cannot be calculated. Therefore, we can rewrap unwrapped phase image I to (π,π] and calculate the residual wrapped map between rewrapped phase image and wrapped phase image. Residual wrapped map at pixel location (i,j) can be expressed as Display Formula

E(i,j)=W[I(i,j)]φ(i,j).(23)

In general, if the rewrapped phase image is close to the original wrapped phase image or this residual wrapped map is close to 0, it means that the unwrapping algorithm has a good performance. Thus, RMSE between rewrapped phase image and original wrapped phase image is defined as Display Formula

RMSE2=i=1mj=1nEi,j2m×n.(24)

In cases where we do not have access to true phase image, we can use RMSE2 to evaluate the unwrapping performance.

Execution time

Execution time should be considered in the implemented programming and are based on the average value of 10 runs. We need this time to be as short as possible.

RMSE1, NAD, and execution time are used to evaluate performance against simulated images. RMSE2 and execution time are used to test the unwrapping results of OCT phase-resolved Doppler image.

To demonstrate the performance of unwrapping method, we present the results against simulated images. In this section, network programming and our proposed method are tested on simulated images with random impulse noise as this noise has the large prevalence in phase images of OCT system.43 In the simulation experiment, three evaluation indexes (RMSE1, NAD, and execution time) are used to compare unwrapping results of two methods. For the simulation experiment, all unwrapping process are conducted through MATLAB R2013 in a computer with the Intel® core i7-4790 processer @ 3.6 GHz and 8GB RAM.

First, 0% to 90% random-valued impulse noise is added to Fig. 2(a) of true phase image with the multiple phase regions. Then, true phase image with noise and wrapped phase image is shown in Figs. 2(b) and 2(c). From Fig. 2(d), network programming method can only produce acceptable results for wrapped phase image with 0% or 10% random impulse noise. For the wrapped phase image with higher noise, many visible lines artifacts are generated in the unwrapped phase image as the assumption of error values is integer multiple of 2π. Figure 2(e) shows that our proposed method can achieve good unwrapping results for noise level less than 70%. Even if wrapped phase image is seriously corrupted by 60% random impulse noise, meaningful unwrapped image can still be obtained. Moreover, we are also surprised to find that unwrapped phase images have lower noise level than the corresponding true phase image with noise. As for the reason of noise reduction property, we have analyzed that in Sec. 2.1.

Graphic Jump Location
Fig. 2
F2 :

Unwrapping results of network programming and our proposed method for the simulated image with different random-valued impulse noise levels: (a) true phase image, (b) true phase image with noise, (c) wrapped phase image, (d) unwrapped results using network programming method, and (e) unwrapped results using our proposed method.

Then, RMSE1, NAD values, and execution time are used to evaluate the performance of unwrapping method as shown in Fig. 3. When the random impulse noise level is 10%, RMSE1 and NAD values of the network programming are small. But, RMSE1 and NAD values increase for 20% to 90% noise levels. For every noise level, the proposed method consistently obtains the smaller RMSE1 and NAD values than network programming. What is more, the NAD values of our proposed method are less than 1 all the time, which shows that the proposed method has an extra-advantage of weakening noise. For the time cost, the proposed method is similar to network programming. From these compared results, we can clearly see that the proposed method has better unwrapping performance than network programming method.

Graphic Jump Location
Fig. 3
F3 :

Comparison of (a) RMSE1, (b) NAD, and (c) average execution time of 10 runs using network programming and our proposed methods for images with different random-valued impulse noise levels.

In the previous section, network programming and our proposed method are tested for simulated images and the comparison results show that the proposed method can always obtain the best unwrapping results. To show the unwrapping ability on phase-resolved Doppler OCT images, we select two groups of wrapped phase images (phantom and mouse artery images) from two different OCT systems. Similar to simulation process, we still compare the unwrapping results of network programming and our proposed method. In the evaluation process, phase unwrapping metrics (RMSE2 and execution time) are used to analyze the unwrapping results. For the program environment, computer configuration of all the unwrapping process is the same as simulation experiments.

Unwrapping Results of Phantom Images

Phantom comes from transparent plastic tubes (parameters: 0.5-mm inner diameter, 0.9-mm outer diameter) with flowing milk pumped by the precision pump. The OCT system for phantom imaging is a home-built spectral domain OCT system, which operates at a central wavelength of 1300 nm with a bandwidth of 60 nm, 70 fps with 1000 A-scans per frame and has a measured axial resolution of 14  μm, imaging range of 6.7 mm. The sensitivity and phase stability of the OCT system are 92 dB and 70 mrad, respectively.

To illustrate the unwrapping performance of two methods for phantom with different phase wrapping degrees, Fig. 4 shows seven groups phantom images with different flux levels. Figures 4(a) and 4(b) are the intensity and phase images, respectively. Figures 4(c) and 4(d) are the unwrapped results of network programming and our proposed method for Fig. 4(b). It is clear from unwrapped phase image that network programming cannot get satisfactory unwrapped phase images. The unwrapped result of our proposed method is satisfactory since unwrapped phase image achieves the excellent continuity and contrast. Thus, these visual comparison results indicate that our proposed method has the outstanding unwrapping performance compared to the network programming for OCT phase images.

Graphic Jump Location
Fig. 4
F4 :

Phase unwrapping result comparison for phantom images with different flux levels (251×841  pixels). (a) Intensity image, (b) phase-resolved Doppler image, (c) unwrapped phase images with network programming method, and (d) unwrapped phase images with our proposed method. (Scale bar: 250  μm.)

Figure 5 shows the RMSE2 and time consumption of two unwrapping methods. From the RMSE2 curve, we can see that our proposed method still consistently achieves the minimum unwrapping error for different phantom images. For the execution time, we can see there is a large rise up near the boundary where phase wrapping is about to occur (Q=1100  μl/min and Q=1150  μl/min). The reason might be that the algorithm requires more time to find the optimal solution for pixels near the boundary as they can be either random noise or phase wrapping points. For other flux conditions, the time cost of our proposed method is less than the network programming method.

Graphic Jump Location
Fig. 5
F5 :

Comparison of (a) RMSE2 and (b) average execution time of 10 runs using network programming and our proposed method for phantom images with different flux levels.

Unwrapping Results of Vessel Images

Eventually, we applied network programming method and our proposed method to four groups of the real biological sample images, which are mouse femoral artery images, shown in Fig. 6. These images were acquired from a separate OCT system reported in Ref. 44 to demonstrate the applicability of our proposed method.

Graphic Jump Location
Fig. 6
F6 :

Phase unwrapping results using network programming and our proposed method for vessel phase images with different wrapping degrees (200×999  pixels). (a) Vessel intensity image, (b) phase-resolved Doppler image, (c) unwrapped phase image using network programming method, and (d) unwrapped phase image using our proposed method. (Scale bar: 250  μm.)

Figures 6(a) and 6(b) show the intensity image and Doppler phase image of mouse vessels, respectively. Note that phase images of Fig. 6(b) present different wrapping degrees, which are caused by different blood flow of mouse artery. Figures 6(c) and 6(d) are unwrapped phase images of network programming and our proposed method. From the result of Fig. 6(c), unwrapped phase image of network programming is still poor, thus network programming shows limitation when facing large random noise in OCT phase images. However, the proposed method can tolerate noise well. Figure 6(d) shows that the proposed method obtains attractive unwrapped results no matter which group of phase images.

Similarly, Fig. 7 shows RMSE2 values and time cost of these two unwrapping methods. From the RMSE2 curve, network programming has the small RMSE2 value for the third vessel phase image, but we can observe from Fig. 6(c) that this RMSE2 value of network programming is not a sign of good unwrapping results. After the RMSE2 comparison, our proposed method still has minimum unwrapping error compared to network programming. For time cost, our proposed method has lower cost than network programming.

Graphic Jump Location
Fig. 7
F7 :

Comparison of (a) RMSE2 and (b) average execution time of 10 runs using network programming and our proposed method for vascular images with different wrapping degrees.

We have demonstrated the capability of our proposed method to solve the 2π phase ambiguity issue for phase-resolved Doppler images in the Fourier domain Doppler OCT systems. It is observed from the experimental test that the proposed method has excellent performance of unwrapping phase and weakening noise. It can provide accurate unwrapped phase images.

In summary, the first advantage of our proposed method is prefiltering free. In other common unwrapping methods such as Goldstein method, least-squares with fast Fourier transform methods, and path independent with total variation denoising method, prefiltering process is of great interest to the whole unwrapping process. However, it is difficult to design perfect filtering algorithms suitable for unwrapping process. Because on one hand noises need to be removed as much as possible. On the other hand, accurate phase image gradient needs to be preserved, which is key to post true phase reconstruction. From this perspective, we can see that the proposed method does not need the prefiltering process and also achieves stable and good results. From the unwrapped phase images of phantom and mouse artery, we can clearly see that unwrapped phase images of the proposed method have good continuity and low noise level. General filtering can be applied to improve the image quality after using our unwrapping method. The second advantage of our proposed method is the property of weakening noise. The simulation experiment shows that NAD values of the proposed method are less than 1, which means the proposed approach can not only tolerate the influence of noise but also reduce noise, which comes from the procedure of building the error relationship and solving the error values in the optimization process.

Due to the requirement of optimizing process, the major limitation of current proposed method is relative high time cost. Our future work will focus on the reduction of the time cost.

The authors have no relevant financial interests in the paper and no other potential conflicts of interest to disclose.

The authors gratefully acknowledge funding of National Natural Science Foundation of China (NSFC, Project Number 61505006) and Startup Plan for Young Teachers from Beijing Institute of Technology.

Zebker  H. A., and Lu  Y., “Phase unwrapping algorithms for radar interferometry: residue-cut, least-squares, and synthesis algorithms,” J. Opt. Soc. Am. A. 15, (3 ), 586 –598 (1998). 0740-3232 CrossRef
Goldstein  R. M., , Zebker  H. A., and Werner  C. L., “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci.. 23, (4 ), 713 –720 (1988).CrossRef
Bagher-Ebadian  H., , Jiang  Q., and Ewing  J. R., “A modified Fourier-based phase unwrapping algorithm with an application to MRI venography,” J. Magn. Reson. Imaging. 27, (3 ), 649 –652 (2008).CrossRef
Chavez  S., , Xiang  Q. S., and An  L., “Understanding phase maps in MRI: a new cutline phase unwrapping method,” IEEE Trans. Med. Imag.. 21, (8 ), 966 –977 (2002). 0278-0062 CrossRef
Li  S.  et al., “Two-dimensional wavelet transform for reliability-guided phase unwrapping in optical fringe pattern analysis,” Appl. Opt.. 51, (12 ), 2026 –2034 (2012). 0003-6935 CrossRef
Xu  Y.  et al., “Recovery of absolute height from wrapped phase maps for fringe projection profilometry,” Opt. Express. 22, (14 ), 16819 –16828 (2014). 1094-4087 CrossRef
Perciante  C. D.  et al., “Wrapping-free phase retrieval with applications to interferometry, 3D-shape profiling, and deflectometry,” Appl. Opt.. 54, (10 ), 3018 –3023 (2015) 0003-6935 CrossRef
Cheng  Z.  et al., “Practical phase unwrapping of interferometric fringes based on unscented Kalman filter technique,” Opt. Express. 23, (25 ), 32337 –32349 (2015). 1094-4087 CrossRef
Ma  L. H.  et al., “Fast algorithm for reliability-guided phase unwrapping in digital holographic microscopy,” Appl. Opt.. 51, (36 ), 8800 –8807 (2012). 0003-6935 CrossRef
Iglesias  I., “Phase estimation from digital holograms without unwrapping,” Opt. Express. 22, (18 ), 21340 –21346 (2014). 1094-4087 CrossRef
Khodadad  D., “Phase-derivative-based estimation of a digital reference wave from a single off-axis digital hologram,” Appl. Opt.. 55, (7 ), 1663 –1669 (2016). 0003-6935 CrossRef
Jaedicke  V.  et al., “Multiwavelength phase unwrapping and aberration correction using depth filtered digital holography,” Opt. Lett.. 39, (14 ), 4160 –4163 (2014). 0146-9592 CrossRef
van den Doel  L. R., and van Vliet  L. J., “Temporal phase-unwrapping algorithm for dynamic interference pattern analysis in interference-contrast microscopy,” Appl. Opt.. 40, (25 ), 4487 –4500 (2001). 0003-6935 CrossRef
Yaqoob  Z.  et al., “Improved phase sensitivity in spectral domain phase microscopy using line-field illumination and self phase-referencing,” Opt. Express. 17, (13 ), 10681 –10687 (2009). 1094-4087 CrossRef
Watanabe  E., , Hoshiba  T., and Javidi  B., “High-precision microscopic phase imaging without phase unwrapping for cancer cell identification,” Opt. Lett.. 38, (8 ), 1319 –1321 (2013). 0146-9592 CrossRef
Morofke  D.  et al., “Wide dynamic range detection of bidirectional flow in Doppler optical coherence tomography using a two-dimensional Kasai estimator,” Opt. Lett.. 32, (3 ), 253 –255 (2007). 0146-9592 CrossRef
Hendargo  H. C.  et al., “Doppler velocity detection limitations in spectrometer-based versus swept-source optical coherence tomography,” Biomed. Opt. Express. 2, (8 ), 2175 –2188 (2011). 2156-7085 CrossRef
Sun  C.  et al., “In vivo feasibility of endovascular Doppler optical coherence tomography,” Biomed. Opt. Express. 3, (10 ), 2600 –2610 (2012). 2156-7085 CrossRef
Huang  Y.  et al., “Microvascular anastomosis guidance and evaluation using real-time three-dimensional Fourier-domain Doppler optical coherence tomography,” J. Biomed. Opt.. 18, (11 ), 111404  (2013). 1083-3668 CrossRef
Wang  Y.  et al., “Retinal blood flow measurement by circumpapillary Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt.. 13, (6 ), 064003  (2008). 1083-3668 CrossRef
Doblhoff-Dier  V.  et al., “Measurement of the total retinal blood flow using dual beam Fourier-domain Doppler optical coherence tomography with orthogonal detection planes,” Biomed. Opt. Express. 5, (2 ), 630 –642 (2014). 2156-7085 CrossRef
Srinivasan  V. J.  et al., “Quantitative cerebral blood flow with optical coherence tomography,” Opt. Express. 18, (3 ), 2477 –2494 (2010). 1094-4087 CrossRef
Leitgeb  R. A.  et al., “Doppler optical coherence tomography,” Prog. Retin. Eye Res.. 41, , 26 –43 (2014).CrossRef
Goldstein  R. M., , Zebker  H. A., and Werner  C. L., “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci.. 23, (4 ), 713 –720 (1988).CrossRef
Bone  D. J., “Fourier fringe analysis: the two-dimensional phase unwrapping problem,” Appl. Opt.. 30, (25 ), 3627 –3323 (1991). 0003-6935 CrossRef
Quiroga  J. A., , Gonzalez-Cano  A., and Bernabeu  E., “Phase unwrapping algorithm based on adaptive criterion,” Appl. Opt.. 34, (14 ), 2560 –2563 (1995). 0003-6935 CrossRef
Pritt  M. D., and Shipman  J. S., “Least-squares two-dimensional phase unwrapping using FFT’s,” IEEE Trans. Geosci. Remote Sens.. 32, (3 ), 706 –708 (1994). 0196-2892 CrossRef
Ghiglia  D. C., and Romero  L. A., “Robust two-dimensional weighted and un-weighted phase unwrapping that uses fast transforms and iterative methods,” J. Opt. Soc. Am. A.. 11, (1 ), 107 –117 (1994). 0740-3232 CrossRef
Costantini  M., “A novel phase unwrapping method based on network programming,” IEEE Trans. Geosci. Remote Sens.. 36, (3 ), 813 –821 (1998). 0196-2892 CrossRef
Gass  J., , Dakoff  A., and Kim  M. K., “Phase imaging without 2π ambiguity by multiple-wavelength digital holography,” Opt. Lett.. 28, (13 ), 1141 –1143 (2003). 0146-9592 CrossRef
Parshall  D., and Kim  M. K., “Digital holographic microscopy with dual wavelength phase unwrapping,” Appl. Opt.. 45, (3 ), 451 –459 (2006). 0003-6935 CrossRef
Khmaladze  A.  et al., “Dual-wavelength linear regression phase unwrapping in three-dimensional microscopic images of cancer cells,” Opt. Lett.. 36, (6 ), 912 –914 (2011). 0146-9592 CrossRef
Hendargo  H. C.  et al., “Synthetic wavelength based phase unwrapping in spectral domain optical coherence tomography,” Opt. Express. 17, (7 ), 5039 –5051 (2009). 1094-4087 CrossRef
Xia  X. G., and Wang  G., “Phase unwrapping and a robust Chinese remainder theorem,” IEEE Signal Process. Lett.. 14, (4 ), 247 –250 (2007). 1070-9908 CrossRef
Tang  S., , Zhang  X., and Tu  D., “Micro-phase measuring profilometry: its sensitivity analysis and phase unwrapping,” Opt. Lasers Eng.. 72, , 47 –57 (2015).CrossRef
Karout  S. A.  et al., “Two-dimensional phase unwrapping using a hybrid genetic algorithm,” Appl. Opt.. 46, (5 ), 730 –743 (2007). 0003-6935 CrossRef
Ghiglia  D. C., , Mastin  G. A., and Romero  L. A., “Cellular-automata method for phase unwrapping,” J. Opt. Soc. Am. A. 4, (1 ), 267 –280 (1987). 0740-3232 CrossRef
da Silva Maciel  L., and Albertazzi  A. G., “Swarm-based algorithm for phase unwrapping,” Appl. Opt.. 53, (24 ), 5502 –5509 (2014). 0003-6935 CrossRef
Huang  H. Y.  et al., “Path-independent phase unwrapping using phase gradient and total-variation (TV) denoising,” Opt. Express. 20, (13 ), 14075 –14089 (2012). 1094-4087 CrossRef
Loffeld  O.  et al., “Phase unwrapping for SAR interferometry—a data fusion approach by Kalman filtering,” IEEE Trans. Geosci. Remote Sens.. 46, (1 ), 47 –58 (2008). 0196-2892 CrossRef
Xie  X. M., and Pi  Y. M., “Multi-baseline phase unwrapping algorithm based on the unscented Kalman filter,” IET Radar Sonar Navig.. 5, , 296 –304 (2011).CrossRef
Xie  X., and Li  Y., “Enhanced phase unwrapping algorithm based on unscented Kalman filter, enhanced phase gradient estimator, and path-following strategy,” Appl. Opt.. 53, (18 ), 4049 –4060 (2014). 0003-6935 CrossRef
Xia  S.  et al., “Adaptive anisotropic diffusion for noise reduction of phase images in Fourier domain Doppler optical coherence tomography,” Biomed. Opt. Express. 7, (8 ), 2912 –2926 (2016). 2156-7085 CrossRef
Huang  Y.  et al., “MEMS-based handheld Fourier domain Doppler optical coherence tomography for intraoperative microvascular anastomosis imaging,” PLoS One. 9, (12 ), e114215  (2014). 1932-6203 CrossRef

Shaoyan Xia received his BS degree from Binzhou University, China, in 2014. He is currently a master’s student at the School of Optoelectronics at Beijing Institute of Technology, China. His research interests mainly focus on image processing.

Yong Huang received his PhD from the Department of Electrical and Computer Engineering from Johns Hopkins University, USA, in 2013. He is currently an associate professor at the School of Optoelectronics at Beijing Institute of Technology, China. His research interests include intraoperative imaging using optical coherence tomography, biomedical imaging, and image processing.

Shizhao Peng received his BS degree from Beijing Institute of Technology in 2015. Currently, he is a master’s student at the School of Optoelectronics at Beijing Institute of Technology, China.

Yanfeng Wu received his BS degree from Jilin University, China, in 2015. He is currently a master’s student at the School of Optoelectronics at Beijing Institute of Technology, China.

Xiaodi Tan received his PhD from Tokyo University in 2001. He is a professor at the School of Optoelectronics at Beijing Institute of Technology, China. His research interests include holographic data storage, holographic display, and biophotonics.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation

Shaoyan Xia ; Yong Huang ; Shizhao Peng ; Yanfeng Wu and Xiaodi Tan
"Robust phase unwrapping for phase images in Fourier domain Doppler optical coherence tomography", J. Biomed. Opt. 22(3), 036014 (Mar 29, 2017). ; http://dx.doi.org/10.1117/1.JBO.22.3.036014


Figures

Graphic Jump Location
Fig. 1
F1 :

Gradient analysis of simulated phase image with no noise. (a) True phase image, (a-1) gradient of true phase image in the x direction, (a-2) gradient of true phase image in the y direction, (b) wrapped phase image, (b-1) gradient of wrapped phase image in the x direction, (b-2) gradient of wrapped phase image in the y direction, (b-3) gradient in the x direction after unwrapping operator, (b-4) gradient in the y direction after unwrapping operator, (c) gradient along the y direction of row 20 in (a-2) and (b-2). Red line denotes true phase image gradient and black line denotes wrapped phase image gradient with phase jumps. Black arrow and red arrow denote the phase jump and unwrapping operation process.

Graphic Jump Location
Fig. 2
F2 :

Unwrapping results of network programming and our proposed method for the simulated image with different random-valued impulse noise levels: (a) true phase image, (b) true phase image with noise, (c) wrapped phase image, (d) unwrapped results using network programming method, and (e) unwrapped results using our proposed method.

Graphic Jump Location
Fig. 3
F3 :

Comparison of (a) RMSE1, (b) NAD, and (c) average execution time of 10 runs using network programming and our proposed methods for images with different random-valued impulse noise levels.

Graphic Jump Location
Fig. 4
F4 :

Phase unwrapping result comparison for phantom images with different flux levels (251×841  pixels). (a) Intensity image, (b) phase-resolved Doppler image, (c) unwrapped phase images with network programming method, and (d) unwrapped phase images with our proposed method. (Scale bar: 250  μm.)

Graphic Jump Location
Fig. 7
F7 :

Comparison of (a) RMSE2 and (b) average execution time of 10 runs using network programming and our proposed method for vascular images with different wrapping degrees.

Graphic Jump Location
Fig. 6
F6 :

Phase unwrapping results using network programming and our proposed method for vessel phase images with different wrapping degrees (200×999  pixels). (a) Vessel intensity image, (b) phase-resolved Doppler image, (c) unwrapped phase image using network programming method, and (d) unwrapped phase image using our proposed method. (Scale bar: 250  μm.)

Graphic Jump Location
Fig. 5
F5 :

Comparison of (a) RMSE2 and (b) average execution time of 10 runs using network programming and our proposed method for phantom images with different flux levels.

Tables

References

Zebker  H. A., and Lu  Y., “Phase unwrapping algorithms for radar interferometry: residue-cut, least-squares, and synthesis algorithms,” J. Opt. Soc. Am. A. 15, (3 ), 586 –598 (1998). 0740-3232 CrossRef
Goldstein  R. M., , Zebker  H. A., and Werner  C. L., “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci.. 23, (4 ), 713 –720 (1988).CrossRef
Bagher-Ebadian  H., , Jiang  Q., and Ewing  J. R., “A modified Fourier-based phase unwrapping algorithm with an application to MRI venography,” J. Magn. Reson. Imaging. 27, (3 ), 649 –652 (2008).CrossRef
Chavez  S., , Xiang  Q. S., and An  L., “Understanding phase maps in MRI: a new cutline phase unwrapping method,” IEEE Trans. Med. Imag.. 21, (8 ), 966 –977 (2002). 0278-0062 CrossRef
Li  S.  et al., “Two-dimensional wavelet transform for reliability-guided phase unwrapping in optical fringe pattern analysis,” Appl. Opt.. 51, (12 ), 2026 –2034 (2012). 0003-6935 CrossRef
Xu  Y.  et al., “Recovery of absolute height from wrapped phase maps for fringe projection profilometry,” Opt. Express. 22, (14 ), 16819 –16828 (2014). 1094-4087 CrossRef
Perciante  C. D.  et al., “Wrapping-free phase retrieval with applications to interferometry, 3D-shape profiling, and deflectometry,” Appl. Opt.. 54, (10 ), 3018 –3023 (2015) 0003-6935 CrossRef
Cheng  Z.  et al., “Practical phase unwrapping of interferometric fringes based on unscented Kalman filter technique,” Opt. Express. 23, (25 ), 32337 –32349 (2015). 1094-4087 CrossRef
Ma  L. H.  et al., “Fast algorithm for reliability-guided phase unwrapping in digital holographic microscopy,” Appl. Opt.. 51, (36 ), 8800 –8807 (2012). 0003-6935 CrossRef
Iglesias  I., “Phase estimation from digital holograms without unwrapping,” Opt. Express. 22, (18 ), 21340 –21346 (2014). 1094-4087 CrossRef
Khodadad  D., “Phase-derivative-based estimation of a digital reference wave from a single off-axis digital hologram,” Appl. Opt.. 55, (7 ), 1663 –1669 (2016). 0003-6935 CrossRef
Jaedicke  V.  et al., “Multiwavelength phase unwrapping and aberration correction using depth filtered digital holography,” Opt. Lett.. 39, (14 ), 4160 –4163 (2014). 0146-9592 CrossRef
van den Doel  L. R., and van Vliet  L. J., “Temporal phase-unwrapping algorithm for dynamic interference pattern analysis in interference-contrast microscopy,” Appl. Opt.. 40, (25 ), 4487 –4500 (2001). 0003-6935 CrossRef
Yaqoob  Z.  et al., “Improved phase sensitivity in spectral domain phase microscopy using line-field illumination and self phase-referencing,” Opt. Express. 17, (13 ), 10681 –10687 (2009). 1094-4087 CrossRef
Watanabe  E., , Hoshiba  T., and Javidi  B., “High-precision microscopic phase imaging without phase unwrapping for cancer cell identification,” Opt. Lett.. 38, (8 ), 1319 –1321 (2013). 0146-9592 CrossRef
Morofke  D.  et al., “Wide dynamic range detection of bidirectional flow in Doppler optical coherence tomography using a two-dimensional Kasai estimator,” Opt. Lett.. 32, (3 ), 253 –255 (2007). 0146-9592 CrossRef
Hendargo  H. C.  et al., “Doppler velocity detection limitations in spectrometer-based versus swept-source optical coherence tomography,” Biomed. Opt. Express. 2, (8 ), 2175 –2188 (2011). 2156-7085 CrossRef
Sun  C.  et al., “In vivo feasibility of endovascular Doppler optical coherence tomography,” Biomed. Opt. Express. 3, (10 ), 2600 –2610 (2012). 2156-7085 CrossRef
Huang  Y.  et al., “Microvascular anastomosis guidance and evaluation using real-time three-dimensional Fourier-domain Doppler optical coherence tomography,” J. Biomed. Opt.. 18, (11 ), 111404  (2013). 1083-3668 CrossRef
Wang  Y.  et al., “Retinal blood flow measurement by circumpapillary Fourier domain Doppler optical coherence tomography,” J. Biomed. Opt.. 13, (6 ), 064003  (2008). 1083-3668 CrossRef
Doblhoff-Dier  V.  et al., “Measurement of the total retinal blood flow using dual beam Fourier-domain Doppler optical coherence tomography with orthogonal detection planes,” Biomed. Opt. Express. 5, (2 ), 630 –642 (2014). 2156-7085 CrossRef
Srinivasan  V. J.  et al., “Quantitative cerebral blood flow with optical coherence tomography,” Opt. Express. 18, (3 ), 2477 –2494 (2010). 1094-4087 CrossRef
Leitgeb  R. A.  et al., “Doppler optical coherence tomography,” Prog. Retin. Eye Res.. 41, , 26 –43 (2014).CrossRef
Goldstein  R. M., , Zebker  H. A., and Werner  C. L., “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci.. 23, (4 ), 713 –720 (1988).CrossRef
Bone  D. J., “Fourier fringe analysis: the two-dimensional phase unwrapping problem,” Appl. Opt.. 30, (25 ), 3627 –3323 (1991). 0003-6935 CrossRef
Quiroga  J. A., , Gonzalez-Cano  A., and Bernabeu  E., “Phase unwrapping algorithm based on adaptive criterion,” Appl. Opt.. 34, (14 ), 2560 –2563 (1995). 0003-6935 CrossRef
Pritt  M. D., and Shipman  J. S., “Least-squares two-dimensional phase unwrapping using FFT’s,” IEEE Trans. Geosci. Remote Sens.. 32, (3 ), 706 –708 (1994). 0196-2892 CrossRef
Ghiglia  D. C., and Romero  L. A., “Robust two-dimensional weighted and un-weighted phase unwrapping that uses fast transforms and iterative methods,” J. Opt. Soc. Am. A.. 11, (1 ), 107 –117 (1994). 0740-3232 CrossRef
Costantini  M., “A novel phase unwrapping method based on network programming,” IEEE Trans. Geosci. Remote Sens.. 36, (3 ), 813 –821 (1998). 0196-2892 CrossRef
Gass  J., , Dakoff  A., and Kim  M. K., “Phase imaging without 2π ambiguity by multiple-wavelength digital holography,” Opt. Lett.. 28, (13 ), 1141 –1143 (2003). 0146-9592 CrossRef
Parshall  D., and Kim  M. K., “Digital holographic microscopy with dual wavelength phase unwrapping,” Appl. Opt.. 45, (3 ), 451 –459 (2006). 0003-6935 CrossRef
Khmaladze  A.  et al., “Dual-wavelength linear regression phase unwrapping in three-dimensional microscopic images of cancer cells,” Opt. Lett.. 36, (6 ), 912 –914 (2011). 0146-9592 CrossRef
Hendargo  H. C.  et al., “Synthetic wavelength based phase unwrapping in spectral domain optical coherence tomography,” Opt. Express. 17, (7 ), 5039 –5051 (2009). 1094-4087 CrossRef
Xia  X. G., and Wang  G., “Phase unwrapping and a robust Chinese remainder theorem,” IEEE Signal Process. Lett.. 14, (4 ), 247 –250 (2007). 1070-9908 CrossRef
Tang  S., , Zhang  X., and Tu  D., “Micro-phase measuring profilometry: its sensitivity analysis and phase unwrapping,” Opt. Lasers Eng.. 72, , 47 –57 (2015).CrossRef
Karout  S. A.  et al., “Two-dimensional phase unwrapping using a hybrid genetic algorithm,” Appl. Opt.. 46, (5 ), 730 –743 (2007). 0003-6935 CrossRef
Ghiglia  D. C., , Mastin  G. A., and Romero  L. A., “Cellular-automata method for phase unwrapping,” J. Opt. Soc. Am. A. 4, (1 ), 267 –280 (1987). 0740-3232 CrossRef
da Silva Maciel  L., and Albertazzi  A. G., “Swarm-based algorithm for phase unwrapping,” Appl. Opt.. 53, (24 ), 5502 –5509 (2014). 0003-6935 CrossRef
Huang  H. Y.  et al., “Path-independent phase unwrapping using phase gradient and total-variation (TV) denoising,” Opt. Express. 20, (13 ), 14075 –14089 (2012). 1094-4087 CrossRef
Loffeld  O.  et al., “Phase unwrapping for SAR interferometry—a data fusion approach by Kalman filtering,” IEEE Trans. Geosci. Remote Sens.. 46, (1 ), 47 –58 (2008). 0196-2892 CrossRef
Xie  X. M., and Pi  Y. M., “Multi-baseline phase unwrapping algorithm based on the unscented Kalman filter,” IET Radar Sonar Navig.. 5, , 296 –304 (2011).CrossRef
Xie  X., and Li  Y., “Enhanced phase unwrapping algorithm based on unscented Kalman filter, enhanced phase gradient estimator, and path-following strategy,” Appl. Opt.. 53, (18 ), 4049 –4060 (2014). 0003-6935 CrossRef
Xia  S.  et al., “Adaptive anisotropic diffusion for noise reduction of phase images in Fourier domain Doppler optical coherence tomography,” Biomed. Opt. Express. 7, (8 ), 2912 –2926 (2016). 2156-7085 CrossRef
Huang  Y.  et al., “MEMS-based handheld Fourier domain Doppler optical coherence tomography for intraoperative microvascular anastomosis imaging,” PLoS One. 9, (12 ), e114215  (2014). 1932-6203 CrossRef

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

Advertisement
  • Don't have an account?
  • Subscribe to the SPIE Digital Library
  • Create a FREE account to sign up for Digital Library content alerts and gain access to institutional subscriptions remotely.
Access This Article
Sign in or Create a personal account to Buy this article ($20 for members, $25 for non-members).
Access This Proceeding
Sign in or Create a personal account to Buy this article ($15 for members, $18 for non-members).
Access This Chapter

Access to SPIE eBooks is limited to subscribing institutions and is not available as part of a personal subscription. Print or electronic versions of individual SPIE books may be purchased via SPIE.org.