A novel method is presented to extract phase distribution from phase-shifted interferograms with unknown tilt phase shifts. The proposed method can estimate the tilt phase shift between two temporal phase-shifted interferograms with high accuracy, by extending the regularized optical flow method with the spatial image processing and frequency estimation technology. With all the estimated tilt phase shifts, the phase component encoded in the interferograms can be extracted by the least-squares method. Both simulation and experimental results have fully proved the feasibility of the proposed method. Particularly, a flat-based diffractive optical element with quasi-continuous surface is tested by the proposed method with introduction of considerably large tilt phase shift amounts (i.e., the highest estimated tilt phase shift amount between two consecutive frame reaches 6.18λ). The phase extraction result is in good agreement with that of Zygo’s MetroPro software under steady-state testing conditions, and the residual difference between them is discussed. In comparison with the previous methods, the proposed method not only has relatively little restrictions on the amounts or orientations of the tilt phase shifts, but also works well with interferograms including open and closed fringes in any combination.
©2013 Optical Society of America
The temporal phase-shifting interferometry (PSI) has been generally accepted as the most accurate technique for wave-front reconstruction based on automated interferogram analysis . PSI electronically records a series of interferograms while the reference phase of the interferometer is changed, so as to extract the phase information encoded in the variations in the intensity pattern of the recorded interferograms . However, the accuracy of PSI is limited by the phase-shifting uncertainties resulting from the imperfect conditions in the real testing environments, such as the miscalibration, non-linear responses, aging of the phase shifter , and the mechanical vibration .
As a result, much effort has been made to improve the PSI methods, to alleviate the influence of phase-shifting uncertainties on phase extraction accuracy [1, 3–15]. The papers [3–10] assume the phase steps are constant for all pixels in one interferogram but typically different between interferograms. However, sometimes such an assumption may be invalid, for instance, due to the unbalanced piezoelectric effect in the phase shifter or instability of the optical platform, the phase shifter may probably introduces non-negligible orientation errors during the shift [1, 11–15]; or as the extreme situations discussed in the paper , the strong environmental vibration will induce quite large tilt phase shift due to the relative motions between the test surface and the interferometer. In all those cases, the phase steps related to a specific interferogram are not longer constant for all pixels but vary with a tilt function across the field.
Several methods have been proposed to compensate the tilt-shift errors [11–15]. The method given in the paper  is based on a first-order Taylor series expansion of the phase shift errors, which iteratively update the parameters related to the tilt phase shifts until numerical convergence and the required accuracy are achieved. However, as the ratio of the tilt-shift errors to the unknown phase steps increases the first-order approximation will gradually lose the accuracy, then the compensation performance of the tilt-shift errors will get worse. In the papers [12, 14] the interferograms are divided into a set of small blocks, and the phase steps within an individual block are viewed as an unknown constant. In the paper  these constants are solved by using calculated contrast maps, while in the paper  they are optimized by the least-squares method . The size of the block in the papers [12, 14] should be reasonably defined beforehand. The size of the block would be very small if the tilt-shift amount is large. In that case, the estimation of tilt-shift parameters will be much coarse and even untruthful. The paper  proposes a method to extract the tilt-shift parameters from the set of zero-crossing points of the difference between the fringe intensities. However, the amounts and orientations of the tilt phase shifts are crucial to the performance of that method (for instance, if the amount of the tilt phase shift is less than 1λ, it is impossible to estimate the related tilt-shift parameters by that method), which may not be carefully controlled in the real testing environments. In the papers [1, 15] each interferogram is processed with the Fourier transform method (spatial carriers are assumed in the interferograms) so as to solve the related tilt phase shift parameters, and then the modulation phase can be extracted with the least-squares method. The paper  even has implemented experiments using a Mach-Zehnder interferometer, where the phase-shifts are introduced by induced vibrations without the need of any phase-shifter device. As the paper  includes affine registration preprocessing of the interferograms regarding for the camera movement due to strong vibration, it is expected to be applicable in a very hostile testing environment. However, if the interferograms are composed of closed fringes or the carrier frequency is inadequate, the methods given in the papers [1, 15] will be erroneous.
In this paper a novel method is proposed to cope with unknown tilt phase shifts. Firstly, a set of interferograms are captured, and the coarse wrapped phase shift distributions across the whole field between the selected interferograms are estimated in order, based on the regularized optical flow method . Secondly, those phase shift distributions are subsequently filtered and unwrapped, the parameters related to the tilt phase-shifting planes are determined with the help of an efficient 2-D frequency estimation method , and the sign ambiguity of the tilt phase shifts among interferograms is eliminated by comparing the coarse phase demodulation results obtained with the regularized optical flow method. Finally, accurate demodulation phase can be extracted from all the interferograms by the least-squares method, with all the estimated tilt phase shifts. The proposed method has some more flexibility than the previous methods [1, 11–15], i.e., it not only has relatively little restrictions on the amounts or orientations of the tilt phase shifts, but also works well with interferograms including open and closed fringes in any combination.
It should be pointed out that the extracted phase by the proposed method would have indetermination in the global sign if no prior knowledge about the phase shifts is available. However, such indetermination is not a particular problem of our method, because it is common to any asynchronous approach, if no information is given about the phase-shifts .
The rest of the paper is organized as follows: the principle of the proposed method is provided in detail in section 2, simulation and experimental results are reported in section 3 and section 4, respectively. Finally, conclusions are drawn in Section 5.
2. Principle13]. For instance, with the interferograms, the wrapped phase can be solved as :
In practice, if possible, it is always preferable to take more than three interferograms for the phase extraction. The reasons for it can be easily understood. On one hand, for pixels where one or both of the phase shifts are equal or close to integer numbers of , Eq. (2) will be invalid and the corresponding phase extraction result will be untruthful, as the equivalent number of interferograms is reduced for those pixels; On the other hand, as a rule of thumb, the PSI methods’ resistance ability to noise can be improved as more interferograms are involved.
Figure 1 shows the flow chart of the proposed method in this paper.
Obviously, the accurate estimation of all the phase-shifting plane parameters [see Eq. (1)] would be crucial to the performance of the phase extraction procedure. As we will see below, the proposed method would provide an accurate and fast solution for it.
Next, we will first briefly review the regularized optical flow method for the two-step interferometry, and then further introduce the procedures of the proposed method step by step. To be simple and intuitive, the descriptions in sections 2.1-2.4 are based on the case of two interferograms.
2.1 The regularized optical flow method for the two-step interferometry
Supposing there are two interferograms,, and their high-pass filtered background-suppressed versions are denoted as and . Then we have:, , where represents the phase shift distribution between the and , which is located at the phase-shifting plane: .
The paper  has put forward a regularized optical flow method to extract the wrapped demodulation phase from two interferograms with an arbitrary unknown but constant phase shift inside the range, with a single exception of , i.e. and are assumed. First, the fringe direction is obtained with a regularized optical flow method in an iterative way
2.2 Computation of the wrapped phase shift distribution between two interferograms
Actually, the phase extraction method given in the paper  can provide us with more information than the authors revealed.
On one hand, the phase shift between the two interferograms can be easily obtained with minor extension to that method, if its distribution is uniform across the field (i.e.). We find that the shape of the direction map is almost invariant as the phase shift changes, but the global sign of it is dependent on the range value of the phase shift. Specifically, if we denote the resultant direction map as (), and (), respectively, then we will have. As a result, we further haveEqs. (3)-(6). Meanwhile, from previous analysis, it is easy to deduce that: . Thus we can further get the following equationEqs. (9) and (10), we can conclude that by adding the two extracted phase results obtained with the regularized optical flow method, an estimation of the phase shift can be obtained which is specially wrapped inside the approximate range . The accuracy of the phase shift estimation can be improved by averaging it across the whole field, i.e., letting, where the functionis related to the average operation. Additionally, at some pixels the estimated phase shiftwill fall outside the range, mainly due to the noise in the interferograms and the truncation effect of the modulo operation. To compensate for it, we will add to these pixels where .
On the other hand, the phase shift distribution can also be analyzed as described above even if it is non-uniform across the whole field, since the regularized optical flow method is implemented in a local way. However, if the phase shift is equal to somewhere in the field, the computed wrapped phase shift distribution would include some flip lines. In this case, it should be correctly unwrapped inside the approximate range, which can be implemented as follows.
2.3 Unwrap the phase shift distribution
First, the map is filtered with a mean filter (the size of filter adopted in this paper is pixels) to generate a map denoted as, then a binary map denoted ascan be further produced as17, 18] are successively applied, to remove the spurious pixels and to fill small holes in the map. Then the map can be divided into several sub-regions according to the value of . In each sub-region of the map the extreme-value points are searched by scanning along the row or column direction (the selection of scanning direction is dependent on the orientation of the sub-regions), and then the corresponding flip line can be fixed by linear fitting with the coordinates of the extreme-value points. The estimated flip line related to the border sub-region would be rechecked to ensure its validity, by referring to the fitting residual error and comparing the slope of it with the counterparts in other sub-regions. Supposing totally estimated flip lines are found, the whole field will be divided intosub-regions, which will be labeled asin order by their spatial locations. Then a new binary map can be produced by letting
2.4 Estimation of the tilt phase shift parameters
From the expression related to the phase-shifting plane, we can deduce that the data in any row (column) of the complex signal is single-tone, with frequency (). Thus, we can estimate the phase-shifting plane parameters,from using an efficient 2-D frequency estimation method given in the paper . If the estimated results are denoted as and, then the estimated parameter d can be further obtained as
2.5 Elimination of the sign ambiguity among the estimated phase-shifting planes
Supposing there are totallyphase-shifted interferograms, then all the phase-shifting planes can be estimated in order with the procedures introduced in section 2.1-2.4. However, it can be noted that as to the specific phase-shifting plane, the selection of Eq. (12) or Eq. (13), will give rise to two different estimation results with the same phase-shifting distribution but with an opposite global sign. As a result, the estimated phase-shifting planes will have at mostpossible arrangements , and the majority of them will give erroneous phase extraction results, i.e., only the two arrangement which are closest to the true phase-shifting arrangement or its opposite version , will achieve the accurate estimation of the modulation phase or its opposite version . Here we put forward a simple way to find the interested arrangement. Considering the case of three phase-shifted interferograms , where the estimated coarse phase distributions [see Eq. (6)], the intermediate binary maps [see Eqs. (12) or (13)], and the estimated phase-shifting planes from the two-step phase-shifted interferograms , are denoted as, , and , respectively. Then we can define
2.6 Phase extraction with the least-squares method
Equation (1) can be rewritten as1, 4, 15] as followsEq. (20), . is the phase extraction result; while the estimated contrast parametersare dependent both onand the phase shifts related to the specific pixel, which can be taken as the reliability measure of the phase extraction result at that pixel.
3. Simulation results
As the accurate estimation of the phase-shifting plane parameters is crucial to the performance of the phase extraction procedure, a series of computation simulations have been carried out to verify the effectiveness of phase-shifting plane estimation by the proposed method. Some results will be given in the following paragraphs.
In the first simulation, the expressions for the two background-suppressed inteferograms are and , where , , and ; the units of and are both in radians, and the image sizes of and are both , i.e., and; andrepresent zero-mean additive Gaussian white noise, the standard deviation of which are both equal to 6; while is a built-in function of MATLAB, which is obtained by translating and scaling Gaussian distributions. The simulated phase-shifting plane is shown in Fig. 2(c), the tilt-shift amount across the whole field is 1.40λ(8.78rad); the estimated phase-shifting plane is shown in Fig. 2(i), and the residual PV (peak to valley) error and RMS (root mean square) error are as small as 0.019λ(0.012rad) and 0.004λ(0.026rad), respectively [see Fig. 2(j)]. The total processing time of this simulation isusing alaptop and MATLAB, and the main allocations of the runtime are as follows: aboutis spent associated with the regularized optical flow method (that method is run two times to estimate the phase shift); aboutis spent associated with the spatial imaging processing of the intermediate binary maps and , including the connected-component labeling, the morphological operations, and division of sub-regions operations; and about is spent on solving the parameters , by the frequency estimation method. We noted that the residual error distribution of the phase-shifting plane is linear proportional to the estimation errors of the parameters,and , thus the residual error will accumulate to its highest values at the border region of the field. If we divide the field into blocks and solve the parameter [see Eq. (15)] locally for each block, then the resultant residual PV error and RMS error related to this simulation will be decreased by 33% and 38%, respectively. However, such improvement ratios are case-dependent, so we will not have further discussions on it.
We also have evaluated the performance of our method by testing some other phase-shifting planes with different tilt-shift amounts and different orientations. The simulation parameters are all the same as in Fig. 2, except for the simulated phase-shifting plane. Specifically, four different phase-shifting planes are considered, where,, , and, i.e., the corresponding tilt-shift amounts across the field are ,,and, respectively. Actually, and are obtained by rotating anticlockwise in the x-y coordinate system with and, respectively. Herein, the definition ofis the same as shown in Fig. 2(c). The residual errors of the estimation results in PV and RMS measures are shown in Table 1. As can be seen from Table 1, all the residual PV errors and RMS errors are withinand, respectively.
The above simulations demonstrate that the proposed method can well estimated the tilt phase-shifting planes, and the estimation performance is robust to noise, the orientation and the amplitude of the phase-shifting plane. As explained before, the phase-shifting plane estimation result will face a global sign ambiguity problem. To reasonably evaluate the estimation performance of the proposed method, the simulation results given in Fig. 1 and Table 1 are all assigned with the correct signs.
4. Experimental results
For further verification of the feasibility of the proposed method, we have applied it to the experimental interferograms. Two set of experimental results will be provided below.
In the first experiment, the data of interferograms are obtained from the paper , which include totally nine real phase-shifted interferograms with unknown but constant phase shifts (the image size of them is). In this case, we only need to estimate the parameters. In Fig. 3 we show the reference phase obtained by the AIA method  using all the nine interferograms (b), as well as the phase extraction results by the AIA method but with only eight interferograms (c) and the proposed method (d). As the extracted phase by the AIA method would include a trivial constant bias, the residual errors of the latter two results after removing the bias would bein PV measure, in RMS measure [see Fig. 3(e)], andin PV measure, in RMS measure [see Fig. 3(f)], respectively. Particularly, we think the residual errors shown in Fig. 3(e) can be taken as an indirect measure of the non-ideality of the testing environment. Then we can conclude that in this experiment the proposed method has similar accuracy to the AIA method, as the residual errors shown in Fig. 3(e) and Fig. 3(f) are comparable in PV and RMS measures. The typical “double-frequency” fringe error shown in Fig. 3(f) is due to the discrepancies in phase shift estimations between the AIA method and the proposed method, i.e., as to the AIA method, the phase shifts are estimated in an iterative way and the information of different interferograms will be coupled in the least-squares sense during the iteration process, while for the proposed method, the estimation of phase shift between two interferograms is independent of the information encoded in other interferograms.
In another experiment, a circular flat-based diffractive optical element (DOE) with aperture of 100mm is tested, by measuring its transmitted phase. This DOE has quasi-continuous surface fabricated by the ion-beam etching technology, and it was made on the substrate of K9 glass. Firstly, the measurements are performed with a ZYGO Fizeau interferometer with the active vibration isolation workstation turned off. The recorded interferograms include moderate amount of tilt phase shift induced by environmental vibration and considerably large tilt phase shift by purposely rotating the reflective mirror (located behind the DOE) in the test arm around the axis of the PZT. As the rotation operation is implemented manually, in our experiment the interferograms are captured with a time interval of about 3s, to eliminate the contrast reduction in them due to the sudden “rotation” actions. We randomly pick out six consecutively captured frames to extract the phase component, where two of the frames are shown in Figs. 4 (a) and 4(b), and the wrapped phase extraction result with the proposed method is given in Fig. 4(c). On the other hand, we have also tested the DOE by the Zygo interferometer under vibration-isolated conditions with a calibrated PZT. The phase of the test surface extracted by the Zygo’s MetroPro software is shown in Fig. 4(d). Since noticeable lighting scattering would take place at the staircase locations of the DOE surface, the phase data are missing somewhere with the Zygo’s software. The difference between the phase extraction results by the proposed method and the Zygo’s software is presented in Fig. 4(e), which amounts toin RMS measures with the piston and tilt components removed; it must be pointed out that from the sense of surface test the relative difference between them would decrease nearly by half, as all the phase extraction results given in the following figures relate to the “double-pass” transmitted wave-front. In addition, the estimated results of the phase-shifting planes with the proposed method are shown in Fig. 5. It can be seen that the tilt-shift amounts between consecutive interferograms reach,,,,across the field, respectively; as to the accumulated phase shifts relative to the first interferogram, the highest tilt-shift amount across the field would be , corresponding to the phase shift distribution between the fifth and the first interferograms.
The residual error as shown in Fig. 4 is mainly due to the following factors: (1) the retrace errors; (2) the errors accompanied with the low estimated contrast parameters; (3) the fluctuations of the laser power; and (4) the dynamic variations of the testing environment, such as the air turbulence, which will be discussed in more detail. Firstly, as the wave-front of the DOE is far from the ideal plane wave, the introduction of tilt phase shifts will also give rise to non-negligible retrace errors. To demonstrate it, the DOE is measured at different status under vibration-isolated conditions. As shown in Fig. 6(c), the difference between the phase extraction results is mainly composed of tilt component, which equals to in the PV measure, and the resultant relative retrace error between them is shown in Fig. 6(d), which equals toin the RMS measure. As the tilt shift amount increases, the retrace error will be larger, and vice versa. As to our experiment, different retrace errors have been introduced in the interferograms we used, and the phase extraction result by the Zygo’s software also include some retrace error itself. Therefore, the residual error as shown in Fig. 4(e) will inevitably contain the contribution of retrace errors. Secondly, as shown in Fig. 4(e) and (f) the areas with high residual error values are found to be related with extremely low estimated contrast parameters. It is easy to understand that in such areas the accuracy of phase extraction result is much more susceptible to the noise in the interferograms, as well as the residual estimation errors of the phase-shifting planes. This problem can be alleviated by using more interferograms for the phase extraction. Thirdly, as the interferograms used in this experiment were captured during a relative long time (about 16s), the parameters,[see Eq. (1)] will be time dependent as a result of the fluctuations in laser power, which will introduce some error into the phase extraction result. Finally, the dynamic variations of the testing environment, such as the air turbulence will also make some differences between measurements.
As to the second experiment, since the tilt phase shift amounts related to the phase-shifting planes are considerably large, the AIA method will be failed; and the block-wise methods given in the papers [12, 14] also can hardly retrieve the correct phase, as in this case the interferograms should be divided into extremely many blocks with small sizes (probably more than one thousand blocks are required), so that the estimated phase shifts related to each small block are prone to be inaccurate. Besides, the methods given in the papers [1, 15] also can hardly be applied in this experiment, which rely on introduction of adequate spatial carriers into the interferograms. Additionally, we also found that if the tilt-shift amount across the field get too large (for instance, larger than), the direct estimation of the phase-shifting plane by the proposed method is prone to be inaccurate, due to the probable failure in the correct sub-region division of the intermediate binary map (as to the definition of the map, please referring to the section 2.3).
We have proposed a novel method for extracting phase distribution from interferograms with unknown tilt phase shifts. The proposed method can estimate the unknown tilt phase shift between two temporal phase-shifted interferograms, by extending the regularized optical flow method provided in the paper  with the spatial image processing and frequency estimation technology. With all the estimated tilt phase shifts, we can further obtain the phase extraction result with the least-squares method. Both simulation and experimental results have proved the feasibility of the proposed method. The proposed method is expected to be used in a testing environment with low frequency and high amplitude vibration, where costly and accurate phase-shifting devices are not longer required for steady-state measures.
Additionally, the proposed method has shown some more flexibility in comparison with the previous methods [1, 11–15]. On one hand, the proposed method has relatively little restrictions on the amounts or orientations of the tilt phase shifts. The methods proposed in [11, 12, 14] are mainly adaptable to the phase-shifted interferograms with relatively low tilt-phase shift amounts (the typical amounts reported in those papers are much less than 1λ); and the method proposed in  is quite sensitive to the amounts and orientations of the unknown tilt phase shifts among interferograms, particularly, if the amount of the tilt phase shift is less than 1λ, that method will fail. While, the method proposed in this paper can handle with the unknown tilt phase shifts with amount up to several wavelength, no matter the orientations of them. On the other hand, unlike the methods proposed in the papers [1, 15], the method proposed in this paper requires no spatial carriers in the interferograms, i.e., it can work well with interferograms including open and closed fringes in any combination.
We are grateful to Professor Javier Vargas at Centro Nacional de Biotecnología-CSIC, Spain, for his kind help about the phase-shifting methods and providing us with partial experimental data of interferograms used in this paper. The research was partially supported by the National Basic Research Program of China under grant No. 2011CB706701. Thanks also go to the anonymous reviewers for their valuable comments and suggestions.
References and links
2. D. Malacara, Optical Shop Testing, 3rd ed. (Wiley, 2007), Chap. 14, pp. 547–666.
3. R. Juarez-Salazar, C. Robledo-Sánchez, C. Meneses-Fabian, F. Guerrero-Sánchez, and L. M. Arévalo Aguilar, “Generalized phase-shifting interferometry by parameter estimation with the least squares method,” Opt. Lasers Eng. 51(5), 626–632 (2013). [CrossRef]
5. X. F. Xu, L. Z. Cai, Y. R. Wang, X. F. Meng, W. J. Sun, H. Zhang, X. C. Cheng, G. Y. Dong, and X. X. Shen, “Simple direct extraction of unknown phase shift and wavefront reconstruction in generalized phase-shifting interferometry: algorithm and experiments,” Opt. Lett. 33(8), 776–778 (2008). [CrossRef] [PubMed]
8. J. Vargas and C. O. S. Sorzano, “Quadrature Component Analysis for interferometry,” Opt. Lasers Eng. 51(5), 637–641 (2013). [CrossRef]
10. J. Vargas, J. A. Quiroga, C. O. S. Sorzano, J. C. Estrada, and J. M. Carazo, “Two-step interferometry by a regularized optical flow algorithm,” Opt. Lett. 36(17), 3485–3487 (2011). [CrossRef] [PubMed]
15. J. Xu, Q. Xu, and L. Chai, “Tilt-shift determination and compensation in phase-shifting interferometry,” J. Opt. A, Pure Appl. Opt. 10(7), 075011 (2008). [CrossRef]
16. S. Ye and E. Aboutanios, “Two dimensional frequency estimation by interpolation on Fourier coefficients,” in Proc. of IEEE Int. Conf. on Acoustics, Speech and Signal Processing, 3353–3356 (2012). [CrossRef]
17. R. C. Gonzalez, R. E. Woods, and S. L. Eddins, Digital Image Processing Using MATLAB (Prentice Hall, 2004).
18. L. G. Shapiro and G. C. Stockman, Computer Vision (Prentice Hall, Upper Saddle River, New Jersey, USA, 2001).