## Abstract

Focus on the phase reconstruction from three phase-shifting interferograms with unknown phase shifts, an advanced principal component analysis method is proposed. First, use a simple subtraction operation among interferograms, two intensity difference images are obtained easily. Second, set the center region of the data of intensity difference images to zero, and then construct a covariance matrix to obtain a transformation matrix. Third, two principal components of interferograms can be determined by the Hotelling transform and then phase can be calculated from the two normalized principal components by an arctangent function. By means of the simulation calculation and the experimental research, it is proved that the phase with high precision can be obtained rapidly by the proposed algorithm.

© 2015 Optical Society of America

## 1. Introduction

Phase-shifting interferometry (PSI) [1], which is a high precision optical phase measurement technique, has been widely used in various scientific fields, such as optical surface testing [2], deformation measurement, etc. Generally, to retrieve the measured phase, it is needed to capture at least three phase-shifting interferograms with the known phase shifts in PSI. Recently, several algorithms for phase retrieval and phase shift extraction from the phase shifting interferograms have been reported [3–22].

In [3], an advanced iterative algorithm (AIA) is proposed to determine the phase shift and the measured phase, while it is greatly time-consuming because of iterative operation. In [4], based on searching for the maximum and the minimum intensity of each pixel of interferograms in one phase shifting period, the corresponding background and modulation amplitude of interferogram can be determined, and then the phase can be retrieved by an arccosine function. However, the method consumes a lot of time, because the background and the modulation amplitude of interferogram are determined by the pixel by pixel calculation. In [5], a fast and simple phase shift extraction algorithm based on Euclidean matrix norm is presented, but it requires the number of interferograms is more than half phase shifting period, which means a lot of interferograms are needed. In [6], a phase shift extraction algorithm based on 1-norm is presented, while this method is also needed lots of interferograms.

Usually, in the temporal domain, the accuracy of phase retrieval is influenced by the surrounding environmental interference. Moreover, the quality of the captured interferogram is greatly related with the performance of phase shifter, which means in the practical PSI, it is still required to perform the phase retrieval from a few phase-shifting interferograms with unknown phase shifts. Recently, several two-step phase demodulation algorithms from two phase-shifting interferograms have been reported [7–10], and using these algorithms, the influence of environmental factors, such as the mechanical vibration, air turbulence, etc., can be reduced effectively. A two-step phase demodulation algorithm based on Fourier transform is introduced in [7], in which the measured phase can be directly retrieved by Fourier transform, but the accuracy of phase extraction is sensitive to the noise. In addition, a regularized optical-flow algorithm [8] is employed to retrieve phase. This method is robust against the additive noise and suitable for the phase extraction from interferograms with arbitrary phase shift. In [9], the Gram–Schmidt orthonormalization algorithm is used to extract the phase with high precision and rapid calculating speed. And in [10], it presents a phase retrieval method based on the extreme values of interference. This method is simple, fast and with high precision. However, all the two-step phase demodulation algorithms mentioned above, require that the background of the interferogram should be eliminated in advance by the filtering operation [11]. Therefore, eliminating the background of interferogram is an important factor to improve the precision of phase. Usually, the background of interferogram with slow changing can be eliminated effectively by filter in frequency domain, while this filtering method cannot work well when the background intensity with high frequency change.

Recently, Vargas *et al* [12–16] proposed a new method based on principal component analysis (PCA), which is able to extract the phase rapidly from randomly phase-shifting interferograms. Later, Xu et al extended the PCA technique to multiple beam interferometry and proposed a method which combining the PCA and the least squares method [18] to retrieve phase with low residual phase error. More recently, Zhang *et al* [19] extended the PCA technique to simultaneous dual-wavelength phase-shifting interferometry. However, in presenting PCA methods, the phase shifts of phase shifting interferograms should be well distributed in the range of [0, 2π] [12], because of which, the background of the interferograms can be eliminated effectively and Eq. (18) in [12] can be well satisfied. Furthermore, lots of phase shifting interferograms are needed in this PCA method. Besides, some algorithms based on intensity difference of phase shifting interferograms [6, 20–22] are presented to retrieve phase or phase shifts. These algorithms can eliminate the background of the interferogram effectively to retrieve phase with higher precision.

Clearly, among those methods mentioned above, some methods are time-consuming due to the iterative calculation [3] or the determination of background and modulation amplitude of interferogram by the pixel by pixel calculation [4], some methods have relatively low precision due to the unsatisfied background elimination of interferograms [7–10], and some methods are required lots of interferograms [5,6,12,13]. To solve these problems, we propose an advanced principal component analysis (APCA) method. In the APCA method, only three phase shifting interferograms with unknown phase shifts are needed. Generally, if the phase shift interval is not too small, the APCA method can be obtained good result. Moreover, there is no limitation to the background of interferogram.

## 2. Principle

In PSI, the intensity distribution of one arbitrary pixel of interferogram can be expressed as

where*n*and

*k*respectively represent the sequence number of phase-shifting interferograms and the pixel position, and in our case,

*n*= 1, 2, 3. For simplicity, we define ${\theta}_{1}=0$;

*a*and

_{k}*b*respectively represent the background and the modulation amplitude;

_{k}*φ*and

_{k}*θ*respectively represent the measured phase and the phase shift between interferograms; and

_{n}*K*represents the total pixel number of interferogram.

Here, we define the intensity difference images between the first interferogram and the *n*th interferogram as

According to the probability theory, the covariance matrix of two intensity difference images can be expressed as

*is transpose operator, and the covariance matrix*

^{T}*C*is a square and symmetric matrix, whose size is 2 × 2. For any real symmetric matrice $A=\left[\begin{array}{cc}a& c\\ c& b\end{array}\right]$,

*A*can be diagonalized by a matrice of $U=\left[\begin{array}{cc}1& 0\\ -c/a& 1\end{array}\right]$, namely,

*B*=

*U*･

*A*･

*U*, where

^{T}*B*is a diagonal matrix. Similarly, the covariance matrix C can be diagonalized aswhere

Then, the principal components can be determined by the Hotelling transform

*y*is a matrix of size 2 ×

*K*and two principal components

*y*

_{1}and

*y*

_{2}are corresponding to the first row and the second row of the

*y*matrix respectively. Substituting Eqs. (2), (3) and (8) into Eq. (10),

*y*

_{1}and

*y*

_{2}can be expressed as

To retrieve the measured phase distribution, *y*_{1} and *y*_{2} should be normalized by

Therefore, the measured phase can be retrieved directly by an arctangent function as

Importantly, to obtain the phase with higher precision, it is required that the approximations in Eqs. (9) and (15) can be well satisfied. On one hand, if there are more fringes in the interferogram, the approximations in Eqs. (9) and (15) can be better satisfied. On the other hand, the modulation amplitude of interferogram usually subjects to Gaussian distribution, thus the center region of the modulation amplitude is large, and the edge region is small. Moreover, for spherical wavefront, the phase change is faster comparing with that of the plane wavefront. If the center region of data of the interferogram is set to zero (the modulation amplitude is equivalent to zero), the approximations in Eqs. (9) and (15) can be better satisfied. In this paper, this characteristic is employed to retrieve the phase with higher precision.

The steps of the APCA method are as follows. First of all, use Eqs. (2) and (3) to obtain two intensity difference images${D}_{1,k}$ and${D}_{2,k}$. Second, set the center region of the data of the intensity difference images to zero, and expressed as${{D}^{\prime}}_{1,k}$and${{D}^{\prime}}_{2,k}$, then obtain a matrix ${Q}^{\prime}$, and subsequently construct a covariance matrix${C}^{\prime}$to obtain a transformation matrix${U}^{\prime}$. Third, two principal components of interferograms can be determined by the Hotelling transform, namely,$y={U}^{\prime}\cdot Q$. Finally, phase can be calculated from the two normalized principal components by an arctangent function.

## 3. Simulated results

To prove the effectiveness of the proposed APCA method, we have tested it with the simulated fringe patterns. As shown in Fig. 1, three simulated fringe patterns are generated, in which the background and the modulation amplitude are respectively set as *a*(*x*, *y*) = 100*exp*[-0.25(*x*^{2} + *y*^{2})] and *b*(*x*, *y*) = 120*exp*[-0.25(*x*^{2} + *y*^{2})]. The size of the fringe pattern is 512 × 512 pixels, and the phase is set as spherical wavefront *φ*(*x*, *y*) = 1.5*π*(*x*^{2} + *y*^{2})-<1.5*π*(*x*^{2} + *y*^{2})>, where <･> is the average operator, and −2.56 ≤ *x*, *y* ≤ 2.56. The phase shifts between fringe patterns are respectively 1 rad and 3 rad, and the noise is set as an additive Gaussian distribution with the signal-to-noise ratio (SNR) of 3%.

After obtaining two intensity difference images *D*_{1}* _{,k}* and

*D*

_{2}

*by Eqs. (2) and (3), we set the center region of the data of intensity difference images to zero, and expressed as${{D}^{\prime}}_{1,k}$and ${{D}^{\prime}}_{2,k}$. If we display*

_{,k}*D*

_{1}

*and ${{D}^{\prime}}_{1,k}$by the Matlab software, Figs. 2(a) and 2(b) are obtained. In Fig. 3, Fig. 3(a) shows the theoretical wrapped phase map, and Fig. 3(b) presents the corresponding wrapped phase map with the proposed APCA algorithm. In contrast, the wrapped phase map with the PCA algorithm is presented in Fig. 3(c), in which 30 phase-shifting interferograms are employed. In addition, the wrapped phase maps extracted by the GS3 algorithm [20] are shown in Fig. 3(d). Moreover, to compare the errors of the extracted phase and the calculating speed with different algorithms, the root mean square (RMS) errors of the difference between the theoretical phase and the extracted phase, as well as the processing time by the APCA, PCA and GS3 algorithms are respectively shown in Table 1.*

_{,k}In Table 1, we show the different RMS errors and the processing time which are obtained by a 3.2 GHz desktop. The RMS errors of the APCA, PCA and GS3 methods are respectively 0.038 rad, 0.008 rad and 0.042 rad. Meanwhile, the processing time of the APCA, PCA and GS3 methods is respectively 0.048 s, 0.482 s and 0.038 s. Obviously, from Table 1, it can be found that the RMS error of the proposed APCA algorithm is larger than that of the PCA algorithm, while it is smaller than that of the GS3 algorithm. And the processing time of the proposed APCA algorithm is much less than that of the PCA algorithm, and it is slightly slower than that of the GS3 algorithm.

Besides, we test the proposed method in the case of the interferogram with a significant higher noise. In our simulation, the noise is set as the additive Gaussian distribution with the SNR of 15%, and other parameters are the same as mentioned above. As shown in Fig. 4, three simulated fringe patterns with higher noise are generated, and the calculated RMS errors of the APCA, PCA and GS3 methods are respectively 0.216 rad, 0.039 rad and 0.217 rad. These results indicate that the APCA and GS3 method have almost the same precision, while the PCA approach will return better result, because a lot of frames are used and the noise can be reduced efficiently. In other words, the proposed method can perform well in the case of the interferograms with low noise and only requires 3 frames with unknown phase shifts.

Finally, we study the effects induced by different phase distribution on the measurement accuracy. As shown in Fig. 5, two sequences of simulated fringe patterns with different wavefronts are generated according to Eq. (1). Each sequence includes three frames with the size of 512 × 512 pixels, in which the background and the modulation amplitude are the same as mentioned above. The phase distributions, corresponding to the two sequences, are respectively set as the plane wavefront with *φ*(*x*, *y*) = 2*πx*-<2*πx*> and the complex wavefront with *φ*(*x*, *y*) = 2*π*Peaks(512)-<2*π*Peaks(512)>, where “Peaks” is the *peaks* function in Matlab. The noise is also set as the additive Gaussian distribution with the SNR of 3%. The calculated results are as follows: (1) for plane wavefront, the calculated RMS errors of the APCA, PCA and GS3 methods are respectively 0.037 rad, 0.008 rad and 0.037 rad; (2) for complex wavefront, the calculated results of RMS errors of the APCA, PCA and GS3 methods are respectively 0.037 rad, 0.008 rad and 0.037rad. The results indicate that the proposed method is suitable for different wavefronts. Furthermore, for both cases, the precision of the proposed method and the GS3 method is almost the same.

By analyzing the results, we find an interesting phenomenon. For spherical wavefront, the proposed method can be obtained with higher precision than GS3 method, while for plane wavefront and complex wavefront, the proposed method almost have the same precision as GS3 method. The main reason is that the modulation amplitude of interferogram usually subjects to Gaussian distribution, thus the center region of the modulation amplitude is large, and the edge region is small. Moreover, for spherical wavefront, the phase change is faster comparing with that of the plane wavefront. If the center region of the data of the interferogram is set to zero (the modulation amplitude is equivalent to zero), the approximations in Eqs. (9) and (15) can be better satisfied. The result indicates that the proposed is more suitable for the spherical wavefront.

## 4. Experimental results

Furthermore, the proposed APCA algorithm is also employed to perform the phase extraction with the experimental phase-shifting interferograms. As shown in Fig. 6, three phase-shifting interferograms are used for the phase extraction by the proposed APCA algorithm. In Fig. 7, Fig. 7(a) shows one of two intensity difference images, and Fig. 7(b) shows the intensity difference image whose center square region is set to zero, and the size of the center region is 200 × 200 pixels.

The wrapped phase map obtained by the conventional four-step method is used as the reference phase, as shown in Fig. 8(a). In addition, the wrapped phase maps with the APCA, PCA and GS3 algorithms are respectively shown in Figs. 8(b)–8(d). The corresponding RMS errors by the APCA, PCA and GS3 algorithms, as well as the processing time of phase extraction, are shown in Table 2. The RMS errors of the APCA, PCA and GS3 methods are respectively 0.036 rad, 0.018 rad and 0.040 rad. Meanwhile, the processing time of the APCA, PCA and GS3 methods is respectively 0.047 s, 0.504 s and 0.037 s. Clearly, like the simulation analysis, the RMS error with the proposed APCA algorithm is larger than that of the PCA algorithm, while the processing time of the proposed APCA algorithm is much less than that of the PCA algorithm. Moreover, the accuracy of phase extraction by the proposed APCA algorithm is higher than that of the GS3 algorithm, while the processing time is slightly slower than that of the GS3 algorithm.

Besides, phase shifting interferograms with phase distribution of plane wavefront are used to test the feasibility of the proposed method. As shown in Fig. 9, Fig. 9(a) shows one of interferograms with phase distribution of plane wavefront. The wrapped phase map obtained by the conventional four-step method is used as the reference phase, and shows in Fig. 9(b). The wrapped phase maps with the APCA, PCA and GS3 algorithms are respectively shown in Figs. 9(c)–9(e), and the corresponding RMS errors of the APCA, PCA and GS3 algorithms are respectively 0.166 rad, 0.040 rad and 0.166 rad. It is easy to find that the precision of the proposed APCA method and the GS3 method is almost the same.

In the proposed method, it is needed to set the center region of the data of the difference images to zero. Therefore, how to choose the size of the center region of the difference images is an important problem. Here we give an empirical approach. For spherical wavefront, the size of the center region of the difference image is near to the first center cycle fringes. For plane and complex wavefronts, the size of the center region of the difference image is very small, or even it is not necessary to set the data of center region of the difference images to zero.

## 5. Conclusion

In this paper, an APCA method is proposed to perform the phase reconstruction from three phase-shifting interferograms with unknown phase shifts. First, use the simple subtraction operation among interferograms, two intensity difference images without background can be obtained easily. Second, set the center region of the data of the difference images to zero and then construct a covariance matrix C to obtain a transformation matrix U. Third, two principal components of interferograms can be determined by the Hotelling transform. Finally, phase can be calculated from the two normalized principal components by an arctangent function. By means of the simulation calculation and the experimental research, it is proved that the measured phase can be obtained rapidly with high precision. Comparing the proposed APCA method with the PCA method, the PCA method is of higher precision, but it requires that the phase shifts of phase shifting interferograms should be well distributed in the range of [0, 2π]. Furthermore, lots of phase shifting interferograms are needed to retrieve phase. However, APCA method has no this limitation, only three phase shifting interferograms with unknown phase shifts are needed. Generally, if the phase shift interval is not too small, the APCA method can be obtained good result. Besides, the processing time with the proposed APCA algorithm is much less than that of the PCA algorithm. Comparing the proposed APCA method with the GS3method, the proposed APCA method is of the higher precision, while the processing time is slightly slower than that of the GS3 algorithm. Besides, the proposed APCA method is suitable for different wavefronts, while has a clear limitation with noisy interferograms. In a word, this proposed APCA method offers a powerful tool for the phase retrieval with three unknown phase shift interferograms.

## Acknowledgments

This work was supported by National Nature Science Foundation of China grants (51402148) and Shenzhen Nanshan Innovation Project (KC2014JSQN0011A).

## References and links

**1. **D. Malacara, M. Servín, and Z. Malacara, *Interferogram Analysis for Optical Testing* (the Chemical Rubber Company, 2005).

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

**3. **Z. Wang and B. Han, “Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms,” Opt. Lett. **29**(14), 1671–1673 (2004). [CrossRef] [PubMed]

**4. **Q. Hao, Q. Zhu, and Y. Hu, “Random phase-shifting interferometry without accurately controlling or calibrating the phase shifts,” Opt. Lett. **34**(8), 1288–1290 (2009). [CrossRef] [PubMed]

**5. **J. Deng, L. Zhong, H. Wang, H. Wang, W. Zhang, F. Zhang, S. Ma, and X. Lu, “1-norm character of phase shifting interferograms and its application in phase shift extraction,” Opt. Commun. **316**, 156–160 (2014). [CrossRef]

**6. **J. Deng, H. Wang, D. Zhang, L. Zhong, J. Fan, and X. Lu, “Phase shift extraction algorithm based on euclidean matrix norm,” Opt. Lett. **38**(9), 1506–1508 (2013). [CrossRef] [PubMed]

**7. **T. M. Kreis and W. P. O. Juptner, “Fourier transform evaluation of interference patterns: demodulation and sign ambiguity,” Proc. SPIE **1553**, 263–273 (1992). [CrossRef]

**8. **J. Vargas, J. A. Quiroga, C. O. 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]

**9. **J. Vargas, J. A. Quiroga, C. O. Sorzano, J. C. Estrada, and J. M. Carazo, “Two-step demodulation based on the gram-schmidt orthonormalization method,” Opt. Lett. **37**(3), 443–445 (2012). [CrossRef] [PubMed]

**10. **J. Deng, H. Wang, F. Zhang, D. Zhang, L. Zhong, and X. Lu, “Two-step phase demodulation algorithm based on the extreme value of interference,” Opt. Lett. **37**(22), 4669–4671 (2012). [CrossRef] [PubMed]

**11. **J. Antonio Quiroga and M. Servin, “Isotropic n-dimensional fringe pattern normalization,” Opt. Commun. **224**(4-6), 221–227 (2003). [CrossRef]

**12. **J. Vargas, J. A. Quiroga, and T. Belenguer, “Analysis of the principal component algorithm in phase-shifting interferometry,” Opt. Lett. **36**(12), 2215–2217 (2011). [CrossRef] [PubMed]

**13. **J. Vargas, J. A. Quiroga, and T. Belenguer, “Phase-shifting interferometry based on principal component analysis,” Opt. Lett. **36**(8), 1326–1328 (2011). [CrossRef] [PubMed]

**14. **J. Vargas, C. Sorzano, J. Estrada, and J. Carazo, “Generalization of the principal component analysis algorithm for interferometry,” Opt. Commun. **286**, 130–134 (2013). [CrossRef]

**15. **J. Vargas and C. Sorzano, “Quadrature component analysis for interferometry,” Opt. Lasers Eng. **51**(5), 637–641 (2013). [CrossRef]

**16. **J. Vargas, J. Carazo, and C. Sorzano, “Error analysis of the principal component analysis demodulation algorithm,” Appl. Phys. B **115**(3), 355–364 (2014). [CrossRef]

**17. **J. Xu, L. Sun, Y. Li, and Y. Li, “Principal component analysis of multiple-beam fizeau interferograms with random phase shifts,” Opt. Express **19**(15), 14464–14472 (2011). [CrossRef] [PubMed]

**18. **J. Xu, W. Jin, L. Chai, and Q. Xu, “Phase extraction from randomly phase-shifted interferograms by combining principal component analysis and least squares method,” Opt. Express **19**(21), 20483–20492 (2011). [CrossRef] [PubMed]

**19. **W. Zhang, X. Lu, C. Luo, L. Zhong, and J. Vargas, “Principal component analysis based simultaneous dual-wavelength phase-shifting interferometry,” Opt. Commun. **341**, 276–283 (2015). [CrossRef]

**20. **H. Wang, C. Luo, L. Zhong, S. Ma, and X. Lu, “Phase retrieval approach based on the normalized difference maps induced by three interferograms with unknown phase shifts,” Opt. Express **22**(5), 5147–5154 (2014). [CrossRef] [PubMed]

**21. **H. Guo and Z. Zhang, “Phase shift estimation from variances of fringe pattern differences,” Appl. Opt. **52**(26), 6572–6578 (2013). [CrossRef] [PubMed]

**22. **C. S. Guo, B. Sha, Y. Y. Xie, and X. J. Zhang, “Zero difference algorithm for phase shift extraction in blind phase-shifting holography,” Opt. Lett. **39**(4), 813–816 (2014). [CrossRef] [PubMed]