1. Introduction

In phase-shifting interferometry, miscalibration of phase shifter is a major factor decreasing measurement accuracy, and a number of algorithms have been developed for solving this problem. For example, the least-squares methods allow minimizing the measurement error induced by random phase shift errors [1,2], and averaging techniques and some optimized techniques can also compress the influence of phase shift errors [3–5]. A kind of more effective methods is to retrieve phases and phase shifts simultaneously, so that the errors caused by miscalibrations of phase shifters can be eliminated [6–11]. With these algorithms, however, capturing at least five interferogarms has been demonstrated to be the sufficient and necessary condition under which the phases and phase shifts are uniquely determined [11]. In recent years, some attempts have been made in order to break through this limitation, and several algorithms using three or four interferogarms have been proposed. These algorithms are useful in the situation that only fewer than five fringe patterns are available, for example in spatial phase-shifting interferometry [12] or when measuring a moveable profile so that the image capturing time is very tight. Among these algorithms, the algorithm of Cai et al. [13] exploits the statistics of fringes and is effective when the amplitudes of beams have been known prior to measurements. Also proposed by Cai et al. [14], the technique assumes the amplitude of the reference wave to be a constant over fringe patterns. The methods proposed by Xu et al. [15] and by Gao et al. [16] are also based on the same assumption that the amplitudes of the reference waves are uniform. The least-squares iterative algorithm of Wang and Han [17] can give satisfactory results when background intensities and modulations are evenly distributed. All these algorithms allow estimating phase shifts from only three interferograms, but the supplementary constraints in them, such as those concerning uniformities of the beam amplitudes or background intensities, imply that the number of unknowns decreases. These algorithms share a drawback that the supplementary restraints are usually very strong and not always satisfied in measurement practice, typically when the beam is not perfectly shaped or when directly a Gaussian beam is used.

Without any supplementary constraints, it is not easy to determine the unknowns such as phases and phase shifts when only three fringe patterns are captured, because the aforementioned sufficient and necessary condition is valid. In other words, the number of independent equations in this case is smaller than that of unknowns, thus making the equation system under-determined. A practically effective solution for this problem was proposed by Goldberg and Bokor [18]. It uses Fourier transforms to determine the phase shift between two interferograms with a carrier frequency. In the technique of Soloviev and Vdovin [19], properly tilting the reference mirror and using Hough transforms of the image differences can estimate the phase shifts between consecutive frames. The algorithm in Ref [20]. permits blindly estimating phase shifts by minimizing the cross power spectrum between the background and fringes. However, it is not stable enough when the problem is poorly conditioned. The reason is that the crosstalk introduced by phase shift errors between the calculated background and fringes is not linear but complicatedly nonlinear; and it is impossible to fully describe such nonlinear interactions by simply using a second order correlation (e.g. cross power spectrum). This reason makes the reliability of this algorithm significantly decreased.

From the aforementioned facts, we know that deriving a self-calibrating algorithm depending on only three interferograms in the absence of any preconditions or knowledge about interferograms still remains substantially challenging. To solve this problem, this paper presents, to the best of our knowledge, a novel blind self-calibrating algorithm, that minimizes the high order correlations between phases and amplitudes of the reference wave.

2. Principle

2.1 Reconstructions of phase and beam amplitudes

In phase-shifting interferometry, the intensities of the kth interferogram is often described with a function of the form

When K = 3, by defining c ₀ = |r|² + |o|², c ₁ = 2|r||o|cosφ, and c ₂ = −2|r||o|sinφ, Eq. (1) is restated as

Through Eqs. (3-8), the errors in phase shifts will introduce distortions in the reconstructed φ, |r|, and |o|. See Fig. 1 for example, (a) simulates a phase map and (b) is one of three interferograms generated with phase shifts being δ_k = {0, 1.7, 2.9} and the amplitude of the reference wave being Gaussian-shaped (with a 70% decrease in magnitude at the corners). The size of the interferogram is 256 × 256 pixels. Because the true phase shifts are unknown in measurement practice, we use their nominal values δ_k = {0, π/2, π} to recover φ and |r|. The resulting errors are displayed in (c) and (d), respectively, where the maximum phase error is 0.1676 radians and the RMS (root-mean-square) phase error is 0.1098 radians. From (c) and (d) we can observe that the shapes of artifacts in φ and |r| are similar to each other and to that of the fringes in (b), but they may have different spatial frequencies. This phenomenon typically reveals the effect of a nonlinear interaction between φ and |r| in the presence of phase shift errors. In other words, the phase shift errors introduce special correlations between φ and |r|. Because |r| is mainly determined by the light source and the optics in the reference arm, whereas φ depends on the optical path difference, naturally there should be no correlations between their distributions. For this reason, the actual phase shifts can be estimated by minimizing these correlations, but a problem arisen is how these correlations can be measured.

 

Fig. 1 Numerical simulation results. (a) Phase map. (b) One of three interferograms. (c) The errors in the reconstructed amplitudes of the reference wave and (d) the errors in the reconstructed phase map by using nominal phase shifts.

Download Full Size | PDF

2.2 Cross-bispectrum as a measure of correlations

As mentioned in Section 1, the algorithm from Ref [20] employs a cross power spectrum to estimate the correlations between the background intensities and fringes. This algorithm is not stable enough, because the cross power spectrum can only measure the second-order correlations between two signals and is blind to higher order ones introduced by nonlinear interactions [For example, there are two signals, i.e. f ₁(t) = cos(ωt) and its square f ₂(t) = cos²(ωt). The frequency of the first one is ω. The second signal can be restated as f ₂(t) = [1 + cos(2ωt)]/2, so its frequencies are 0 and 2ω. The spectra of f ₁(t) and f ₂(t) are not overlapped, thus their cross power spectrum equals 0. This fact means that the cross power spectrum cannot describe the nonlinear correlations between these two signals].

We overcome the above problem by using high-order correlations. Considering the cis function of φ

Continuing the simulation, 10 cross-sections shown in Fig. 2(a) are selected. They are uniformly distributed over the interferogram, thus being able to adequately represent its fringe features. Because it is possible that all fringes in an interferogram are parallel to a direction, we prefer to select cross-sections in both vertical and horizontal directions. For a foursquare interferogram like the one in Fig. 2(a), the diagonals, which have the same sizes in pixels with a vertical or horizontal cross-section, can be chosen supplementarily. Increasing the number of cross-sections may improve the accuracy of the estimated cross-bispectrum at the expense of computational burden. In addition, the sets of points randomly selected from the interferograms can also be utilized instead of these cross-sections.

 

Fig. 2 Numerical simulation results. (a) Selected cross-sections in the fringe patterns. (b) Magnitude of cross-bispectrum when (δ ₁, δ ₂) = (1.7, 2.9). (c) In the search regions, β as a function of (δ ₁, δ ₂). (d) The contracted search regions.

Download Full Size | PDF

When (δ ₁, δ ₂) = (1.7, 2.9), the magnitudes of $\widehat{B}$ are illustrated in Fig. 2 (b). Although $\widehat{B}$ is a measure for the third-order correlations, it is not convenient to compare. To solve this problem, we reasonably define a single-valued objective function of the form

2.3 Search procedure

It is almost impossible to minimize the objective function by using a parameterized method because of its prohibitive complexity, so we have to employ a search procedure to do it, calculating the objective function value for each pair of phase shifts and seeking the phase shifts by which the objective function has a minimum.

However, the possible range of each phase shift is the open interval (0, 2π) [or (-π, π)]. The computational burden of implementing a brute-force search procedure over the entire search ranges is not trivial. We must make some efforts to improve the efficiency. First, the inequality Σc ₀ ²≥Σ(c ₁ ² + c ₂ ²) [i.e. Σ(|r|² + |o|²)≥Σ2|r||o|] can be used as a criterion to exclude the phase shifts with very large errors before calculating the corresponding $\widehat{B}$. Figure 2(d) exemplarily demonstrates that the search regions are significantly contracted when those points dissatisfying the inequality just mentioned have been excluded. Second, $\widehat{B}$ is symmetrical regarding the line ω ₁ = ω ₂ as shown in Fig. 2(b). Thus we need only to compute the half-plane of Fig. 2(b) below ω ₁ = ω ₂. Third, noting that β(δ ₁, δ ₂) is symmetrical regarding the point (π, π), searching only the lower triangle when δ ₁<δ ₂ (or the upper one when δ ₁>δ ₂) can decrease the computational time to a half. Finally, we employ a cascade search strategy, with which several search cycles with increasing resolutions are performed in the cascade way. For example, in the first cycle we seek the minimum over the whole search region by using a search step of 0.02π, and in the second cycle we search for the minimum over a very small region centered at the previous searching result with a resolution 0.01π (halving the previous search step). In this way, the resolution can achieve 10⁻⁴ level after 7 cycles, and the computational time for one simulation is less than a minute.

3. Experiments

Experiments were carried out for further demonstrating the idea of the proposed algorithm. Figure 3 shows the measurement results of a silicon plate with a Fizeau interferometer. In this interferometer, the wavelength of the laser is 632.8nm, and an 8-bit CCD camera was employed to capture the fringe patterns. Figure 3(a) is one of three captured interferograms with the nominal phase shifts being δ_k = {0, 2π/3, 4π/3} and the size of each interferogram being 256 × 256 pixels. Figure 3(b) illustrates the calculated β(δ ₁, δ ₂), from which we see that the search regions have been contracted. In the low triangle, the minimum is located at the position (δ ₁, δ ₂) = (2.0044, 3.9617), which deviate from the nominal phase shifts.

 

Fig. 3 Measurement results of a silicon plate. (a) One of three interferograms. (b) Search regions. (c) The amplitudes of the reference wave and (d) the phase map, reconstructed by using the nominal phase shifts. (e) and (f) are similar to (c) and (d) but the used phase shifts are estimated with the proposed algorithm.

Download Full Size | PDF

The second row of Fig. 3 illustrates the reconstructed results by using the nominal phase shifts, with the left panel being the amplitude distribution of the reference wave and the right one being the phase map (the tilt has been removed). From both the reference wave amplitudes and the phases, we see the ripple-like artifacts typically induced by phase shift errors. In the reconstructed phase map, the RMS and PV (peak-to-valley) phase values are 0.1435 and 2.0150 radians, respectively. The third row is similar to the second row but is obtained by using the phase shifts estimated with the proposed algorithm. The artifacts appearing in the second row have been significantly compressed in the third row with the proposed algorithm. Correspondingly, the RMS and PV values of the phases are reduced to 0.1315 and 1.9003 radians, respectively.

Another example is presented in Fig. 4 , in which a grating is measured. Compared with the flat plate in Fig. 3, it has a more complex profile. Figure 4(a) is one of three interferograms generated using a Linnik interference microscope (The wavelength of the light source is 530nm and the camera has an 8-bit quantization resolution). The size of each interferogram is 512 × 512 pixels. The nominal phase shifts for the three interferograms are δ_k = {0, 2π/3, 4π/3}, but the phase shifts estimated with the proposed algorithm are (δ ₁, δ ₂) = (2.0593, 4.4251) which minimize the objective function β as shown in Fig. 4(b). Reviewing Section 2.3, the search regions are contracted by excluding the phase shift pairs that may make the recalculated intensities negative. In this example, however, the fringe modulation is relatively small because the measured grating has a much lower reflectance than the reference plate, and the contracted search regions in Fig. 4(b), as a result, are much larger than that in Fig. 3(b).

 

Fig. 4 Measurement results of a grating profile. (a) One of three interferograms. (b) Search regions. (c) The amplitudes of the reference wave and (d) the phase map, reconstructed by using the nominal phase shifts. (e) and (f) are similar to (c) and (d) but the used phase shifts are estimated with the proposed algorithm.

Download Full Size | PDF

Similarly to Fig. 3, the second and third rows of Fig. 4 illustrate the reconstructed results by using the nominal phase shifts and the estimated ones, respectively. In them, the left panels are the amplitude distributions of the reference wave and the right ones are the recovered phase maps (the tilt has been removed). The artifacts in the second row induced by the phase shift errors have also been significantly compressed in the third row.

4. Discussions

There are issues, for example the stability, worth discussion regarding this algorithm. When true phase shifts being close to multiples of 2π, very small phase shift errors will introduce large errors in the calculated phases, and the problem becomes poorly conditioned. The performance of this algorithm in this case have been evaluated with a numerical simulation, in which true phase shifts (δ ₁, δ ₂) = (0.01π, 1.99π), the estimated results are (0.0314, 6.2518) and the errors are very small; whereas the method from Ref [20] fails in estimating phase shifts because of the reason aforementioned. In fact, the proposed algorithm is effective even if only a fraction of a fringe is present when measuring a flat with a small curvature, but the accuracy may decrease. Tilting the reference mirror can introduce more fringes into the interferograms thus guaranteeing the accuracy.

The bispectrum methods are generally known to be insensitive to independent random noise. For the interferograms simulated in Section 2, the true phase shifts are (δ ₁, δ ₂) = (1.7, 2.9), and the estimated phase shifts in the presence of noise are listed in Table 1 , where the noise SD 0.001 corresponds to 0.255 gray levels in a 8-bit image and so forth. From this table we see that large noise (e.g. noise SD≥0.01) inversely affects the accuracy of the proposed algorithm, because the noise added to the fringe patterns, after a sequence of nonlinear computational operations with Eqs. (3-9), will not be independent and unbiased. For severely noisy interferograms, we can filter them before implementing the proposed algorithm, since the bispectrum methods are unaffected by a linear transform like low-pass filtering.

Table 1. Estimated phase shifts in the presence of noise

View Table

Compared with other self-calibrating algorithms, this method seems of a higher computational complexity, but it can deal with a much more difficult situation that other methods cannot (i.e. only three interferograms are available and any additional assumptions, like those concerning uniformities of beam amplitude or background intensities, are not satisfied). In addition, using a more advanced search technique may further improve the computational efficiency.

5. Conclusion

In conclusion, we have proposed a blind self-calibrating algorithm, which estimates phase shifts by minimizing the high-order correlations, between the reconstructed phases and the reconstructed amplitudes of the reference wave, introduced by phase shift errors. This algorithm has several advantages over others. First, it allows estimating the random phase shifts from only three fringe patterns. Second, supplementary assumptions (like uniform amplitude of beam) and knowledge about interferograms are not required. Third, it offers a high accuracy and a high stability, as well as a satisfactory computational efficiency.