1.

Introduction

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,^5–8 digital holography,^9–12 microscopic imaging,^13–15 and phase-resolved Doppler optical coherence tomography (DOCT).^16–18

DOCT is a noninvasive high-resolution three-dimensional imaging modality and has been widely used for various flow studies.^19–23 However, phase-resolved Doppler images often suffer from wrapping issue that limits phase values in the range of ($-\pi ,\pi$]. Deviation from real phase profile will occur if measured phase range is over $2\pi$. 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,^24–26 minimum-norm method,^27–29 multiwavelength method,^30–33 and other algorithms.^34–39 Usually preprocess and in-process noise reduction techniques are combined with these methods.^40–42

For path-following method, Goldstein et al.’s²⁴ 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-squares²⁷ 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.²⁸ In the weighted minimum-norm methods, network flow algorithm is the most attractive, which was proposed by Costantini²⁹ 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\pi$. 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.^30–33 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 microphase^34,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,³⁶ cellular-automata,³⁷ and swarm intelligence,³⁸ are introduced for unwrapping optimization. Path-independent phase unwrapping using total-variation denoising adopts noise reduction techniques during the phase unwrapping process.³⁹

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\pi$. 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.

2.

Method

2.1.

Acquisition of Unwrapped Phase Gradients

Due to the $2\pi$ ambiguity of phase-resolved Doppler image, the relationship between true phase image $\varphi$ and wrapped phase image $\phi$ for every pixel can be expressed as

Eq. (1)

$$\phi (i,j)=W[\varphi (i,j)]=\varphi (i,j)+2\pi k(i,j)(i=1,\dots ,m;j=1,\dots ,n),$$

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 $\mathrm{mod}(\varphi +\pi ,2\pi )-\pi$ (mod is modulo operation) to realize the effect of wrapping operator $W$. After wrapping operation, phase is limited to the range ($-\pi ,\pi$]. 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\pi$ 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)

Eq. (2)

$$G[\phi (i,j)]=G[\varphi (i,j)+2\pi k(i,j)].$$

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

Eq. (3)

$$G[\phi (i,j)]=G[\varphi (i,j)]+G[2\pi k(i,j)].$$

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

Eq. (4)

$$G[\phi (i,j)]=G[\varphi (i,j)]+2\pi k_1(i,j),$$

where $k_1(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.)

Eq. (5)

$$U\{G[\phi (i,j)]\}=U\{G[\varphi (i,j)]+2\pi k_1(i,j)\}.$$

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

Eq. (6)

$$U\{G[\phi (i,j)]\}=\{G[\varphi (i,j)]+2\pi k_1(i,j)\}+2\pi k_2(i,j),$$

where $k_2(i,j)$ is an integer.

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

Eq. (7)

$$U\{G[\phi (i,j)]\}=G[\varphi (i,j)]+2\pi [k_1(i,j)+k_2(i,j)].$$

In Eq. (7), the value range of $U\{G[\phi (i,j)]\}$ is ($-\pi ,\pi$] due to the use of unwrapping operator $U$. In general, true phase image $\varphi$ is continuous, thus the value range of $G[\varphi (i,j)]$ is also ($-\pi ,\pi$]. So, to keep the same value range for both sides of Eq. (7), we can conclude that

Eq. (8)

$$k_1(i,j)+k_2(i,j)=0,\forall (i,j).$$

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

Eq. (9)

$$G[\varphi (i,j)]=U\{G[\phi (i,j)]\}.$$

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 $\square {\phi }_x(i,j)=\mathrm{\Delta }{\varphi }_x(i,j)$, $\square {\phi }_y(i,j)=\mathrm{\Delta }{\varphi }_y(i,j)$, where $\square {\phi }_x,\square {\phi }_y$ denotes the derivative of wrapped phase image after using unwrapping operator $U$ and $\mathrm{\Delta }{\varphi }_x,\mathrm{\Delta }{\varphi }_y$ denotes the derivative of true phase image. Specifically, their derivatives are given by

Eq. (10)

$$\square {\phi }_x(i,j)=U[\phi (i+1,j)-\phi (i,j)],$$

Eq. (11)

$$\square {\phi }_y(i,j)=U[\phi (i,j+1)-\phi (i,j)],$$

Eq. (12)

$$\mathrm{\Delta }{\varphi }_x(i,j)=\varphi (i+1,j)-\varphi (i,j),$$

Eq. (13)

$$\mathrm{\Delta }{\varphi }_y(i,j)=\varphi (i,j+1)-\varphi (i,j).$$

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.

Fig. 1

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

Eq. (14)

$$G[\varphi (i,j)]=U\{G[\phi (i,j)]\}+e(i,j).$$

So the error amounts between $\square {\phi }_{x,y}$ and $\mathrm{\Delta }{\varphi }_{x,y}$ can be expressed as

Eq. (15)

$$\begin{gathered} e_x(i,j)=\mathrm{\Delta }{\varphi }_x(i,j)-\square {\phi }_x(i,j), \\ {e}_y(i,j)=\mathrm{\Delta }{\varphi }_y(i,j)-\square {\phi }_y(i,j). \end{gathered}$$

Due to the uncertain property of noise, error amounts $e_x,e_y$ can be any values, instead of integer multiple of $2\pi$. 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

Eq. (16)

$$\begin{gathered} e_x(i,j+1)-e_x(i,j)-e_y(i+1,j)+e_y(i,j)=[\mathrm{\Delta }{\varphi }_x(i,j+1)-\mathrm{\Delta }{\varphi }_x(i,j)-\mathrm{\Delta }{\varphi }_y(i+1,j)+\mathrm{\Delta }{\varphi }_y(i,j)] \\ -[\square {\phi }_x(i,j+1)-\square {\phi }_x(i,j)-\square {\phi }_y(i+1,j)+\square {\phi }_y(i,j)]. \end{gathered}$$

From Eqs. (12) and (13), we can calculate that both $\mathrm{\Delta }{\varphi }_x(i,j+1)-\mathrm{\Delta }{\varphi }_x(i,j)$ and $\mathrm{\Delta }{\varphi }_y(i+1,j)-\mathrm{\Delta }{\varphi }_y(i,j)$ are equal to $\varphi (i+1,j+1)-\varphi (i+1,j)-\varphi (i,j+1)+\varphi (i,j)$. Thus $\mathrm{\Delta }{\varphi }_x(i,j+1)-\mathrm{\Delta }{\varphi }_x(i,j)-\mathrm{\Delta }{\varphi }_y(i+1,j)+\mathrm{\Delta }{\varphi }_y(i,j)$ equals 0 for every pixel. Hence, Eq. (16) can be simplified as

Eq. (17)

$$e_x(i,j+1)-e_x(i,j)-e_y(i+1,j)+e_y(i,j)=R(i,j),$$

where $R(i,j)$ is $-[\square {\phi }_x(i,j+1)-\square {\phi }_x(i,j)-\square {\phi }_y(i+1,j)+\square {\phi }_y(i,j)]$.

As the error values $e_x,e_y$ is limited, thus process of solving $e_x,e_y$ can be treated as the following optimization problem:

Eq. (18)

$$\begin{gathered} \min \sum _{i,j}{c}_x(i,j)|e_x(i,j)|+\sum _{i,j}{c}_y(i,j)|e_y(i,j)| \\ \mathrm{s}.\mathrm{t}{e}_x(i,j+1)-e_x(i,j)-e_y(i+1,j)+e_y(i,j)=R(i,j), \end{gathered}$$

where $c_x,c_y$ are weighted coefficients to balance the objective function,²⁹ and $R(i,j)$ is $-[\square {\phi }_x(i,j+1)-\square {\phi }_x(i,j)-\square {\phi }_y(i+1,j)+\square {\phi }_y(i,j)]$. To solve this problem, we can turn Eq. (18) into linear programming problem and use certain optimizing tool to obtain $e_x,e_y$ values. Equation (18) uses the error terms $e_x(i,j),e_y(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.

3.

Results of Simulated Images

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.⁴³ 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\pi$. 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.

Fig. 2

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.

Fig. 3

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.

4.

Results of Phase-Resolved Doppler OCT Images

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.

4.1.

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\mu \mathrm{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.

Fig. 4

Phase unwrapping result comparison for phantom images with different flux levels ($251\times 841\text{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\mu \mathrm{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\mu \mathrm{l}/\min$ and $Q=1150\mu \mathrm{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.

Fig. 5

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.

4.2.

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.

Fig. 6

Phase unwrapping results using network programming and our proposed method for vessel phase images with different wrapping degrees ($200\times 999\text{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\mu \mathrm{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.

Fig. 7

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.

5.

Conclusion

We have demonstrated the capability of our proposed method to solve the $2\pi$ 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.