## Abstract

We present an asynchronous phase-shifting demodulation approach based on the principal component analysis demodulation method that is robust to typical problems as turbulence, vibrations, and temporal instabilities of the optical setup. The method brings together a two-step and a phase-shifting asynchronous demodulation method to share their benefits while reducing their intrinsic limitations. Thus, the proposed approach is based on a two-fold process. First, the modulating phase is estimated from a two-step demodulation approach. Second, this information is used to compute weights to each phase-shifted pattern of the interferogram sequence, which are used in a novel weighted principal component demodulation approach. The proposed technique has been tested with simulated and real interferograms affected by turbulence and vibrations providing very satisfactory results in challenging cases.

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

## 1. Introduction

Phase-shifting interferometry (PSI) is an accurate optical metrology technique for measuring the modulating phase of interferograms. This information is very useful in optical metrology and can be applied, for example, to estimate the optical quality of an optical setup or device [1], or to image and analyze small biological objects as red blood cells [2,3], among other applications.

In standard PSI, a sequence of N interferograms is obtained having known and synchronous phase-shifts [4]. However, more recently have been proposed asynchronous demodulation methods as for example the iterative least-squares [5], PCA [6–8], the Euclidean Matrix Norm (EMN) [9], the UV factorization approach [10] or combination of these approaches [11] among others. These demodulation methods do not require to know the phase-shifts between interferograms as this information is determined internally by the approaches. In [5], the modulating phase and the phase-shifts are computed through a non-linear optimization process trigged by a first guess of the phase-shifts, which are typically randomly initialized. Moreover, the PCA demodulation method is a popular algorithm that provides a direct solution for both the modulating phase and the phase-shifts through the first two estimated principal components of the phase-shifted interferogram sequence. This approach is very fast and does not require any optimization process, initial guess of the phase-shifts, or additional constrains such as the background and modulation terms to be constant along the field of view. However, the PCA demodulation method requires the phase-shifts to be approximately equally distributed along the [0-π] range. To solve this limitation, recently, we proposed the UV factorization approach [10]. This method is based on an iterative factorization approach and non-linear boosting that expresses the interferogram sequence as the multiplication of two matrices, one depending on the phase-shifts and the other on the modulating phase. This method is robust to cases where the phase-shifted interferograms are unequally distributed in the [0-π] range and show spatially non-constant background illumination.

A typical practical problem in PSI, which affects to all approaches described above, is to determine the optimum number of phase-shifted interferograms. The use of a large number of patterns provides robustness against noise, while the use of fewer phase-shifted interferograms increases the robustness to uncontrolled mechanical vibrations, air turbulence, temperature changes or temporal instabilities of the optical setup. Thus, there is a tradeoff between the precision and accuracy of the reconstructed modulating phase caused by competing effects of the signal to noise ratio of the obtained modulating phase and the stability of the optical setup. In general, at least two phase-shifted interferograms are required to obtain the modulating phase. This fact motivates the different two-step demodulation approaches that have been proposed in the past [12–18]. These methods provide strong robustness to instabilities of the optical setup, as two interferograms can be obtained very fast. However, the reconstructed modulating phase could result very noisy and difficult (or even impossible) to unwrap and interpret. More general demodulation approaches focus on determining the ‘tilt-shifts’ to correct possible tilt errors caused by vibrations in the optical setup [19]. Other methods can determine the modulating phase in the presence of spatially and temporally varying phase-shifts, though, these approaches require a spatial carrier [20] or cannot retrieve unequivocally the phase map, which can be affected by a general non-linear function unless prior knowledge is available about at least one reference phase-shift [21].

In this work, we propose an approach that simultaneously takes advances of the two-step and multi-step demodulation approaches in terms of robustness to noise and instabilities of the optical setup, while is not affected by their intrinsic limitations. The method is based on providing weights to all different interferograms according to their similarity to a first estimated modulating phase obtained by a two-step demodulation approach. This phase map is assumed to be not affected by thermal and mechanical instabilities of the optical setup. Then, the weights and the phase-shifted interferograms are processed by a robust weighted principal component demodulation approach to compute an improved and robust to noise modulating phase.

## 2. Proposed method

The proposed method is based on a two-step process. First the approach estimates the modulating phase using a two-step demodulating approach. Any two-step method can be used, but in this work, we utilized the Kreis demodulation approach [12] because of its fast processing capacity and robustness to the background illumination. The resulting wrapped phase $({W\{{\tilde{\mathrm{\Phi }}({\mathbf r} )} \}} ){\; }$can be used to estimate weights of the *N* observed interferograms at the light of this estimated phase map. These weights are calculated projecting each interferogram into the two quadrature components computed from the wrapped phase $({W\{{\tilde{\mathrm{\Phi }}({\mathbf r} )} \}} ){\; }$determined by the two-step approach as

*X*, we compute the weighted interferogram sequence matrix ($wX$) as,

*W*a diagonal matrix given by

*D*is a diagonal matrix, and

*A*is the transformation matrix. This diagonalization process is practically performed by singular value decomposition. Then, the principal components are obtained directly by the Hotelling transform as We are only interested in the first two principal components with the biggest eigenvalues (first and second columns of

*Y*) that will correspond to a better estimation of ${Q_c}({\mathbf r} )$ and ${Q_s}({\mathbf r} )$ and the phase can be estimated from the interferogram sequence as

## 3. Results

We tested the proposed method with fringe patterns coming from simulations and real interferograms. In all cases, we compared the results obtained from a two-step demodulation approach (Kreis demodulation approach [12]), the PCA demodulation method [6] and the proposed weighted demodulation approach (wPCA). As we show in Code 1 (Ref. [22]), the source code for the presented methods reproducing all the results and figures shown here are freely available under the terms of an open-source software license.

#### 3.1 Simulations

We tested the performance of the proposed method with different simulations. In the first experiment, we generated 300 phase-shifted fringe patterns with random phase-shifts affected by additive white noise resulting into a signal to noise ratio (SNR) of ∼ 0.8. The background and modulation maps are $A({\mathbf r} )= x/50$ and $B({\mathbf r} )= \textrm{exp}({ - 0.5({{x^2} + {y^2}} )/{{10}^4}} )$, respectively, where *x* and *y* correspond to the normalized distance of every pixel to the image center and range between [-1, +1]. The image size is 450 × 450 px. The phase-shifts are randomly assigned and each interferogram is affected by a tilt aberration with random amplitudes obtained from a Gaussian distribution with standard deviation of 15 rad, which simulates the temporal instabilities of the optical setup. In Fig. 1(a) we show the first 20 interferograms used in the simulation.

As it can be seen from Fig. 1(a), the effect of the tilt aberration is in some case very strong, compromising severely the capacity of classical phase-shifting demodulation methods to compute an accurate estimation of the modulating phase. In Fig. 1(b) and 1(d) we show the reference and obtained wrapped and unwrapped phases using the Kreis, wPCA and PCA approaches. To reconstruct the unwrapped phases, we used the approach described in [23] because of its fast processing capacity. These images show that the PCA approach is completely unable to obtain an accurate phase map when the interferograms are affected by a strong varying tilt aberration. Indeed, from Fig. 1(d) can be seen that the unwrapped phase map cannot be reconstructed from the PCA wrapped map. Moreover, Fig. 1(b) shows that the Kreis approach can provide a correct phase map, although this wrapped phase is severely corrupted by noise. On the other hand, the result provided by our proposed wPCA approach provides an accurate and precise phase map reconstruction. The obtained root-mean-square error (rms) between the reference and the obtained unwrapped phase maps are 1.23, 0.36 and 7.85 rad for the Kreis, wPCA and PCA approaches, respectively. These results are also provided in Table 1. In Fig. 1(c) we show the obtained and used weights by our proposed wPCA approach. As can be seen from this figure, most of the interferograms are not taken into consideration in the phase demodulation as only 9% of the fringe patters show a weight larger than 0.2. The required processing time by the different approaches when running in a laptop using an Intel Core i7-8750H processor are given in Table 1.

In the second simulation, we generated similar interferograms than in the previous case with the only exception that here we reduced considerably the amplitude of the tilt aberration. Thus, in this case, the tilt aberration random amplitudes come from a Gaussian distribution with standard deviation of 2 rad. In Fig. 2 we show the results obtained. Figure 2 shows that the simulated interferograms are much less affected by the varying tilt aberration and consequently the weights are much higher in this case when compared with the previous simulation. Now 65% of the fringe patterns show a weight larger than 0.2. Additionally, from the phase maps shown in Figs. 2(b) and 2(d) we see that maps obtained from wPCA and PCA are both very similar to the reference map. The obtained *rms* and the required processing times computed when running in a laptop using an Intel Core i7-8750H processor are provided in Table 1. Note that even in this case, slightly affected by a temporal tilt aberration, the proposed wPCA can provide better results than the PCA demodulation method.

We also tested our proposed approach over phase-shifted fringe patterns corrupted by a highly varying defocus aberration, which amplitude is determined randomly from a standard normal distribution with standard deviation of 10 rad. The first 20 interferograms are given in Fig. 3(a). As before, we provide in Fig. 3(b) the obtained wrapped and unwrapped phases by the Kreis, wPCA and PCA approaches. Again, this figure shows that the PCA cannot compute a correct phase maps while the Kreis method provides a very noise reconstruction. Oppositely, our proposed wPCA approach provides accurate phase reconstructions. The obtained *rms* and the required processing times are provided in Table 1. In this case, the wPCA demodulation method provides weights larger than 0.2 to 17% of the used interferograms.

Finally, we tested our proposed approach over more complex phase-shifted fringe. To this end, we generated corresponding interferograms obtained from a modulating phase affected by astigmatism, spherical and comma aberrations, and in the same conditions than in the case shown before. The respective results are shown in Fig. 4 and Table 1. These results are equivalent to the ones shown in previous simulations.

#### 3.2 Experimental results

We tested our proposed methods with experimental interferograms obtained from a Mach-Zehnder interferometer in two different situations. The different datasets used here are available from [24, 25]. In the first case, the acquisition was affected by small instabilities while in the second case the optical setup was disturbed by more severe instabilities. As we show in Dataset 1 (Ref. [24]), the interferogram set is composed by 40 phase-shifted fringe patterns of size 1024 × 1024 px that are shown in Fig. 5(a). The interferograms have circular fringes obtained by the interference between plane and spherical waves. The spherical waves were obtained inserting a lens with large focal length in one of the arms of the interferometer. In Fig. 5(b) we show the obtained wrapped phases by the two-step Kreis approach and the wPCA and PCA phase-shifting demodulation methods.

From this figure we can see that all wrapped phases are similar, although both wPCA and PCA methods provide a phase reconstruction less affected by noise. Additionally, from Fig. 5(b) it is difficult to see differences between the wrapped phases obtained by wPCA and PCA as the modulating phase is large compared to the spurious phase coming from the instabilities of the optical setup. To compare both phase reconstructions, we have computed the Zernike decomposition of the wPCA and PCA obtained unwrapped phases. In Fig. 5(d), we show the subtraction of absolute values of corresponding Zernike coefficients using the Wyant notation by

*A*and

*i*refers to the fitted amplitude and Zernike index. As can be seen from Fig. 5(d), the larger difference between Zernike coefficients is for the Zernike coefficient corresponding to tilt aberration along the X axis (Zernike index #1 in the figure). In addition, this value is positive, indicating that the value provided by the PCA approach is more affected by tilt aberration. Moreover, the Zernike coefficients with index 8 and 15, corresponding to primary spherical and secondary spherical aberrations are larger for the wPCA phase, which is in good agreement with our previous knowledge about the shape of the interfering wavefronts. In Table 1, we provide the processing times required by the different approaches to process this data set.

In the second experiment, we modified the lens power and captured 300 phase-shifted fringe patterns with size 582 × 782 px. As we show in Dataset 2 (Ref. [25]), the first 40 interferograms of this image set are shown in Fig. 6(a). This phase-shifting interferograms were used to obtain wrapped phases by the two-step Kreis, the wPCA and PCA phase-shifting demodulation methods. As before, these phases were unwrapped using the approach described in (Zhao, Zhang et al. 2018). The results are shown in Figs. 6(b) and 6(c). As can be seen from these figures, in this case, both the Kreis and PCA methods provide reconstructed phases affected by artifacts, while the wPCA approach delivers an accurate absolute phase. Finally, in Fig. 6(d) we show the weights assigned to the different interferograms by the wPCA approach. In this case, the wPCA demodulation method provides weights larger than 0.2 to 23% of the used interferograms.

Finally, for the sake of completeness, we also compared all the obtained phase maps throughout the different simulations and experimental results, with the results provided by the Optical flow and the Gram-Schmidt two-step demodulation approaches. The results are shown in Table 1 and in Fig. 7.

## 4. Conclusions

In this paper, we introduce a method to improve the demodulation of phase-shifting interferograms affected by uncontrolled mechanical vibrations, air turbulence, temperature changes or temporal instabilities of the optical setup. A key advance of the proposed approach is that here we do not need to decide in advance about the number of patterns to be captured/processed, which have to be collected very carefully to avoid instabilities of the optical setup. In this case, the strategy is different and consists of collecting many interferograms assuming that a large number of them will be of low quality and will not be used by the demodulation approach.

The proposed method is based on combining the strengths of two-step and multi-frame phase-shifting demodulation methods. The former is robust against instabilities of the optical setup while the latter is robust to noise. A two-step approach (Kreis method in this case) is initially used to determine a first estimation of the modulation phase and quadrature components. These quadrature components are then applied to compute normalized interferogram weights by projecting the fringe patterns into them. The weights are then employed in a weighted principal component analysis demodulation approach (wPCA), where interferograms that are highly consistent with the previously estimated phase map contribute more than others that are not in good agreement. An obvious requirement of the proposed wPCA approach is the phase map determined by the two-step demodulation method to be accurate, while it can be severely affected by noise.

We test the wPCA with different simulations and experimental results. These results and the source code of the used methods are publicly available from [22] so all the results shown here can be reproduced independently. In the simulations, we test the performance of the wPCA with fringe patterns affected by a varying tilt and defocus aberrations. The experimental results show the applicability of wPCA to experimental interferograms affected by small instabilities of the optical setup. From our experiments, we show that the wPCA approach improves both the accuracy and precision of the computed phase maps when compared with the Kreis, PCA, GS and OF demodulation methods. In fact, in all the cases reported here our proposed approach clearly outperforms these approaches. Finally, wPCA shows a good trade-off between accuracy and required processing time.

## Funding

Ministerio de Ciencia, Innovación y Universidades (PID2019-108850RA-I00, Ramón y Cajal 2018 program (RYC2018-024087-I)).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

The data from this work is available as Dataset 1 [24], Dataset 2 [25].

## References

**1. **J. Vargas, N. Uribe-Patarroyo, J. A. Quiroga, A. Alvarez-Herrero, and T. Belenguer, “Optical inspection of liquid crystal variable retarder inhomogeneities,” Appl. Opt. **49**(4), 568–574 (2010). [CrossRef]

**2. **L. Xue, J. Vargas, S. Y. Wang, Z. H. Li, and F. Liu, “Quantitative interferometric microscopy cytometer based on regularized optical flow algorithm,” Opt. Commun. **350**, 222–229 (2015). [CrossRef]

**3. **Q. Wei, Y. Li, J. Vargas, J. Wang, Q. Gong, Y. Kong, Z. Jiang, L. Xue, C. Liu, F. Liu, and S. Wang, “Principal component analysis-based quantitative differential interference contrast microscopy,” Opt. Lett. **44**(1), 45–48 (2019). [CrossRef]

**4. **M. Servin, J. C. Estrada, and J. A. Quiroga, “The general theory of phase shifting algorithms,” Opt. Express **17**(24), 21867–21881 (2009). [CrossRef]

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

**6. **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]

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

**8. **J. Vargas and C. O. S. Sorzano, “Quadrature Component Analysis for interferometry,” Opt. Laser Eng. **51**(5), 637–641 (2013). [CrossRef]

**9. **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]

**10. **M. A. Escobar, J. C. Estrada, and J. Vargas, “Phase-shifting VU factorization for interferometry,” Opt. Laser Eng. **124**, 105797 (2020). [CrossRef]

**11. **J. Vargas, C. O. S. Sorzano, J. C. Estrada, and J. M. Carazo, “Generalization of the Principal Component Analysis algorithm for interferometry,” Opt. Commun. **286**, 130–134 (2013). [CrossRef]

**12. **T. Kreis and W. Jueptner, * Fourier transform evaluation of interference patterns: demodulation and sign ambiguity*, San Diego, ‘91 (SPIE, 1992), Vol. 1553.

**13. **J. Vargas, J. A. Quiroga, T. Belenguer, M. Servin, and J. C. Estrada, “Two-step self-tuning phase-shifting interferometry,” Opt. Express **19**(2), 638–648 (2011). [CrossRef]

**14. **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]

**15. **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]

**16. **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]

**17. **M. Trusiak and K. Patorski, “Two-shot fringe pattern phase-amplitude demodulation using Gram-Schmidt orthonormalization with Hilbert-Huang pre-filtering,” Opt. Express **23**(4), 4672–4690 (2015). [CrossRef]

**18. **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]

**19. **F. Zeng, Q. Tan, H. Gu, and G. Jin, “Phase extraction from interferograms with unknown tilt phase shifts based on a regularized optical flow method,” Opt. Express **21**(14), 17234–17248 (2013). [CrossRef]

**20. **J. Vargas, J. A. Quiroga, A. Alvarez-Herrero, and T. Belenguer, “Phase-shifting interferometry based on induced vibrations,” Opt. Express **19**(2), 584–596 (2011). [CrossRef]

**21. **R. Juarez-Salazar, C. Robledo-Sanchez, F. Guerrero-Sanchez, and A. Rangel-Huerta, “Generalized phase-shifting algorithm for inhomogeneous phase shift and spatio-temporal fringe visibility variation,” Opt. Express **22**(4), 4738–4750 (2014). [CrossRef]

**22. **J. Vargas, S. Wang, J. A. Gomez-Pedrero, and J. C. Estrada, “Code for Robust weighted principal components analysis demodulation algorithm for phase-shifting interferometry,” github, 2020https://github.com/1aviervargas/WPCA.

**23. **Z. Zhao, H. Zhang, Z. Xiao, H. Du, Y. Zhuang, C. Fan, and H. Zhao, “Robust 2D phase unwrapping algorithm based on the transport of intensity equation,” Meas. Sci. Technol. **30**(1), 015201 (2019). [CrossRef]

**24. **J. Vargas, S. Wang, J. A. Gomez-Pedrero, and J. C. Estrada, “Dataset 1 for Robust weighted principal components analysis demodulation algorithm for phase-shifting interferometry,” figshare2020 (https://doi.org/10.6084/m9.figshare.14045567.

**25. **J. Vargas, S. Wang, J. A. Gomez-Pedrero, and J. C. Estrada, “Dataset 2 for Robust weighted principal components analysis demodulation algorithm for phase-shifting interferometry,” figshare2020 (https://doi.org/10.6084/m9.figshare.14045585).