Fringe Order Correction for Number-Theoretical Phase Unwrapping Method via Maximum Likelihood Principle

The number-theoretical phase unwrapping method has recently been widely applied in fringe projection profilometry. But fringe order errors may occur due to noise or distortion, leading to errors in the unwrapped phase map, and eventually affecting the accuracy of the reconstructed object surface. In this paper, we propose a novel fringe order correction method based on the maximum likelihood principle. The direct cause of fringe order error is the deviation of an intermediate variable which in theory should be an integer, and the ground truth of the integer stays unchanged within a valid neighborhood. By modeling the calculated intermediate variable as an observed sample from the normal distribution of the unknown ground truth integer, we can determine a valid neighborhood relative to the observed pixel. Then the ground truth integer can be calculated by maximizing the likelihood function and then the fringe order error is corrected. The simulation results and experimental comparisons have verified the feasibility, robustness, and superiority of the proposed method in contrast with other fringe order correction methods.


I. INTRODUCTION
Fringe projection profilometry (FPP) is an important technique for three-dimensional measurement due to its accuracy and high efficiency [1]. The basic FPP system is composed of a camera and a projector. The principle of the FPP technique is to project fringe patterns onto the object surface, calculate phase information of the images captured by the camera in another direction, and then the height of the object surface is obtained by phase-height calibration [2]- [4]. The phase calculated by either phase-shifting algorithm [5], [6] or a transform-based method [7], [8], however, is wrapped in the range of -π to π, and thus the process of phase unwrapping must be carried out to obtain the continuous phase map by adding multiple integer numbers of 2π. The number of 2π for each pixel is called fringe order, and the main challenge of phase unwrapping is calculating the fringe order of each pixel accurately and efficiently.
The associate editor coordinating the review of this manuscript and approving it for publication was Zhenbao Liu . Dozens of phase unwrapping methods have been proposed in recent years, which can be generally divided into two main categories: spatial methods [9]- [11]and temporal methods [12], [13]. The spatial methods unwrap the phase of a pixel by referring to phase information of the spatial neighboring pixels, they tend to fail when the measured objects contain discontinuities or separations. The temporal methods are more robust owing to their pixel-bypixel phase unwrapping procedure and are widely used in high-precision measurement scenarios. The temporal phase unwrapping method can be categorized into gray-code [13], multifrequency [14], multiwavelength [15], and numbertheoretical [16]- [22] approaches. Among these approaches, the number-theoretical method shows better reliability [12] with only bi-frequency fringe patterns needed.
However, fringe order errors (FOE) may occur due to noise or distortion [23], leading to errors in the unwrapped phase. Most previous studies concentrated on reducing noise/distortion [24]- [26] or increasing the tolerance of noise [27]- [29] while other studies focused on eliminating VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the FOEs directly [30]- [33]. Ding et al. [30] applied the majority rule method to determine the most likely order within the interval, but this method is unstable when intervals are small. Assuming that the fringe order changes smoothly, the least square method was applied by referring to the fringe orders of adjacent pixels [31], but this assumption is untenable when the measured surface is complex. Lu et al. [32] used a fringe order gradient map to find false positive pixels and exclude small regions, but it is hard to distinguish the pixel where the FOE occurs from the pixel where the fringe order is correct. Kam et al. [33] employed the k-nearestneighbor search method and spatial comparison method, but the false-positive correspondences may increase because of ambiguity in the phase. Unlike the above fringe order correction methods, this paper proposes a method that focused on the intermediate variable ψ that determines the fringe orders. The ground truth of the observed variable stays unchanged within the neighborhood in most cases even though the surface is complex, inspiring an idea to combine the spatial method with the temporal method: the fringe order for each pixel can be corrected by maximizing the likelihood function related to the neighborhood, yielding relatively better robustness and feasibility.

II. PRINCIPES OF NUMBER-THEORETICAL PHASE UNWRAPPING A. PHASE-SHIFTING ALGORITHM
The digital projector is usually employed as the light source to project the sinusoidal fringe patterns in FPP. The analytical expression of the sinusoidal fringe pattern in the N-step phase-shifting algorithm can be expressed by: where x P , y P is the coordinate of the projected pixels, a P is the direct current component of the intensity, b P is the amplitude, λ P 0 is the spatial wavelength of the sinusoidal signal, n is the phase-shift index, and N is the total number of phase-shift steps.
After projecting these patterns onto the object surface, the phase signals are modulated, and the images captured by the camera can be expressed as: where A (x, y) is the average intensity of the image, B (x, y) is the intensity modulation, and (x, y) is the phase containing the height information of the measured object. From the captured images, the wrapped phase φ (x, y) can be retrieved as follows: Since the phase is calculated pixel-by-pixel, the pixel coordinate index (x, y) is removed from equations henceforth in this paper.

B. BI-FREQUENCY NUMBER-THEORETICAL PHASE UNWRAPPING METHOD
The output of the four-quadrant inverse tangent function in (3) ranges from -π to π, and thus the phase solved by (3) is wrapped. To unwrap the phase using temporal phase unwrapping methods, multiple phases derived from fringe images in different frequencies are needed. For the bi-frequency method, the wrapped phase derived from high-frequency and low-frequency phase-shifting images are represented as φ h and φ l respectively, and the corresponding unwrapped phases are h and l , which can be calculated by: where k h and k l represent the fringe orders. The essential task of the phase unwrapping algorithm is to determine the fringe orders accurately and efficiently. According to the numbertheoretical phase unwrapping method [16], there exists a unique mapping from fringe order pair (k h , k l ) to wrapped phase information if the following inequality holds: where W is the total length along the sinusoidal signal (i.e., the horizontal resolution of the projector) and LCM() is the function that outputs the least common multiple of input variables. The unwrapped phase h and l have the relationship as follows: Define p h = LCM (λ h , λ l ) λ h , p l = LCM (λ h , λ l ) λ l , and then (6) can be rewritten as: Combining (4) and (7), we have: where k h p l − k l p h is an integer, and thus ψ in theory should also be the same integer. The results of the k h p l − k l p h can be pre-calculated and restored as entries in a lookup table (LUT) [17]. The closest integer of ψ, is uniquely mapped to an orders pair; examining the closest integer of ψ in the table, the fringe order pair is determined by (9), and the unwrapped phase can be calculated by (4) eventually.

C. FRINGE ORDER ERROR
In the number-theoretical phase unwrapping method, ψ is rounded to the closest integer. However, if the unwrapped phase φ h or φ l is corrupted, the error of ψ may exceed 0.5, and the corruption will lead to the wrong integer and then the wrong pair of fringe order. These phase errors may result from sensor noise, gamma distortion, or projector/camera defocus.
In the works of predecessors [12], [34]- [36], the Gaussian model is widely adopted to describe the error of phase. Considering all the effect of phase error as a perturbation on the true phase, the variance of the wrapped phase φ h and φ l can be calculated as: where σ 2 I is the variance of the intensity of the fringe image captured by the camera. According to (8), the variance of ψ can be calculated as: Combined with statistical analysis, the n-sigma criterion of reliability theory is applied to estimate the maximum error of ψ: ψ max = nσ ψ (99.73% confidence level when n = 3). The FOE occurs when | ψ max | > 0.5, thus, we can find the upper boundary of wrapped phase variance by (13):

III. FRINGE ORDER CORRECTION VIA MAXIMUM LIKELIHOOD PRINCIPLE A. THE VALID NEIGHBORHOOD
Under the circumstance that the wrapped phase is affected by noise and distortion, for each pixel, we characterize ψ as a random variable subjected to the normal distribution: whereψ is the ground truth integer, and ψ * is the observed value of ψ. Without loss of generality, we suppose a fringe pattern with 600 pixels in horizontal resolution, and the wavelengths of bifrequency fringes are λ h = 16 pixels and λ l = 39 pixels. The ideal distribution ofψ along the horizontal direction can be calculated by (8), and is shown in Fig. 1. We can see thatψ is distributed as several stairs along with the horizontal axis, and pixels of each stair share the sameψ, which is different from any other stairs. In other words, if the neighboring pixels share the same fringe order pair k h and k l , according to (8), they should share the sameψ. For a certain observed value ψ * x 0 ,y 0 of the pixel (x 0 , y 0 ), the observed group ψ of neighboring pixels is denoted as: where m is the total number of pixels of the neighborhood, ε x and ε y are the coordinate limits of the neighborhood. When all neighborhood pixels belong to the same fringe order k h in high-frequency images and k l in low-frequency images, the group ψ can be considered as repeated observation of ψ * x 0 ,y 0 . In this paper, we call ψ the valid neighborhood of ψ * x 0 ,y 0 . However, when the neighborhood covers pixels from other fringe orders, there will be at least one jumping change in the group ψ, and in this case the ψ is the invalid neighborhood of ψ * x 0 ,y 0 . As for the example given in Fig. 1, the difference of adjacent pixels is shown in Fig. 2. For the 599 pixels in horizontal direction, the difference is zero in most cases (547 pixels), while the difference of some pixels has jumping change (p h in 14 pixels, −p l in 37 pixels and p h − p l in 1 pixel). The neighborhood is invalid if it covers these jumping pixels. According to the property of chi-square distribution, if the neighborhood is valid, the following probability formula is satisfied: where s 2 is the sample variance of the neighborhood, α is the significance, in this paper we use α = 0.001 to achieve 99.9% confidence. We use this condition to determine valid neighborhood: the neighborhood is valid only if the following inequality holds: Otherwise, the neighborhood is invalid. When the neighborhood is invalid, there exist two or more ground truth integers ofψ in the neighborhood, or outliers VOLUME 8, 2020 caused by noise. K-means clustering method can divide ψ into clusters according to the condition of minimum sample variance. As with the prior example, the difference between ground truth integers is n h p h − n l p l , (n h , n l ∈ N). Then, by changing the values of one or more clusters, we can turn the invalid neighborhood into a valid neighborhood. The detail of the procedure is as follows: Step 1: Make cluster analysis of ψ in minimum-variancecriterion, and ψ is divided into ψ 1 , ψ 2 , ψ 3 . . . .
Step 2: Find the clusters that contain only one element or do not satisfy (16), and then exclude them from the ψ.
Step 3: Treat the cluster which includes the currently observed value ψ * x 0 ,y 0 as the target cluster. If the target cluster is excluded in Step 2, replace it with the cluster whose mean value is closest to ψ * x 0 ,y 0 .
Step 4: Change the value of elements of other clusters by adding n h p h − n l p l , (n h , n l ∈ N). If the union cluster of the post-change cluster and target cluster satisfies (16), the changed cluster is retained. Otherwise, the changed cluster is discarded. Therefore, the reunion of remaining clusters and target clusters is a valid neighborhood of ψ * x 0 ,y 0 . To elaborate on the procedure, we add Gaussian noise (σ 2 ψ = 0.143) to the previous distribution of ψ in Fig.1. We choose the pixel 296, which is on the edge of jumping change, as the observed pixel. The 1 × 5 sized neighborhood of ψ * 296 is: ψ = (ψ * 294 , ψ * 295 , ψ * 296 , ψ * 297 , ψ * 298 ) = According to (13), the values in the valid neighborhood can be considered as samples from the normal distribution with unknown meanψ. The probability density function relative to observed ψ * is defined as The likelihood function is defined as the product of all individual probability of the values in the neighborhood: Since the unknown meanψ should be an integer, the likelihood function is discrete, and the problem can be simplified as:ψ whereψ is the most likely integer that maximizes the likelihood function; ψ 1 ,ψ 2 ,ψ 3 . . . are several integers that are close to round (ψ * ), and can be considered as possible candidates of ground truth value. Ifψ is not equal to round (ψ * ), the fringe orders previously calculated by (9) are considered incorrect, and then fringe orders are recalculated by replacing round (ψ * ) in (9) withψ. As for the example given in Section III(A), the procedure of order correction is detailed as follows: Fringe order pair of pixel 296 is previously determined by LUT (7) because round ψ * 296 = 7. The candidate ground truth integers of ψ * 296 are (6,7,8). The likelihoods for candidate integers are calculated by (18). Since L (8) is larger than L (6) and L (7), the fringe order pair of pixel 296 is corrected by LUT (8).

IV. SIMULATIONS A. SELECTION OF NEIGHBORHOOD SIZE
In the proposed method, the fringe orders are calculated pixel-by-pixel. For each pixel, the values of ψ * in the valid neighborhood are used to estimate the most likely fringe order. Thus, the selection of neighborhood size is an essential algorithmic component of the proposed method, which is determined by ε x and ε y in (14). A small-size neighborhood may be susceptible to noise or distortion, while a large-size neighborhood will increase computational complexity and is likely to cover pixels from other fringe orders. The neighborhood size should be selected as the trade-off of accuracy and efficiency.
The simulation platform is programmed by MATLAB and the computation processor is Inter Core i7-7700 CPU. The simulated surface of peaks distributed at 600 × 400 pixels is generated by the computer and considered as the measured object. We add noise with standard deviation 12 to the simulated images, and we set λ h = 16 pixels, and λ l = 39 pixels. Since fringe order calculated by noise-free fringe images is absolutely correct, success rate (SR) of fringe order calculation under noise conditions can be obtained by comparison. The fringe orders are calculated using different neighborhood sizes for five times, and the average of calculated SR and time consumption results are shown in Table 1. As can be seen, SR by and large increases as we set neighborhood size  larger, but the SR decreases and time consumption increases (e.g., the 1 × 9 and 7 × 7 neighborhood) when neighborhoods with too many columns may cover pixels from different fringe orders, resulting in trial-and-error determining the valid neighborhood. Through this paper, we choose the 3 × 3 sized neighborhood as the trade-off between SR and the time consumption.

B. ANTI-NOISE PERFORMANCE
According to (11), the smaller p 2 h + p 2 l , the less likely fringe orders are affected by noise. In other words, FOEs are determined not only by the noise of fringe images but also by fringe-frequency-selection strategies.
The simulated measured object is the same as Section IV(A). To verify the anti-noise performance of the proposed approach, fringe order results under two different levels of noise, applying two different fringe-frequencyselection strategies are simulated in this section: level 1: the standard deviation of image noise is 8; level 2: the standard deviation of image noise is 12. Strategy 1: λ h = 16 pixels, λ l = 39 pixels, p h = 39, p l = 16; strategy 2: λ h = 28 pixels, λ l = 88 pixels, p h = 22, p l = 7. For comparison, the fringe order result is corrected by the majority rule method [30], the least square method [31], and the proposed method respectively. The corrected fringe order results are shown in Fig. 3. The simulation environment is different for each row of Fig. 3; the uncorrected results, and the results corrected by the majority rule method, the least square method and the proposed method are shown separately in columns of Fig. 3. Through the comparison of Fig. 3(a)-(b) and Fig. 3(c)-(d), the fringe orders are easily corrupted by noise, but the effect of noise can be reduced using selected fringe frequency. Through the comparison of SR and phase distribution in different columns, the proposed method outperforms the other two methods. As illustrated in Fig. 3(b1)-(b4), the proposed fringe order correction method maintains good performance against other methods even under the high noise level.

C. APPLICABILITY TO COMPLEX SURFACES
Many fringe order correction methods [31] are based on the assumption that fringe orders or unwrapped phases VOLUME 8, 2020 change smoothly. However, complex surfaces may have huge jumps in height, causing discontinuities in fringe orders or unwrapped phases.
A simulated complex surface with multiple steps is considered as the measured object in this section, as shown in Fig. 4(a). The applied noise level 2 and fringe-frequencyselection strategy 1 is detailed in Section IV(B), and other simulation environments are the same as Section IV(B). The fringe order result is corrected by the least square method [31] and the proposed method respectively. By subtracting the original noise-free fringe order result from the corrected fringe order results, the residual FOEs are obtained, as shown in Fig. 4(b) and (c). The SR is 95.61% for the least square method and 99.36% for the proposed method. The result proves that the proposed method is not appreciably affected by surface step edges, and that it is more applicable to complex surfaces.

V. ACTUAL EXPERIMENTS
A. EXPERIMENTAL ENVIRONMENT An experimental FPP system is set up to verify the proposed method, as shown in Fig. 5. The system mainly consists of a digital projector (TENGJU TJ21HB, resolution: 1280 × 720), a CMOS camera (HIKVISION MV-CA013-21UM, resolution: 1280 × 1024), and a computer for image processing. A lens (COMPUTAR M0814-MP2) with a focal length of 12 mm was attached to the camera. The measured object is placed on a reference plane. The algorithms are programmed by MATLAB and the computation processor is Inter Core i7-7700 CPU.
In the experimental implementation, the wrapped phases are calculated by N-step phase-shifting algorithm, to achieve relatively better precision and robustness, and the number of steps is set to 4 unless we specified otherwise. The shadow or low reflective areas of the object usually yield unreliable and useless wrapped phase results, and these areas can be removed by the preset threshold of the intensity modulation and in Eq. (2) [6]. For each pixel, if their intensity modulation is under 20, the phase result will be eliminated.
The unwrapped phase increases along the horizontal axis, so it is not easy to distinguish the phase which is modulated by the measured object from the background. For better illustration in our experiments, the phase map of the background reference plane is calculated in advance, and then the background phase is subtracted from all phase maps of the measured object.

B. ESTIMATION OF WRAPPED PHASE VARIANCE
As we discussed in Section III(A), the variance of ψ is necessary to verify the valid neighborhood. According to (12), the wrapped phase variance σ 2 φ should be pre-calculated. In synthetic experiments presented in Section IV, we can simulate the image noise and calculate σ 2 φ by (10). However, error sources that come from the real-world are unpredictable, and in this section, we present a simple method to estimate the real-world σ 2 φ . For the phase-shifting algorithm introduced in Section II(A), if the number of steps is 6 at least, we can compute two uncorrelated wrapped phase results. The difference between the two results is consistent with a normal distribution: 157722 VOLUME 8, 2020 where φ 1 and φ 2 are the two wrapped phase and phase shift is 2π n. For example, the 8-step phase-shifting images can generate two wrapped phase results, φ 1−3−5−7 (calculated by 4-step phase-shifting algorithm using 1st, 3rd, 5th and 7th images) and φ 2−4−6−8 (calculated by 4-step phase-shifting algorithm using 2nd, 4th, 6th and 8th images). σ 2 φ can be estimated by the variance of the difference between φ 1−3−5−7 and φ 2−4−6−8 − π 4.
In this experiment, three objects with different materials (paper, wood, and metal) are tested in different measuring conditions. The horizontal resolution of the projector is W = 1280 pixels, fringes with different wavelengths (W , W 2, W 4, W 8, W 16, and W 32) are projected onto the tested objects, and then the wrapped phase variances σ 2 φ are estimated, as shown in Fig. 6. It can be easily concluded that σ 2 φ is not significantly affected by the wavelength λ, though vulnerable to the reflection rate of tested material, hardware device parameters (e.g., camera exposure) and measurement strategy (e.g., number of steps of phase-shifting).
In actual measuring situations, poor measurement conditions and unsuitable surfaces may lead to huge σ 2 φ , resulting in massive FOEs. Thus, hardware device parameters and measurement strategy should be adapted to the measured surfaces. However, it is the precision of σ 2 φ estimation instead of σ 2 φ itself that affects the feasibility of the proposed method. A further discussion is covered in the next section.

C. ANALYSIS OF ALGORITHMIC COMPONENT
The efficiency and accuracy of the proposed method are determined by several algorithmic components. In this section, experiments are performed to analyze the influences of several controllable algorithmic components. We choose a metal part as the measured object because the reflective nature of its surface can stress the proposed method in case of high noise levels.
In the first experiment, the influence of the precision of σ 2 φ estimation is analyzed. The fringe-frequency-selection strategy is set as λ h = 30 pixels, λ l = 129 pixels, and the pre-estimated value of σ 2 φ in this experiment is 0.00073. The captured images and phase map results calculated using different σ 2 φ are shown in Fig. 7. The FOEs are counted in the region of interest (marked as a rectangle in Fig. 7(c)) of each map. The counted FOEs for Fig. 7 (h) 2. The results show that the corrected phase map has the least FOEs when the wrapped phase variance is correctly estimated. When wrapped phase variance is under-estimated, the FOE rate rises, while the FOE rate does not rise noticeably when wrapped phase variance is over-estimated. The reason is as follows: the step differences of ψ are far more than 1 in most cases, and thus an   over-estimated wrapped phase variance can still distinguish the invalid neighborhoods. In practice, we prefer an overestimated wrapped phase variance to reduce the influence of wrapped phase variance.
In the second experiment, the influence of the fringefrequency-selection strategy is analyzed. The wrapped phase variance in this experiment is set as 0.0015. The phase map results calculated using different fringe-frequency-selection strategies are shown in Fig. 8. The FOEs are counted in the region of interest (marked as a rectangle in Fig. 8(a)) of each map. The counted FOEs for Fig. 8

D. EXPERIMENTAL COMPARISONS
In this section, we present experimental comparisons of different fringe order correction methods. The fringe-frequency-selection strategy is set as λ h = 30 pixels and λ l = 129 pixels. Several representative objects are placed on the reference plane, including a frosted plastic mouse, a wooden comb, and a dark-color plastic calculator. The fringe order results are corrected by the majority rule method, the least square method, and the proposed method separately. The unwrapped phase map obtained using three-frequency phase unwrapping method can be considered as ground truth. For each pixel, if the absolute difference between the calculated phase and ground truth exceed π, the pixel is marked as FOE pixel.
The phase map results are illustrated in Fig. 9. For a better visual presentation, we set boundaries for the phase of FOE pixels, and thus the FOEs are denoted black dots or white dots in the figures. Fig. 9(b-2) and (c-2) show the incapacity of majority rule method to deal with order steps with small intervals. The least square method shows poor robustness when dealing with pixels close to surface step edges, such as comb teeth shown in Fig. 9(b-3). Fig. 9(a-4), (b-4), and (c-4) show that there exist fewer black dots in phase maps calculated by the proposed method than other methods. The three-dimensional phase distribution results in Fig. 9(a-5), (b-5), and (c-5) are in shape with the measured objects and can be further mapped to height distribution after phase-height calibration.
For a clearer comparison, the FOEs are counted for each method. The histogram of FOEs for the above experiments is shown in Fig. 10.
In the case of mouse, the surface is smooth and quasi-Lambertian. Compared with the uncorrected phase map, all the three methods show effectiveness at the smooth area. However, in the case of comb, the surface is complex with narrow teeth. The FOEs of pixels close to teeth edges are easily spread to neighboring pixels. Since the effect of edge pixels is reduced via maximum likelihood estimation, the proposed method is more applicable in complex surfaces. In the case of calculator, more original FOEs are aroused by noise because of the low-reflection of surface. The FOEs are significantly reduced by the proposed method compared with other methods. In conclusion, the comparison results are able to prove that the proposed method provides better feasibility and robustness for fringe order correction.

VI. CONCLUSION
In this paper, we present a novel fringe order correction method for the bi-frequency number-theoretical phase unwrapping. Specifically, we consider the calculated intermediate variable which determines the fringe orders as an observed sample from the normal distribution with an unknown ground truth integer. Via the pre-estimated wrapped phase variance and the condition of chi-square distribution, the valid neighborhood of the observed pixel is determined. By introducing the maximum likelihood principle to the neighborhood, the most likely correct integer is calculated, and the FOEs are corrected eventually. We performed simulations and actual experiments to stress the robustness and feasibility of the proposed method, with the fringe order errors significantly reduced. The comparison results demonstrate that the proposed method is capable to unwrap phase accurately, even in high-noise and complex-shape environments.