## Abstract

In fringe projection profilometry, the purpose of using two- or multi-frequency fringe patterns is to unwrap the measured phase maps temporally. Using the same patterns, this paper presents a least squares algorithm for, simultaneously with phase-unwrapping, eliminating the influences of fringe harmonics induced by various adverse factors. It is demonstrated that, for most of the points over the measured surface, projecting two sequences of phase-shifting fringe patterns having different frequencies enables providing sufficiently many equations for determining the coefficient of a high order fringe harmonic. As a result, solving these equations in the least squares sense results in a phase map having higher accuracy than that depending only on the fringe patterns of a single frequency. For the other few points which have special phases related to the two frequencies, this system of equations becomes under-determined. For coping with this case, this paper suggests an interpolation-based solution which has a low sensitivity to the variations of reflectivity and slope of the measured surface. Simulation and experimental results verify that the proposed method significantly suppresses the ripple-like artifacts in phase maps induced by fringe harmonics without capturing extra many fringe patterns or correcting the non-sinusoidal profiles of fringes. In addition, this method involves a quasi-pointwise operation, enabling correcting position-dependent phase errors and being helpful for protecting the edges and details of the measurement results from being blurred.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Phase-shifting technique, a high-resolution fringe pattern analysis tool, is popularly used in the fields like interferometry [1–5], holography [6], moiré topography [7–9], fringe projection profilometry [10–16], and fringe deflectometry [17,18]. This technique requires capturing a sequence of fringe patterns having different phase shifts and recovers phases from them by using a pointwise algorithm. In the past decades, there have been a number of phase-shifting algorithms developed for restraining the effects of various adverse factors [19–21]. Usually, these factors induce phase errors by distorting sinusoidal fringe profiles, i.e., by introducing harmonics into the fringe signals.

With phase-shifting fringe projection profilometry [10–16], many factors cause harmonics in fringe patterns. Among them, the camera and projector nonlinearities [22] are recognized as one of the most crucial, especially when a low-cost profilometer is used. These fringe harmonics induce ripple-like artifacts on the calculated phase maps. To solve this problem, the most straightforward method is to perform photometric calibrations to both projector and camera, and actively compensate for their nonlinear response using the calibration results in advance [23,24]. Some methods project sinusoidal fringes onto a white plane and record the patterns for the same purpose [25–27].

Despite being calibrated, the camera and projector in a measurement system can never be perfectly linear. Even if using some commercial profilometers having highly linear devices, we still observe, from their raw data of measurement, tiny high-frequency artifacts parallel to the fringe directions. These artifacts are typically induced by fringe harmonics. In this case, developing a software-based compensating method is helpful for further improving measurement accuracy. For example, once the response of the camera is calibrated to be linear (we calibrate the camera, because it is much easier to calibrate than a projector), the phase errors purely induced by the projector nonlinearity are phase-dependent and can be compensated for by using a look-up table (LUT) [28,29] or an error function [30–33] in the data processing stage. These software-based methods help measurement break through the accuracy limitations of hardware. Note that, when the camera and projector have nonlinearities simultaneously, the overall luminance response of the measurement system, depending not only on the device nonlinearities but also on the reflectivity and slope variations of the measured object surface, is a function of pixels. For this reason, a pointwise compensating algorithm is more preferable for coping with the nonlinearity issue.

Besides the device nonlinearities, fringe harmonics may come from manifold factors. For example, some dynamic measurement techniques project specially-designed binary patterns [34–37] at high rates by taking advantage of digital micromirror device (DMD) projectors. With them, properly defocusing the binary patterns generates approximately sinusoidal fringes in the measurement space, but these defocused fringes may have harmonics depending on the defocus level related to the object depths [38]. Another typical example is that, in color-coded fringe projection techniques, the cross-talk between chromatic channels may induce harmonics in each separated fringe patterns [39]. In addition, fringe distortions caused by camera saturation [13] may occur when measuring a high dynamic range object [40]. These harmonic-induced errors are position-dependent, thus requiring a pointwise compensating algorithm for suppressing their effects.

Phase-shifting technique has potential for solving the aforementioned problems because it involves only pointwise operation. In fact, most phase-shifting algorithms correspond to finite impulse response (FIR) filters for extracting the fundamental frequency component from a temporal signal. Spectral analysis is helpful for developing phase-shifting algorithms robust to harmonics [41]. However, performance of a phase-shifting algorithm in restraining the effects of harmonics strongly depends on the number of phase shifts. For example, when using the synchronous detection algorithm [1] which has equal-spaced phase step over a 2π period, the *N*-step algorithm can remove the influence of the harmonics up to the (*N*-2)th order [6,42]. According to this principle, we can suppress the effects of the fringe harmonics by simply increasing the number of phase shifts, at the expense of time duration for capturing more fringe patterns.

Requiring likewise capturing more fringe patterns, fringe projection profilometry usually uses two- or multi-frequency fringe patterns for determining the absolute fringe orders, thus unwrapping the measured phases [43–45]. Using these fringe patterns, however, the measurement accuracy is determined only by the fringe patterns having the highest frequency; whereas the low-frequency fringes serve just as an aid for providing the reference phases having a wide unambiguous range to help unwrap the phases of high frequency fringes. Note that, with these temporal phase-unwrapping techniques, the fringe patterns having different frequencies are possibly subject to the same adverse factors, for example, to the same nonlinearities of projector and camera or to the same defocus level when projecting binary patterns. In this situation, the harmonics in the fringe patterns of different frequencies may have the same coefficients, implying a possibility of deriving a pointwise algorithm for eliminating effects of such fringe harmonics.

For confirming the presumption above, this paper presents, to the best of our knowledge, a new algorithm of fringe harmonics elimination suitable for multi-frequency phase-shifting fringe projection profilometry. It uses two fringe pattern sequences having different frequencies for, simultaneously with phase-unwrapping, eliminating the influences of fringe harmonics induced by some adverse factors. In implementation, it enables accurately calculating the phase value at an ordinary point on the measured surface by solving a system of equations in the least squares sense. This system of equations is determined by the fringes having different frequencies and includes the coefficient of a high order harmonics as an unknown. For the special points which have certain phase values related to the frequencies, this system of equations becomes under-determined, and an interpolation-based solution is suggested for coping with this case. Simulation and experimental results demonstrate this proposed method to be effective in restraining effects of fringe harmonics without capturing extra many fringe patterns or correcting the non-sinusoidal profiles of fringes. This technique can be used not only in fringe projection profilometry, but also in fringe deflectometry [17,18] as well as multi-wavelength interferometry [3–5] for improving measurement accuracies.

## 2. Measurement system and effects of fringe harmonics

The measurement system, as showed in Fig. 1(a), is mainly composed of a camera, a projector, and a computer. Using *N*-step synchronous detection algorithm [1] which requires equal-spaced phase steps over a 2π period, the *n*th (*n* = 0, 1, …, *N*−1) fringe pattern generated in the computer is described as

*u, v*) denote the pixel coordinates of a point on the image plane of the projector, then

*g*(

_{n}*u, v*) is the gray level at (

*u, v*), and

*p*is the fringe pitch. The positive numbers ${a_0}$ and ${a_1}$, satisfying $0 \le {a_0} \pm {a_1} \le 1$, denote the bias and the contrast of the fringes, respectively.

As mentioned in introduction, using the measurement system in Fig. 1(a), many factors may cause harmonics in fringe patterns. To gain insight into the principle of the proposed technique in restraining effects of fringe harmonics, we presume for the moment that these harmonics are induced by the camera and projector nonlinearities. If the harmonics are caused by other factors, the proposed method is still effective as long as these factors affect the multiple frequency fringe patterns in the same fashion. Figure 1(b) shows the transfer of the fringe patterns in the measurement system. The projector may have nonlinearity, so we reasonably assume that its output luminance is a continuous function of input gray level. In this case, this function, according to the Weierstrass approximation theorem, can be uniformly approximated as closely as desired by a polynomial. Therefore, the patterns casted by the projector driven by the fringes in Eq. (1) is expressed as

with ${c_k}$ being its coefficients. When the fringe patterns in Eq. (2) are casted onto the object, the reflected brightness depends on the local colors, textures, and reflectivity, as well as on slopes of the surface according to the Lambertian model. Simultaneously, it is influenced by ambient illumination. Considering all these factors, the reflected brightness can be simply expressed as where (*x, y*) denotes the point coordinates illuminated by (

*u, v*). $R({x,{\; }y} )$ and $B({x,{\; }y} )$, both are functions of (

*x, y*), denote reflection coefficient and ambient brightness, respectively. The camera is also considered to be a nonlinear device, so its captured intensity is expressed using another polynomial of the form with ${d_k}$ being its coefficients. From Eqs. (1) through (4), it is easy to know the captured patterns are not sinusoidal but contains harmonics, namely,

Figure 2 exemplarily shows the transfer of the fringe patterns in the measurement system. In it, Fig. 2(a) simulates a standard one-dimensional (1D) sinusoidal fringe pattern by use of Eq. (1) with ${a_0} = 0.5$, ${a_1} = 0.4$, and $p = 512/9$ pixels. Assuming the projector has a nonlinearity represented by the curve in Fig. 2(b), its projected fringes become nonsinusoidal as plotted in Fig. 2(c). For simulating the effects of the varied reflectivity and slope over the measured surface, we set the ununiformly distributed reflection coefficients $R({x,{\; }y} )$ to be Gaussian shaped having a 50% decrease at the ends on both sides of the pattern. The additional intensity caused by the ambient illumination is $B({x,{\; }y} )= 0.1$. As a result, the reflected brightness of the surface is shown in Fig. 2(d). Using a camera having a nonlinear response curve in Fig. 2(e), the captured fringe pattern is given in Fig. 2(f). This figure shows the characteristics of the captured fringe patterns when both the camera and projector have nonlinearities. First, the captured fringe patterns are no longer sinusoidal, but have high order harmonics. Second, the amplitudes of the fringe harmonics are not uniform, but functions of the pixel coordinates.

## 3. Phase measuring

#### 3.1 Phase-shifting algorithm

Using synchronous detection algorithm [1], the phases in Eq. (5) are estimated as

*N*-2)th order. Increasing the number of phase shifts is helpful for restraining the effects of fringe harmonics. On the other hand, we prefer to use as few as possible fringe patterns for improving the measurement efficiency. Simultaneously with calculating the phases, we can calculate the background intensities and the fringe modulations as and

#### 3.2 Phase-unwrapping using two-frequency fringe patterns

The phases calculated from Eq. (6) lie within the principal range from -π to π radians because of the arctangent function involved. In fringe projection profilometry, the purpose of phase-unwrapping is not only to recover continuous phase maps, but also to determine absolute fringe orders for each pixel since the relations between the fringe phases and the object depths are not linear. The temporal phase-unwrapping techniques [43], compared with spatial techniques [45], are more suitable for this purpose because of their pointwise operation.

The temporal phase-unwrapping techniques require capturing at least two fringe pattern sequences having different frequencies. Using Eq. (1), we generate two fringe pattern sequences have fringe pitches ${p_\textrm{I}}$ and ${p_{\textrm{II}}}$ with ${p_\textrm{I}}\;<\;{p_{\textrm{II}}}$. For the convenience purpose, we define their ratio as

## 4. Harmonics elimination technique

#### 4.1 Least-squares algorithm

Using the synchronous detection algorithm [1], the *N*-step algorithm can remove the influence of the harmonics up to the (*N*-2)th order [6,42]. Consequently, the measurement result will be affected mainly by the (*N*-1)th order harmonic. It is possible to cancel the influence of the (*N*-1)th order harmonic by solving its amplitude from the fringe patterns. For this purpose, we include the (*N*-1)th order harmonic in the equations, resulting in

*N*. The coordinate notation $({x,{\; }y} )$ is omitted in Eq. (14) for shortening the expressions.

Based on Eq. (14), we have 2*N* equations involving four knowns, i.e., ${b_0}$, ${b_1}$, ${b_{N - 1}}$, and ${\varPhi _\textrm{I}}$, for each pixel $({x,\; y} )$. Note that the terms related with ${b_2}$ through ${b_{N - 2}}$ are not included in Eq. (14), because the *N*-step algorithm is insensitive to the harmonics up to the (*N*-2)th order. These equations are blind to these terms. In other words, their absence does not affect the solutions of the four unknowns just mentioned. The problem is that the system of equations based on Eq. (14) is nonlinear and hence not easy to solve. For overcoming this difficulty, we linearize Eq. (14) by using Taylor series expansion and solve it in an iterative way. Denoting the estimates of ${b_0}$, ${b_1}$, ${b_{N - 1}}$, and ${\varPhi _\textrm{I}}$ as ${\hat{b}_0}$, ${\hat{b}_1}$, ${\hat{b}_{N - 1}}$, and ${\hat{\varPhi }_\textrm{I}}$, respectively, we define their errors as

*T*denotes transpose. Using the estimates of unknowns, we can repeat calculating the fringe intensities using

*N*equations involving four unknowns. Its least squares solution can be calculated by solving the normal system of equations as follows. which contains exactly four equations having four unknowns. Based on Eq. (22), we calculate ${b_0}$, ${b_1}$, ${b_{N - 1}}$, and ${\varPhi _\textrm{I}}$ through the following iterative procedure.

**Step 1**. Determine the initial values.

Here, we use the values calculated in Section 3.1 as the initial values for ${b_0}$ and ${b_1}$, i.e., $ b_0^{(0 )} = {\hat{b}_0}$, $b_1^{(0 )} = {\hat{b}_1}$, use the unwrapped phase map in Section 3.2 as the initial values for ${\varPhi _\textrm{I}}$, i.e., $\varPhi _\textrm{I}^{(0 )} = {\hat{\varPhi }_\textrm{I}}$, and assign $b_{N - 1}^{(0 )} = 0$, where the superscript in the pair of brackets denotes the number of iterations.

**Step 2**. If the *i*th iterative results [$b_0^{(i )}$, $b_1^{(i )},{\; }b_{N - 1}^{(i )}$, $\varPhi _\textrm{I}^{(i )}$] are calculated, use them to estimate the (*i *+ 1)th iterative results.

Substituting [$b_0^{(i )}$, $b_1^{(i )},{\; }b_{N - 1}^{(i )}$, $\varPhi _\textrm{I}^{(i )}$] instead of $[{\hat{b}_0},{\hat{b}_1},{\hat{b}_{N - 1}},{\hat{\varPhi }_\textrm{I}}$] into Eq. (16), we calculate $\hat{I}_\textrm{I}^{(i )}$ and $\hat{I}_{\textrm{II}}^{(i )}$ for $n = 0,{\; }1,{\; }\ldots ,{\; }N - 1$, and further obtain the vector ${{\mathbf \Delta }^{(i )}}$ through Eq. (17). Substituting [$b_0^{(i )}$, $b_1^{(i )},{\; }b_{N - 1}^{(i )}$, $\varPhi _\textrm{I}^{(i )}$] into Eq. (19), we calculate the partial derivatives, and then obtain the coefficient matrix ${{\mathbf \Omega }^{(i )}}$ through Eq. (21). Solve the system of equations

**Step 3**. Repeat implementing Step 2 until the algorithm converges.

Continue the simulation in Fig. 3. Using the least squares algorithm in this subsection, we eliminate phase errors induced by the second order harmonics of the fringes. The result after 20 iterations is shown in Fig. 3(e). Its residual errors obtained by subtracting Fig. 3(a) is shown in Fig. 3(f). From these results, we observe that the large ripple-like errors induced by the fringe harmonics have been significantly suppressed for most points over the whole phase curve. At three special points whose phases plus carrier are -6$\pi $, 0, and 6$\pi $ radians, however, the algorithm fails in eliminating the errors. This fact means that the system of equations we used in this subsection is not always well-determined, and in other words its equations are possibly not independent of one another at some points having special phases. These special points give pitfalls to the least-squares algorithm proposed in this subsection. In the next subsection, we shall cope with this issue.

#### 4.2 Pitfall point recognition

As just mentioned, the above least-squares algorithm may fail in eliminating the errors at some special points. We investigate this phenomenon by checking the coefficient matrix of Eq. (22), which can be explicitly written as

It is known from Eq. (25) that in case

and simultaneously we have ${\omega _1} ={\pm} N$ and ${\omega _2} = {\omega _3} = 0$. For the points having phases satisfying Eq. (27), the second and third equations of Eq. (22) are not independent of each other, making the equation system become singular. Therefore, these points are the pitfall points for the least-squares algorithm proposed in Section 4.1. The pitfall points have special phases related to the pitch ratio ${\lambda }$ defined by Eq. (9) and the number of phase shifts*N*.

In a special situation that the beat period ${p_{\textrm{eq}}}$ is integer multiples of ${p_\textrm{I}}$ and ${p_{\textrm{II}}}$, there are *N* pitfall points appearing in a beat period with their phases being

Noting that the fringe patterns have noise in measurement practice, neither the phases exactly at the pitfall points nor the phases of their neighboring points can be effectively corrected using the least-squares algorithm in Section 4.1. For recognizing them, we set a threshold *T*. If a point has a phase satisfying

*T*is determined empirically, and its proper value is directly related to the noise level in the fringe patterns. Generally, in a high-noise situation, we will set a large threshold for making all the pitfall points reliably recognized.

#### 4.3 Phase error correction at pitfall points

For the pitfall points just recognized, we suggest using an interpolation-based method to correct their phase errors induced by fringe harmonics. For this purpose, we have three choices. Among them, the most straightforward method is to directly interpolate the phase map, but doing so may induce large errors when the phase map containing edges and discontinuities. An alternative solution is to interpolate the coefficient ${b_{N - 1}}$, with which the result may be badly affected by the nonuniform reflectivity and slope variations over the measured surface. By comparison, the preferable choice is to interpolate the amplitude ratio of harmonics defined by

After implementing the interpolation, we correct the fringe harmonics error using an iterative least-squares algorithm almost the same as that in Section 4.1. The difference is that the vector ${\mathbf \delta }$ of the equation system in Eq. (20) contains only three unknowns, i.e.,

We continue the simulation in Fig. 3. Using this interpolation-based technique, the pitfall points in Fig. 3(e) are recognized, and the phase errors at them are effectively corrected. The result is shown in Fig. 3(g) accompanied by its error curve below.

#### 4.4 Implementation procedure

The presence of the pitfall points makes the overall procedure for correcting the harmonics-induced phase errors somewhat complex. For making it clear and easy to follow, we summarize its steps here.

**Step 1**. Ignore the fringe harmonics for the moment, and calculate the background intensities, fringe modulations, and unwrapped phase map using the method in Section 3.

**Step 2.** Set a threshold for the inequality (30), and substitute the phases calculated in Step 1 into its left side in order to recognize the pitfall points.

**Step 3.** For the normal points whose phases not satisfying the inequality (30), correct their phase errors by using the iterative least squares algorithm presented in Section 4.1.

**Step 4.** For the pitfall points, correct their phase errors by using the interpolation-based algorithm suggested in Section 4.3.

Through this procedure, we use two-frequency phase-shifting fringe patterns to suppress the effects of fringe harmonics induced by adverse factors, simultaneously with phase-unwrapping, thus improving the measurement accuracy without capturing extra many fringe patterns or correcting the non-sinusoidal profiles of fringes.

## 5. Numerical simulations

In order to verify the feasibility of the proposed, we performed further numerical simulations. Assuming the number of phase shifts is three, Figs. 4(a) and 4(c) show two deformed fringe patterns having different frequencies (${p_\textrm{I}} = 512/9$ and ${p_{\textrm{II}}} = 64$ pixels) simulated under noise free condition. In them, the background intensities and fringe modulations are not uniform but Gaussian-shaped with a 50% decrease at corners. For introducing fringe harmonics, we assume the projector and camera have nonlinearities the same as those in Fig. 2. Using the synchronous detection algorithm [1], their wrapped phases are recovered and shown in Figs. 4(b) and 4(d). By employing the temporal phase-unwrapping technique in Section 3.2, the absolute phase map is obtained as shown in Fig. 4(e) with its carrier having been subtracted out. From it, subtracting the predefined phases gives the errors as illustrated in Fig. 4(f). In it, the maximum error is 0.1368 radians and the root mean square (RMS) error is 0.0783 radians. From Fig. 4(f), we observe that the phase errors appear as the ripple-like artifacts having frequencies three times higher than the fringes, being typically induced by the second order fringe harmonics. Before using the proposed technique, we employ the method in [32] as a comparison for correcting the phase errors. As we introduced in Section 1, the method in [32] is suitable for correcting phase-dependent errors purely induced by the projector nonlinearity. Its results are shown in Figs. 4(g) and 4(h) with their maximum and RMS errors being 0.0194 and 0.0045 radians, respectively. From Fig. 4(h), we observe the position-dependent residuals, implying that a pointwise algorithm is more suitable for compensating for such errors. When using the newly proposed technique, the result and its residual errors are shown in Figs. 4(i) and 4(j), respectively. It is evident from them that the ripple-like artifacts have been suppressed significantly by using the proposed technique. The maximum and RMS errors have decreased to 0.0065 and 0.0021 radians, respectively.

Furthermore, we implement the same simulations under different noise conditions, and use three-, four- and five-step phase-shifting algorithms for recovering the phase maps. The simulation results are summarized in Table 1. Using the traditional phase shifting method, increasing the number of phase shifts is helpful for improving the measurement accuracy. Theoretically the variance of the phase errors induced by random noise is independent of the phase value. It is proportional to the variance of noise and inversely proportional to the square of modulations multiplied by the number of phase shifts [31]. Simultaneously, using the *N*-step algorithm enables us to restrain the effects of fringe harmonics up to (*N*-2)th order [6,42]. These performances of the traditional phase-shifting algorithm are demonstrated in Table 1 by its third and fourth columns. With it, however, the accuracy of phase measuring depends only on the single fringe pattern sequence having higher frequency, and the low-frequency fringe pattern sequence does not contribute to the measurement accuracy. As a comparison, the fifth and sixth columns of Table 1 list the results of using the method in [32], which is based on an assumption that the camera is linear and the phase errors are purely induced by the projector nonlinearity. The last two columns show that the newly proposed method, using the same fringe patterns, enables achieving much higher accuracies than the existing methods. The reason is that the newly proposed method makes full use of the two-frequency fringe patterns, enabling restraining effects of the fringe harmonics up to (*N*-1)th order and being more effective than others in compressing the phase errors induced by noise. Simultaneously, the new technique effectively suppresses position-dependent phase errors because of its pointwise operation.

As described in Section 4, this proposed technique corrects the phase error at a point through an iterative procedure, no matter whether the point is an ordinary one processed using the iterative least-squares algorithm in Section 4.1, or it is a pitfall one which should be coped with using the interpolation-based method in Section 4.3. The convergence of the iterative procedure affects the computational efficiency. We investigate, through Fig. 5, the convergence of the new method under different noise conditions when the three-step algorithm is used. In it, the horizontal and vertical axes denote the number of iterations and the RMS values of the residual errors, respectively. It shows that the method quickly converges, after first several iterations, to almost 0 under noise free condition or to certain values depending on the noise levels.

## 6. Experiment and discussion

In this section, we use the proposed approach to measure practical objects for confirming its performance. The measurement system, as shown in Fig. 1, consists of a projector (PHILIPS PPX4010, 1280×720 pixels) and a camera (AVT Stingray F-125B, 1024×768 pixels) with a lens (KOWA LM12JC) having a focal length of 12 mm. In these experiments, we measure objects using three-, four-, and five-step phase-shifting algorithms successively. The used two sequences of fringe patterns have a horizontal orientation. Their fringe pitches are 48 and 45 pixels, respectively, so that their beat period, according to Eq. (10), is 720 pixels, exactly covering the full height of the projector image.

As a typically example, Fig. 6 shows the experimental procedure when measuring a circular object using three-step phase-shifting algorithm. In it, Figs. 6(a) and 6(c) are two fringe patterns having different frequencies, accompanied by their wrapped phase maps in Figs. 6(b) and 6(d), respectively. Figure 6(e) gives the unwrapped phase map obtained using the method in Section 3.2. After removing the carrier, this phase map reveals ripple-like errors induced by the fringe harmonics, and more essentially, by the nonlinearities of projector and camera. Similar to the simulation in Fig. 4, the fourth column of Fig. 6 gives the results of using the compensating method in [32], from which we observe that the artifacts in Fig. 6(f) are partially suppressed, but the residual errors are still noticeable because of the limitation of this method in coping with position-dependent errors. Using the newly proposed method in Section 4, we correct the effects of the fringe harmonics, and Figs. 6(i) and 6(j) shows the resultant phase map, before and after removing the carrier, respectively. It is evident that the ripple-like errors have been suppressed significantly by using the newly proposed technique.

We repeat measuring the same object by using four- and five-step phase-shifting algorithms. Figure 7 compares the measurement results before and after correcting the harmonics-induced errors. In it, the three rows, from top to bottom, are results obtained by using three-, four- and five-step phase-shifting methods, respectively. Here, we have converted the phases into depths according to their mapping relationship, which can be obtained by implementing a system calibration [46,47]. The leftmost column shows the depth maps reconstructed using the conventional technique introduced in Section 3 without correcting the harmonics-induced errors. Neighboring them on the right, the curves are vertical cross-sections along the same pixel column within the region of the plane outside the circular object. From these results, we observe that increasing the number of phase shifts is helpful for restraining the effects of fringe harmonics. If very few, e.g., three or four, fringe patterns are available, however, the reconstructed depth maps are badly ruined by ripple-like artifacts. The panels on the third and fourth columns, which are exactly aligned with those on the left ones, give the corresponding results obtained by using our proposed method. In them, the ripple-like errors have been suppressed significantly. To quantitively investigate the effectiveness of the proposed technique, Table 2 lists the RMS values of deviations of the measured depths along the cross-sections in Figs. 7(2) and 7(4). The data in this table reveal that the proposed technique allows us to compress the phase errors induced by fringe harmonics, especially when few phase shifts are used. Figure 8 investigates the convergence of the proposed method when using different number of phase shifts by plotting the RMS depth errors against the number of iterations. It is shown that, after first several iterations, the method converges to certain values depending on the noise levels and the number of phase shifts. As we mentioned previously, using the conventional phase-shifting technique, the *N*-step algorithm can restrain the effects of fringe harmonics up to (*N*-2)th order. When it is combined with two-frequency technique for the phase-unwrapping purpose, using the proposed method can restrain effects of the harmonics up to (*N*-1)th order. In Fig. 8, the errors with three-step algorithm converge to a value having somewhat large residuals since it is still affected by the third order harmonics. The errors with four- and five- step algorithms converge to values close to each other, meaning that noise in these cases has overwhelmingly large effects on the measurement results over the high order fringe harmonics.

To further verify the performance of the proposed method, we measure the second object which has a relatively more complex shape and a larger height. The measurement results are shown in Fig. 9, which has the same layout as that in Fig. 7. The similar phenomena are observed in Fig. 9, demonstrating that the proposed method works well in this case.

The experimental results have solidly demonstrated that using this newly proposed technique allows us suppress the errors induced by fringe harmonics significantly. As a software-based method, the proposed algorithm can be applied to measurement systems with different precision, helping them further improve measurement accuracy even breaking through hardware limitations. For this reason, the accuracy of this technique cannot be evaluated in isolation from a measurement system. Even though, from the simulation and experimental results, it is evident that the proposed method improves the measurement accuracy of a system by several times and even more, especially when few phase shifts are used.

With this method, however, there still exist some issues worth discussing. In addition to the fringe harmonics whose errors can be restrained by using the proposed technique, noise is the second most crucial error-inducing factor [48]. Random noise in a fringe pattern is induced by complicated physical factors, and can be simply modeled with the Gaussian distribution [49] according to the central limit theorem from probability theory. Usually we implement spatial low-pass filtering to an image, or capture multiple images and then average them temporally, in order to restrain the effects of noise. Note that this proposed method makes full use of two-frequency phase-shifting fringe patterns. It calculates the phase map form double fringe patterns, thus being more effective in compressing the noise-induced errors in comparison with the conventional techniques. Another error source is the illumination fluctuations of the projector, which induce ripple-like phase errors having the same frequencies with fringes. It can be corrected by use of fringe histograms [50]. Besides, the measurement accuracy is tightly related to the phase sensitivity of the measurement system [51]. The phase sensitivity mainly depends on the system geometry. If the system is fixed, the phase sensitivity is associated with fringe pitch and fringe orientations. In measurement practice, we prefer to use fringes having as small as possible fringe pitch in order to achieve a high phase sensitivity. On the other hand, very fine fringes may lead to difficulties in phase unwrapping, especially when using only two-frequency fringe pattern sequences. Using multi-frequency fringe patterns enables overcoming this problem. In fact, the principle of this technique is easy to expand to adapt to multi-frequency fringe patterns, so that the phase errors induced by fringe harmonics and by noise are eliminated by using multiple fringe pattern sequences having different frequencies. This proposed method can also be used for improving the measurement accuracy of other similar techniques like fringe deflectometry [17,18] and multifrequency interferometry [3–5].

Regarding the efficiency, the proposed technique suppresses the harmonics-induced phase errors at the expense of computational simplicity for data processing. The computational time with this technique also depends the available hardware resource. In these experiments, we used a laptop having a 2.80 GHz CPU (Intel Core i5-8400) and an 8GB RAM. Using the traditional phase-shifting algorithm without correcting the effects of fringe harmonics, it took about 1.3 seconds to recover a phase map of 590×510 pixels. Using the proposed technique, the algorithm converged within 10 seconds. Even though, this technique is still an efficient technique, because it does not require capturing extra more fringe patterns and can work well when the fringe pattern sequence of each frequency has few to three fringe patterns. It cancels time-consuming calibration procedure for correcting the non-sinusoidal profiles of fringes, and saves the time duration for image capturing in each measurement.

## 7. Conclusion

This paper has proposed a new algorithm of fringe harmonics elimination for two- or multi-frequency phase-shifting fringe projection profilometry. It makes full use of two fringe pattern sequences having different frequencies, in order to restrain the effects of fringe harmonics on the measurement results, simultaneously with phase-unwrapping. The simulation and experimental results have demonstrated that this technique enables effectively suppressing the ripple-like artifacts caused by the fringe harmonics, which are possibly sourced from various adverse factors, for example, the imperfections of device nonlinearities, insufficient defocus level when projecting binary patterns in dynamic measurement, and the camera saturation when measuring high dynamic range objects. Compared with the existing techniques, the proposed method offers several advantages. First, it directly estimates the coefficient of a high order harmonic from distorted fringe patterns, thus canceling time-consuming calibration procedure for correcting the non-sinusoidal profiles of fringes. Second, it works simultaneously with phase-unwrapping, thus not requiring capturing extra many fringe patterns. Third, it involves a quasi-pointwise operation, thus being suitable for coping with position-dependent phase errors and being helpful for protecting the edges and details of the measurement results from being blurred. As a software-based method, this newly proposed algorithm can be applied to measurement systems with different precision, helping them break through the accuracy limitations of hardware. Besides the fringe projection profilometry, this method can also be combined with other similar techniques like fringe deflectometry and multifrequency interferometry for improving measurement accuracies.

## Funding

National Natural Science Foundation of China (51975345).

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White, and D. J. Brangccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt. **13**(11), 2693–2703 (1974). [CrossRef]

**2. **J. C. Wyant, “Use of an ac heterodyne lateral shear interferometer with real–time wavefront correction systems,” Appl. Opt. **14**(11), 2622–2626 (1975). [CrossRef]

**3. **Y.-Y. Cheng and J. C. Wyant, “Two-wavelength phase shifting interferometry,” Appl. Opt. **23**(24), 4539–4543 (1984). [CrossRef]

**4. **Y.-Y. Cheng and J. C. Wyant, “Multiple-wavelength phase shifting interferometry,” Appl. Opt. **24**(6), 804–807 (1985). [CrossRef]

**5. **C. E. Towers, D. P. Towers, and J. D. C. Jones, “Optimum frequency selection in multifrequency interferometry,” Opt. Lett. **28**(11), 887–889 (2003). [CrossRef]

**6. **K. A. Stetson and W. R. Brohinsky, “Electro-optic holography and its application to hologram interferometry,” Appl. Opt. **24**(21), 3631–3637 (1985). [CrossRef]

**7. **J. J. Dirckx, W. F. Decraemer, and G. Dielis, “Phase shift method based on object translation for full field automatic 3-D surface reconstruction from Moiré topograms,” Appl. Opt. **27**(6), 1164–1169 (1988). [CrossRef]

**8. **T. Yoshizawa and T. Tomisawa, “Shadow moiré topography by means of the phase-shift method,” Opt. Eng. **32**(7), 1668–1674 (1993). [CrossRef]

**9. **Y. B. Choi and S.-W. Kim, “Phase-shifting grating projection moiré topography,” Opt. Eng. **37**(3), 1005–1010 (1998). [CrossRef]

**10. **V. Srinivasan, H. C. Liu, and M. Halioua, “Automated phase-measuring profilometry of 3-D diffuse objects,” Appl. Opt. **23**(18), 3105–3108 (1984). [CrossRef]

**11. **M. Halioua, R. S. Krishnamuthy, H.-C. Liu, and F.-P. Chiang, “Automated 360° profilometry of 3-D diffuse objects,” Appl. Opt. **24**(14), 2193–2196 (1985). [CrossRef]

**12. **S. Zhang, “Recent progresses on real-time 3d shape measurement using digital fringe projection techniques,” Opt. Lasers Eng. **48**(2), 149–158 (2010). [CrossRef]

**13. **H. Guo and B. Lü, “Phase-shifting algorithm by use of Hough transform,” Opt. Express **20**(23), 26037–26049 (2012). [CrossRef]

**14. **C. Zuo, S. Feng, L. Huang, T. Tao, W. Yin, and Q. Chen, “Phase shifting algorithms for fringe projection profilometry: a review,” Opt. Lasers Eng. **109**, 23–59 (2018). [CrossRef]

**15. **R. Zhang and H. Guo, “Depth recovering method immune to projector errors in fringe projection profilometry by use of cross-ratio invariance,” Optics Express. **25**(23), 29272–29286 (2017). [CrossRef]

**16. **J. S. Hyun and S. Zhang, “Enhanced two-frequency phase-shifting method,” Appl. Opt. **55**(16), 4395–4401 (2016). [CrossRef]

**17. **M. C. Knauer, J. Kaminski, and G. Hausler, “Phase measuring Deflectometry: a new approach to measuring specular free-form surfaces,” Proc. SPIE **5457**, 366–376 (2004). [CrossRef]

**18. **H. Guo, P. Feng, and T. Tao, “Specular surface measurement by using least squares light tracking technique,” Opt. Lasers Eng. **48**(2), 166–171 (2010). [CrossRef]

**19. **K. Creath, “Temporal phase measurement method,” in * Interferogram Analysis: Digital Fringe Pattern Measurement*, D. W. Robinson and G. Reid, eds. (IOP, Bristol, UK, 1993), pp. 94–140.

**20. **M. Chen, H. Guo, and C. Wei, “Algorithm immune to tilt phase-shifting error for phase-shifting interferometers,” Appl. Opt. **39**(22), 3894–3898 (2000). [CrossRef]

**21. **H. Guo and M. Chen, “Least-squares algorithm for phase-stepping interferometry with an unknown relative step,” Appl. Opt. **44**(23), 4854–4859 (2005). [CrossRef]

**22. **H. Guo, H. He, and M. Chen, “Gamma correction for digital fringe projection profilometry,” Appl. Opt. **43**(14), 2906–2914 (2004). [CrossRef]

**23. **C. R. Coggrave and J. M. Huntley, “High-speed surface profilometer based on a spatial light modulator and pipeline image processor,” Opt. Eng. **38**(9), 1573–1581 (1999). [CrossRef]

**24. **S. Kakunai, T. Sakamoto, and K. Iwata, “Profile measurement taken with liquid-crystal grating,” Appl. Opt. **38**(13), 2824–2828 (1999). [CrossRef]

**25. **T. Hoang, B. Pan, D. Nguyen, and Z. Wang, “Generic gamma correction for accuracy enhancement in fringe-projection profilometry,” Opt. Lett. **35**(12), 1992–1994 (2010). [CrossRef]

**26. **Z. Li and Y. Li, “Gamma-distorted fringe image modeling and accurate gamma correction for fast phase measuring profilometry,” Opt. Lett. **36**(2), 154–156 (2011). [CrossRef]

**27. **S. Ma, C. Quan, R. Zhu, L. Chen, B. Li, and C. J. Tay, “A fast and accurate gamma correction based on Fourier spectrum analysis for digital fringe projection profilometry,” Opt. Commun. **285**(5), 533–538 (2012). [CrossRef]

**28. **S. Zhang and S.-T. Yau, “Generic nonsinusoidal phase error correction for three dimensional shape measurement using a digital video projector,” Appl. Opt. **46**(1), 36–43 (2007). [CrossRef]

**29. **Z.-W. Li, Y.-S. Shi, C.-J. Wang, D.-H. Qin, and K. Huang, “Complex object 3D measurement based on phase-shifting and a neural network,” Opt. Commun. **282**(14), 2699–2706 (2009). [CrossRef]

**30. **B. Pan, Q. Kemao, L. Huang, and A. Asundi, “Phase error analysis and compensation for nonsinusoidal waveforms in phase-shifting digital fringe projection profilometry,” Opt. Lett. **34**(4), 416–418 (2009). [CrossRef]

**31. **F. Lü, S. Xing, and H. Guo, “Self-correction of projector nonlinearity in phase-shifting fringe projection profilometry,” Appl. Opt. **56**(25), 7204–7216 (2017). [CrossRef]

**32. **S. Xing and H. Guo, “Correction of projector nonlinearity in multi-frequency phase-shifting fringe projection profilometry,” Opt. Express **26**(13), 16277–16291 (2018). [CrossRef]

**33. **S. Xing and H. Guo, “Directly recognizing and removing the projector nonliearity errors from a phase map in phase-shifting fringe projection profilometry,” Opt. Commun. **435**, 212–220 (2019). [CrossRef]

**34. **G. A. Ayubi, J. A. Ayubi, J. M. Di Martino, and J. A. Ferrari, “Pulsewidth modulation in defocused three-dimensional fringe projection,” Opt. Lett. **35**(21), 3682–3684 (2010). [CrossRef]

**35. **Y. Wang and S. Zhang, “Three-dimensional shape measurement with binary dithered patterns,” Appl. Opt. **51**(27), 6631–6636 (2012). [CrossRef]

**36. **W. Lohry and S. Zhang, “Genetic method to optimize binary dithering technique for high-quality fringe generation,” Opt. Lett. **38**(4), 540–542 (2013). [CrossRef]

**37. **B. Li, Y. Wang, J. Dai, W. Lohry, and S. Zhang, “Some recent advances on superfast 3D shape measurement with digital binary defocusing techniques,” Opt. Lasers Eng. **54**, 236–246 (2014). [CrossRef]

**38. **D. Zheng, F. Da, Q. Kemao, and H. S. Seah, “Phase error analysis and compensation for phase shifting profilometry with projector defocusing,” Appl. Opt. **55**(21), 5721–5728 (2016). [CrossRef]

**39. **Z. Zhang, C. E. Towers, and D. P. Towers, “Time efficient color fringe projection system for 3D shape and color using optimum 3-frequency Selection,” Opt. Express **14**(14), 6444–6455 (2006). [CrossRef]

**40. **B. Salahieh, Z. Chen, J. J. Rodriguez, and R. Liang, “Multi-polarization fringe projection imaging for high dynamic range objects,” Opt. Express **22**(8), 10064–10071 (2014). [CrossRef]

**41. **S. Ordones, M. Servin, M. Padilla, A. Muñoz, J. L. Flores, and I. Choque, “Spectral analysis for the generalized least squares phase-shifting algorithms with harmonic robustness,” Opt. Lett. **44**(9), 2358–2361 (2019). [CrossRef]

**42. **H. Guo and M. Chen, “Fourier analysis of the sampling characteristics of the phase-shifting algorithm,” Proc. SPIE **5180**, 437–444 (2003). [CrossRef]

**43. **C. Zuo, L. Huang, M. Zhang, Q. Chen, and A. Asundi, “Temporal phase unwrapping algorithms for fringe projection profilometry: a comparative review,” Opt. Lasers Eng. **85**, 84–103 (2016). [CrossRef]

**44. **S. Xing and H. Guo, “Temporal phase unwrapping for fringe projection profilometry aided by recursion of Chebyshev polynomials,” Appl. Opt. **56**(6), 1591–1602 (2017). [CrossRef]

**45. **D. C. Ghiglia and M. D. Pritt, * Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software* (Wiley, 1998).

**46. **H. Guo, H. He, Y. Yu, and M. Chen, “Least-squares calibration method for fringe projection profilometry,” Opt. Eng. **44**(3), 033603 (2005). [CrossRef]

**47. **S. Xing and H. Guo, “Iterative calibration method for measurement system having lens distortions in fringe projection profilometry,” Opt. Express **28**(2), 1177–1195 (2020). [CrossRef]

**48. **K. Yan, Y. Yu, C. Huang, L. Sui, K. Qian, and A. Asundi, “Fringe pattern denoising based on deep learning,” Opt. Commun. **437**, 148–152 (2019). [CrossRef]

**49. **H. Guo, “A simple algorithm for fitting a Gaussian function,” IEEE Signal Process. Mag. **28**(5), 134–137 (2011). [CrossRef]

**50. **Y. Lu, R. Zhang, and H. Guo, “Correction of illumination fluctuations in phase-shifting technique by use of fringe histograms,” Appl. Opt. **55**(1), 184–197 (2016). [CrossRef]

**51. **R. Zhang, H. Guo, and A. K. Asundi, “Geometric analysis of influence of fringe directions on phase sensitivities in fringe projection profilometry,” Appl. Opt. **55**(27), 7675–7687 (2016). [CrossRef]