A Robust Discontinuous Phase Unwrapping Based on Least-Squares Orientation Estimator

Weighted least-squares (WLS) phase unwrapping is widely used in optical engineering. However, this technique still has issues in coping with discontinuity as well as noise. In this paper, a new WLS phase unwrapping algorithm based on the least-squares orientation estimator (LSOE) is proposed to improve phase unwrapping robustness. Specifically, the proposed LSOE employs a quadratic error norm to constrain the distance between gradients and orientation vectors. The estimated orientation is then used to indicate the wrapped phase quality, which is in terms of a weight mask. The weight mask is calculated by post-processing, including a bilateral filter, STDS, and numerical relabeling. Simulation results show that the proposed method can work in a scenario in which the noise variance is 1.5. Comparisons with the four WLS phase unwrapping methods indicate that the proposed method provides the best accuracy in terms of segmentation mean error under the noisy patterns.

Spatial phase unwrapping (SPU) is a natural and straightforward technique to unwrap wrapped phase images. Generally, SPU methods can be classified into path-following and minimum-norm methods [2]. The basic principle of the path-following algorithm is to perform phase integration along the appropriate path to avoid errors caused by the unreliability of the wrapped phase data [2]. Typical path-following methods include Goldstein's branch-cut algorithm [5], quality-guided algorithms [6][7][8][9][10], Flynn's minimum discontinuity algorithm [11], and so on.
On the other hand, minimum-norm methods [2] have been proposed to unwrap phases by a global function-minimized phase unwrapping. In [2], by employing a fast Fourier and discrete cosine transform, a weighted least-squares method (WLS) was shown to solve phase unwrapping. Considering the effect of noise and discontinuity, the noisy and discontinuous regions are mapped to weight coefficients, and an effective weighted iteration is implemented to accurately solve the phase unwrapping problem. Obviously, the issue is to determine the weight mask, which is used to calculate the weight coefficients of the WLS phase unwrapping. To this end, various phase quality analysis methods, such as correlation maps [2], phase derivative variance (PDV) [12], derivative variance correlation maps (DVCM) [13], etc., have been proposed. It is worth mentioning that a regional reliability evaluation to define weight mask was introduced in [14]. The proposed approach is an acceptable solution for representing phase regions and abnormal edge quality, and it is suitable for specific optical phase-measurement systems with more wrapped phases.
It should be highlighted that the above-mentioned methods focus on continuous fringe patterns. Due to spatially independent surfaces, discontinuous patterns always happen, which yields an error in phase unwrapping. The situation is very challenging in case both noise and discontinuity exist. The difficulty of discontinuity identification is worth highlighting. Recently, several methods have been proposed to solve this problem. Cui et al. [15] introduced a new quality map by combining a phase Laplace gradient and a statistical method which transforms the phase gradient of each pixel into a phase Laplace gradient using its four neighboring pixels. Qian et al. [16] proposed a snake-assisted quality-guided phase unwrapping algorithm, which used a windowed Fourier transform (WFT)-based method to obtain an effective quality measure. It has been shown that this method can obtain good unwrapping results at discontinuities. Based on the wrapping-free phase retrieval method described in [17], Zhou et al. [18] proposed a binary mask of a reverse DVCM as the weight mask. Since the second derivative of the four-step phaseshifting fringe pattern needs to be calculated, these methods are only suitable for particular measurement objects with singularities. Concerning the discontinuous fringe pattern with noise, a local orientation coherence-based fringe segmentation method, which cooperates with boundary-aware coherence enhancing diffusion, was presented in [19]. Furthermore, a complete solution based on [19] was described in [20] to deal with discontinuous phase unwrapping problems. To solve the weak global noise problem, Yan et al. [21] proposed a WLS phase unwrapping algorithm based on a reliability mask. For different patterns, several trials were carried out to determine the appropriate threshold. Inspired by the local orientation coherence of [19], Li et al. [22] presented a least-squares phase unwrapping algorithm using orientation coherence as the weight mask. It is highlighted that Li provides the appropriate thresholds in his research [22].
To deal with the problem, i.e., the weight mask sensitivity to noise and heuristic selection of the threshold, a robust weight mask's calculation method based on the least-squares orientation estimator is presented in this paper. First, the quadratic error norm of the phase image between the gradient vector and the true orientation of each pixel is calculated. Next, a Gaussian kernel function is employed to smooth the incorrect gradients caused by noise. Then, the Lagrange multiplier is used to minimize the error to obtain an orientation map. After that, a bilateral filter [23] is adopted to remove the abnormal direction and preserve the discontinuous edges of the orientation map. Finally, the discontinuous and continuous orientation maps are processed by the SDTS method [24] and numerical relabeling. The experimental results show that the proposed algorithm is superior to the state-of-the-art methods. The contributions of the paper are summarized as follows: (1) A novel orientation estimator, LSOE, is proposed to construct a noise-robust orientation map; (2) A new weight mask solution is developed by post-processing, which comprises the bilateral filter, STDS method, and numerical relabeling instead of the traditional thresholding segmentation methods.
The remainder of this paper is organized as follows. The WLS method for phase unwrapping is briefly introduced in Section 2. In Section 3, we describe the proposed algorithm in detail. Section 4 demonstrates experimental results, followed by the discussion and conclusions in Sections 5 and 6.

WLS Phase Unwrapping
The relationship [2] between the wrapped phase ψ i,j and its unwrapped phase φ i,j on a discrete grid of points is expressed as where the integer multiples c i.j are the wrap counts, and Wrap(·) is the wrapping operator. i and j are the horizontal and vertical pixel coordinates, respectively. M and N are the horizontal and vertical sizes of the phase image, respectively.
To eliminate the influence of invalid pixels on discontinuous phase unwrapping, Ghiglia et al. [2] proposed a phase unwrapping method based on WLS. The aim of WLS phase unwrapping is to find the least-squares solution between the true phase difference and the wrapped phase difference ∆ y i,j in the horizontal and vertical directions as follows: where ∆ x i,j and ∆ y i,j are expressed as U i,j and V i,j are the weights coefficients, denoted as where w i,j is the weight mask. Generally, the weight map is a 0-1 binary image representing the continuous and discontinuous position, which is expressed as where θ i,j is the orientation map of the wrapped phase image, and r is the segmentation threshold. Equation (2) can be simplified to the following discrete Poisson equation with Neumann boundary conditions [2]: where q i,j is the weighted phase Laplacian [2]: The WLS phase unwrapping solution of Equation (7) should be obtained by a rapidly converging algorithm based on a preconditioned conjugate gradient. In the following, a robust orientation map to obtain a weight mask is proposed.

Least-Squares Orientation Estimator
To calculate the weight mask, Wang [25] used the orientation to represent phase quality. The technique to estimate orientation is to look at the structure information in a local neighborhood. Therefore, the structure tensor is used to extend the structure information of each pixel [26]. In Rivera's work [26], the structure tensor T p is used to compute local orientation: where r ∈ p, N denotes the number of pixels in the region p, and I i , I j denotes the horizontal and vertical partial derivatives of the wrapped phase image I, respectively. The eigenvector corresponding to the smaller eigenvalue of the structure tensor is estimated as the fringe orientation. Inspired by Rivera's orientation estimation method, a robust orientation estimator, LSOE, is introduced based on the quadratic error norm for each pixel based on the structure tensor. It should be emphasized that, ideally, the gradient vector g i,j of the fringe point is parallel to its orientation vector v i,j , as shown in Figure 1a, which is perpendicular to the fringe. The error between the gradient vector and the orientation vector is zero. We carefully choose the gradient vector g i,j as the first-order partial derivative of the wrapped phase image.
The eigenvector corresponding to the smaller eigenvalue of the structure tensor is estimated as the fringe orientation.
Inspired by Rivera's orientation estimation method, a robust orientation estimator, LSOE, is introduced based on the quadratic error norm for each pixel based on the structure tensor. It should be emphasized that, ideally, the gradient vector , i j g of the fringe point is parallel to its orientation vector , i j v , as shown in Figure 1a, which is perpendicular to the fringe. The error between the gradient vector and the orientation vector is zero. We carefully choose the gradient vector , i j g as the first-order partial derivative of the wrapped phase image. Now, we consider the situation in which an error exists between the gradient direction and the orientation, which is caused by the different deviation of the vertical direction in the fringe. The error , i j e of the gradient vector with respect to the orientation vector of the point is defined as As observed in Figure 1b   Now, we consider the situation in which an error exists between the gradient direction and the orientation, which is caused by the different deviation of the vertical direction in the fringe. The error e i,j of the gradient vector with respect to the orientation vector of the point is defined as As observed in Figure 1b, e i,j [27] is the projection of g i,j on the vector v i,j . The quadratic error norm measure ε i,j can be expressed as where K σ (·) is the Gaussian kernel function with mean zero and variance σ. Thus, ε can be written as To make g i,j and v i,j as parallel as possible, the quadratic error norm ε i,j should be minimized. The first item on the right side of Equation (11) is a calculated positive value that is not relevant to v i,j . Minimizing the error is equivalent to maximizing Lagrange multipliers are used to solve the maximum, subject to the constraint condition v i,j T v i,j = 1, as given below: where J i,j,σ = g i,j T g i,j K σ is the robust structure tensor, which is a nonlinear operator.
Equation (12) can be expressed as Thus, by differentiating the function and setting the derivative to zero, the result can be expressed as Similarly, the orientation vector v can be viewed as the eigenvector of the characteristic equation, i.e., Equation (14). The Lagrange multiplier λ is interpreted as the eigenvalue. If the phase pattern has a discontinuity, two unequal eigenvalues λ 1 and λ 2 are used to represent two different orientation regions and λ 1 > λ 2 > 0. Accordingly, the smaller value λ 2 is chosen to better indicate the special orientation of each point. Because J i,j,σ is a square matrix, the eigenvector v corresponding to the smaller eigenvalue λ 2 is unique and can be represented as The orientation θ i,j of the fringe pattern can be estimated as follows: The estimated orientation is mainly used to indicate the wrapped phase quality of the structural fringe pattern.

Weight Mask Calculation Based on Post-Processing
As previously mentioned in Section 3.1, the estimated orientation map θ i,j is served as the wrapped phase quality. To obtain the weight mask, it is necessary to further segment the orientation map for the WLS solution to unwrap phases. However, the hard threshold r, as shown in Equation (5), is difficult to determine due to the fact that the gap between two classes, i.e., the phase values on the continuous areas and the discontinuous areas, is small. Therefore, a new post-processing solution is proposed to obtain the weight masks without using a hard threshold as follows: where PF(·) is the post-processed function. The post-processing for obtaining the weight masks includes three steps: bilateral filtering, STDS, and numerical relabeling. In fact, the orientation jump between adjacent periods seriously affects the segmentation of the intrinsic boundary (discontinuity in the phase map), which results in incorrectly unwrapping the phase. To suppress this jumping problem, a bilateral filter is first used to eliminate the abnormal orientation. Then, considering the fact that the useful edges are buried in textures, we employ STDS with an image smoothing ability to extract/remove the edges/textures on the output of the bilateral filter, respectively. The outputs of the STDS are explicitly marked as two regions: regions with 0 at continuous points and regions with 0.9 at discontinuous points. Finally, the weight masks are obtained by relabeling two regions, i.e., the weight masks are 0, corresponding to continuous regions, and the weight masks are 1, corresponding to discontinuous regions. The whole scheme of the proposed method is demonstrated in Figure 2.
ering the fact that the useful edges are buried in textures, we employ STDS with an image smoothing ability to extract/remove the edges/textures on the output of the bilateral filter, respectively. The outputs of the STDS are explicitly marked as two regions: regions with 0 at continuous points and regions with 0.9 at discontinuous points. Finally, the weight masks are obtained by relabeling two regions, i.e., the weight masks are 0, corresponding to continuous regions, and the weight masks are 1, corresponding to discontinuous regions. The whole scheme of the proposed method is demonstrated in Figure 2.

Numerical Simulation
In this section, we report the results of a quantitative evaluation of the weight mask's segmentation performance of the proposed method and compare with the state-of-the-art methods, such as Li's method [22] (orientation coherence, OC), Lu's method [13] (derivative variance correlation map, DVCM), Cui's method [15] (phase Laplace derivative variance, PLDV), and Yan's method [21] (reliability mask, RM). All algorithms were implemented in MATLAB R2018b programming on a laptop equipped with an Intel (R) Core (TM) i5-8250U processor with a clock frequency of 1.6 GHz and 8 GB of RAM. The segmentation means error (SME) of the discontinuity boundaries between the measured and ground truth was used to quantitatively verify the accuracy of the weight mask calculation as below:  

Numerical Simulation
In this section, we report the results of a quantitative evaluation of the weight mask's segmentation performance of the proposed method and compare with the state-of-the-art methods, such as Li's method [22] (orientation coherence, OC), Lu's method [13] (derivative variance correlation map, DVCM), Cui's method [15] (phase Laplace derivative variance, PLDV), and Yan's method [21] (reliability mask, RM). All algorithms were implemented in MATLAB R2018b programming on a laptop equipped with an Intel (R) Core (TM) i5-8250U processor with a clock frequency of 1.6 GHz and 8 GB of RAM. The segmentation means error (SME) of the discontinuity boundaries between the measured and ground truth was used to quantitatively verify the accuracy of the weight mask calculation as below: where P * i is the ideal discontinuity position, and P i is the actual discontinuity coordinate. A fringe pattern (512 × 512) called LINE, containing both linear and irregular circular fringe segments, was simulated with a linear boundary, as shown in Figure 3a Table 1 (first row). The SMEs of the LINE pattern (variance 1.5) for the five algorithms are tabulated in Table 2 (first row). Another fringe pattern (512 × 512) called CIRCLE, containing both linear and regular circular fringe segments, was simulated with a circular boundary, as shown in Figure 4(a-1). Additive random noise of a mean of 0 and variance values of 0, 0.5, 0.8, and 1.5 was then added. The latter two are shown in Figure 4b-1,c-1, respectively. Figure 4a identified boundaries and red lines denoting the ground truth. Quantitatively, for fringe patterns with random noise having variances of 0, 0.8, and 1.5, the SMEs were 0.2072, 1.9633, and 5.0673, respectively. The SMEs of CIRCLE pattern (variance 1.5) for the five algorithms are shown in the second row of Table 2.

Experimental and Results Analysis
Experimental fringe patterns were tested to verify the weight mask's segmentation of the proposed method. A fringe pattern from the phase-shifting electronic speckle pattern shearing interferometer is shown in Figure 5a, where the upper-left corner of the fringe pattern is occupied with random noise. The identified discontinuous boundaries of Figure 5a are shown in Figure 5b-f.
Another deformation measurement experiment was conducted with the fringe pro-

Experimental and Results Analysis
Experimental fringe patterns were tested to verify the weight mask's segmentation of the proposed method. A fringe pattern from the phase-shifting electronic speckle   Table 3. Computational time of five algorithms for real scene (Figure 6a-1). For the wrapped phase with discontinuity, the accurate absolute phase distribution has two characteristics: (1) the phase has a jump at the discontinuity. (2) The wrapped phase after the jump is unwrapped independently of the pre-jump phase. Figure 7 shows the absolute phase of the 1-D wrapped phase (255th in Figure 6a-2) unwrapped by five algorithms. It is observed from Figure 7 that the proposed method has an obvious jump at mark point A/B, which enables the wrapped phase at mark point C to be unwrapped independently of the phase before the jump. However, other algorithms cannot satisfy the above two characteristics simultaneously.  Another deformation measurement experiment was conducted with the fringe projection profilometry based on spatial phase-shifting. The FPP system was composed of a DLP LightCrafter 4500 and JAI GO-5000M-USB camera of Japan to verify the performance of the proposed algorithm. To evaluate the unwrapping phase of the proposed method, in this study, we compare it with other methods, such as WLS-OC [22], WLS-DVCM [13], WLS-PLDV [15], and WLS-RM [21].

WLS-OC WLS-DVCM WLS-PLDV WLS-RM Ours
A real scenario containing a box and a basketball, as shown in Figure 6a-1, was tested in the experiment. It is seen that the scene contains horizontal, vertical, and arcshaped discontinuity. To clearly illustrate the noisy wrapped phase image, we provide the wrapped phase image in Figure 6a Table 3. Computational time of five algorithms for real scene (Figure 6a-1). For the wrapped phase with discontinuity, the accurate absolute phase distribution has two characteristics: (1) the phase has a jump at the discontinuity. (2) The wrapped phase after the jump is unwrapped independently of the pre-jump phase. Figure 7 shows the absolute phase of the 1-D wrapped phase (255th in Figure 6a-2) unwrapped by five algorithms. It is observed from Figure 7 that the proposed method has an obvious jump at mark point A/B, which enables the wrapped phase at mark point C to be unwrapped independently of the phase before the jump. However, other algorithms cannot satisfy the above two characteristics simultaneously. The computational times of the five algorithms for the real scene are tabulated in Table 3. It is seen that the proposed method ranks second. The computation mainly costs on the STDS operation, which employs an iterative strategy for segmentation. The high segmentation performance of the STDS method is achieved at the expense of computational complexity. Table 3. Computational time of five algorithms for real scene (Figure 6a-1). For the wrapped phase with discontinuity, the accurate absolute phase distribution has two characteristics: (1) the phase has a jump at the discontinuity. (2) The wrapped phase after the jump is unwrapped independently of the pre-jump phase. Figure 7 shows the absolute phase of the 1-D wrapped phase (255th in Figure 6a-2) unwrapped by five algorithms. It is observed from Figure 7 that the proposed method has an obvious jump at mark point A/B, which enables the wrapped phase at mark point C to be unwrapped independently of the phase before the jump. However, other algorithms cannot satisfy the above two characteristics simultaneously. Table 3. Computational time of five algorithms for real scene (Figure 6a-1 For the wrapped phase with discontinuity, the accurate absolute phase distr has two characteristics: (1) the phase has a jump at the discontinuity. (2) The w phase after the jump is unwrapped independently of the pre-jump phase. Figure 7 the absolute phase of the 1-D wrapped phase (255th in Figure 6a-2) unwrapped algorithms. It is observed from Figure 7 that the proposed method has an obviou at mark point A/B, which enables the wrapped phase at mark point C to be unw independently of the phase before the jump. However, other algorithms cannot sat above two characteristics simultaneously.

Discussion
In this work, a discontinuous phase unwrapping method is proposed by using LSOE. The numerical results in Figures 3 and 4 show that complete discontinuous segmentation can be achieved even when the noise variance is 1.5. It is observed from Table 2 and Figure 7 that the proposed method achieves the best performance in terms of SME. Moreover, the proposed method is robust due to the following facts: (1) The orientation jump between adjacent periods seriously affects the segmentation of the intrinsic boundaries (discontinuity in the phase map), which yields an incorrect unwrapping phase. To suppress this jumping problem, a bilateral filter in post-processing is first used to eliminate the abnormal orientation.
(2) Considering the fact that the useful edges are buried in textures, the STDS with an image smoothing capability is adopted to extract/remove the edges/textures on the output of the bilateral filter. The experimental results demonstrate that the proposed method is robust to noise on the simulated patterns with different variance noise. It is also shown in Figure 5 that the noise robustness of the proposed method is still valid in real patterns. It is noted that the accuracy of discontinuity segmentation decreases if noise increases, as shown in Figures 3 and 4. Table 3 shows the computational time of the five methods, from which it can be found that the proposed method takes 12 s and achieves the second rank in performance. Improving the algorithm and coding in C++ will be our future work so that the method can be applied to real-time applications.

Conclusions
It is known that the WLS phase unwrapping method depends on the accuracy of segmentation for the weight mask. To this end, this paper proposes a robust weight mask segmentation method based on LSOE. The proposed method employs a quadratic error norm to constrain the distance between gradients and orientation vectors, which minimizes the noise-reduced orientation error. The weight mask is derived from the LSOE by further performing post-processing, which includes a bilateral filter, STDS, and numerical relabeling. The numerical results in Figures 3 and 4 show that complete discontinuous segmentation can be achieved even when the noise variance is 1.5. It is also demonstrated that the noise robustness of the proposed method is still valid in real patterns. The comparisons with other methods, including Li's OC, Lu's DVCM, Cui's PLDV, and Yan's RM, for WLS phase unwrapping indicate that the proposed method offers the best performance in terms of segmentation mean error under the noisy patterns.
In the future, we will use the Gaussian norm instead of the L 2 norm to constrain the distance between gradients and orientation vectors, which may generate a more accurate orientation map. Furthermore, we will optimize the algorithm and code in C++ to reduce the execution time for real-time applications.