*π*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.

## 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^{24} 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^{27} 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 Costantini^{29} 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,^{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\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

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)

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

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

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

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

where ${k}_{2}(i,j)$ is an integer.With the simplified structure, Eq. (6) can be written as

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

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

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

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.

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

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

## Eq. (15)

$${e}_{x}(i,j)=\mathrm{\Delta}{\varphi}_{x}(i,j)-\square {\phi}_{x}(i,j),\phantom{\rule{0ex}{0ex}}{e}_{y}(i,j)=\mathrm{\Delta}{\varphi}_{y}(i,j)-\square {\phi}_{y}(i,j).$$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)

$${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)]\phantom{\rule{0ex}{0ex}}-[\square {\phi}_{x}(i,j+1)-\square {\phi}_{x}(i,j)-\square {\phi}_{y}(i+1,j)+\square {\phi}_{y}(i,j)].$$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

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)

$$\mathrm{min}\sum _{i,j}{c}_{x}(i,j)|{e}_{x}(i,j)|+\sum _{i,j}{c}_{y}(i,j)|{e}_{y}(i,j)|\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}\text{\hspace{0.17em}\hspace{0.17em}}{e}_{x}(i,j+1)-{e}_{x}(i,j)-{e}_{y}(i+1,j)+{e}_{y}(i,j)=R(i,j),$$^{29}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.

## 2.2.

### Unwrapped Phase Restoration

In Sec. 2.1, every pixel can be affected by error amounts ${e}_{x},{e}_{y}$, so estimation of the gradient of true phase image $\mathrm{\Delta}{\varphi}_{x},\mathrm{\Delta}{\varphi}_{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:

## Eq. (19)

$$\varphi (i,j)=\phi (i,j)+\sum _{j}\mathrm{\Delta}{\varphi}_{1}(1,j)+\sum _{\mathrm{i}}\mathrm{\Delta}{\varphi}_{2}(i,j),$$## Eq. (20)

$$\varphi (i,j)=\phi (i,j)+\sum _{i}\mathrm{\Delta}{\varphi}_{2}(i,1)+\sum _{j}\mathrm{\Delta}{\varphi}_{1}(i,j).$$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.

## 2.3.

### Evaluation Metrics for Phase Unwrapping

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

## 2.3.1.

#### Root-mean-square-error

It can be used to describe the closeness of unwrapped phase image $I$ to the true phase image $\varphi $.

## Eq. (21)

$$\mathrm{RMSE}1=\sqrt{\frac{\sum _{\mathrm{i}=1}^{m}\sum _{\mathrm{j}=1}^{n}{({I}_{i,j}-{\varphi}_{i,j})}^{2}}{m\times n}},$$## 2.3.2.

#### Noise amplification degree

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

## Eq. (22)

$$\mathrm{NAD}=\frac{\sum _{\mathrm{i}=1}^{m}\sum _{\mathrm{j}=1}^{n}|{I}_{i,j}-{\varphi}_{i,j}|}{\sum _{\mathrm{i}=1}^{m}\sum _{\mathrm{j}=1}^{n}|{\varphi}_{i,j}-{L}_{i,j}|},$$## 2.3.3.

#### 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 ($-\pi ,\pi $] 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

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

## Eq. (24)

$$\mathrm{RMSE}2=\sqrt{\frac{\sum _{\mathrm{i}=1}^{m}\sum _{j=1}^{n}{E}_{i,j}^{2}}{m\times n}}.$$In cases where we do not have access to true phase image, we can use RMSE2 to evaluate the unwrapping performance.

## 2.3.4.

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

## 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.^{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\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.

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.

## 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\text{\hspace{0.17em}\hspace{0.17em}}\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.

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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{l}/\mathrm{min}$ and $Q=1150\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{l}/\mathrm{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.

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

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.

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

## Disclosures

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

## Acknowledgments

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.

## References

## Biography

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