## Abstract

Principal component analysis phase shifting (PCA) is a useful tool for fringe pattern demodulation in phase shifting interferometry. The PCA has no restrictions on background intensity or fringe modulation, and it is a self-calibrating phase sampling algorithm (PSA). Moreover, the technique is well suited for analyzing arbitrary sets of phase-shifted interferograms due to its low computational cost. In this work, we have adapted the standard phase shifting algorithm based on the PCA to the particular case of photoelastic fringe patterns. Compared with conventional PSAs used in photoelasticity, the PCA method does not need calibrated phase steps and, given that it can deal with an arbitrary number of images, it presents good noise rejection properties, even for complicated cases such as low order isochromatic photoelastic patterns.

© 2016 Optical Society of America

## 1. Introduction

Photoelasticity is a widely employed technique in the field of stress measurement in transparent samples [1]. The basis of this technique lays in the induced birefringence that appears in a transparent material under stress. A 2D slab of material under stress behaves locally as a retarder plate with a retardation proportional to the difference of the principal stresses and also depedent on the orientation of the axes of those principal stresses. If we set the sample under stress inside a polariscope, there will appear an interferogram formed by two fringe patterns modulating each other: the isochromatics, associated with the local retardation and the isoclinics associated with the local orientation of the birrefrinence axes [1]. In principle, as with any other fringe pattern, the photoelastic interferograms are suitable to be analyzed employing one of the wide range of phase sampling algorithms (PSAs) available in the literature [2]. However, the inter-modulation of isochromatics and isoclinics makes necessary the adaptation of the standard PSAs [3]. Typically, the isoclinic interferogram is demodulated first, and then the orientation information obtained is employed in the isochromatic pattern demodulation. In the case of photoelasticity, the phase shifts are produced by the rotation of the polariscope elements, polarizers and/or quarter-wave retardation plates.

A typical PSA employed in photoelasticity is the eight-step technique reported by Quiroga et al [4]. In this case, a 4-step PSA is used for the isoclinics calculation and an 8-step PSA for the isochromatics. Both PSAs are based in the precise placement of the polariscope elements at given angular positions. Deviation of the actual angular positions from the planned ones are translated into detuning errors for the PSA [2–5]. This kind of error can be eliminated if we use self-tunable or wide-band PSAs for the fringe pattern demodulation. In both techniques, it is assumed that the phase step between interferograms is unknown but constant. For the self-tunable PSAs [2] the local phase step is estimated and then passed to a tunable PSA [2] or [6,7]. Interestingly enough, the first reported PSA was a self-tunable method [8]. On the other hand the wide-band PSAs are designed to have low detuning error for a given phase steps range [5,6].

When using PSAs, an important issue is the possibility to use an arbitrary number of images. In general the bigger the number of images used by a PSA the better are the noise rejection properties and the tolerance for other sources of errors as phase step detuning and nonlinarites (interferogram over or under exposure, for example).

In the photoelastic case, isoclinics can be processed using standard PSAs with the number of phase steps depending on the application. However, isochromatics calculation depends on a fixed number of polariscope configurations (typically 6 to 8) that preclude the use of a PSA with an arbitrary number of images.

In this work, we present the application of the Principal Components Analysis Phase Shifting Algorithm (in the following PCA) technique to the case of photoelastic fringe patterns. PCA is a self-calibrating PSA method [9,10] that can demodulate the phase from a set of interferograms with unknown phase steps. In the photoelastic case, the self-calibrating character means that we do not need to set the polariscope elements at precise angular positions. We can use the number of positions (or steps) that are more convenient for our experimental setup, even without knowing the actual angular position with great accuracy. In addition, PCA is inherently resistant to harmonics or/and other issues such as over or under exposure that might arise from image capture. If necessary, the proposed algorithm can be improved using iterative least squares techniques [11,12].

The paper is organized as follows: in next section we present the theoretical foundations of the application of PCA to photoelasticity, afterwards we will present and discuss the results obtained when we have applied this technique to a pair of photoelastic fringe patterns and, finally, we will give the conclusions to end the paper.

## 2. Theoretical foundation

In photoelasticity, the isoclinic and isochromatic fringe patterns modulate each other. Therefore, in order to apply any PSA technique, including PCA, it is necessary to demodulate in two different steps the isoclinics and the isochromatic fringe patterns. Thus, we will present first the procedure to demodulate the isoclinic fringe pattern using the proposed PCA algorithm and, afterwards, we will make use the information obtained from the isoclinics fringe pattern to demodulate the isochromatics one. However, we will first summarize the PCA method for phase demodulation [10].

#### 2.1 PCA analysis of phase shifted interferograms

In phase sampling interferometry (PSI), the main objective is the demodulation of the modulating phase from a set of *K* phase-stepped interferograms given by

*k*-th interferogram and ${I}_{k}\left(r\right)$is the irradiance distribution. The physical interpretation of$\phi $depends on the experimental setup. In the case of photoelasticity, there are two modulating phases, one related with the principal stresses difference and the other with the direction the first principal stress ${\sigma}_{1}.$

We can always rewrite Eq. (1) as

The PCA method consists in three steps. First, every phase-shifted interferogram is stacked columnwise in a$N\times 1$column vector, ${I}_{k}\left(n\right){\text{\hspace{0.17em}}}_{n=1,2,\mathrm{...},N}$, being *N* the total number of pixels in the image. With the vectorized images, we construct a $K\times N$matrix **X** as

**X**is a vectorized interferogram. In the second step, we calculate the$K\times K$covariance matrix for the interferograms aswhere T is the transponse operator and${m}_{X}$a matrix with the same size as

**X**and with all elements of each column equal to the mean value of the respective column of

**X**, ${m}_{X}\left(n,j\right)=\frac{1}{K}{\displaystyle \sum _{k}{I}_{k}\left(n\right)\approx a\left(n\right)\text{\hspace{0.17em}}}\forall j$. Therefore, each element of the columns of ${m}_{X}$is an estimation of the DC term so that its subtraction is, indeed, a DC filtering process. From construction, ${C}_{X}$ is symmetric and it is always possible to find a linear transformation

**A**, (a $K\times K$matrix) that diagonalizes it, so that the transformed covariance matrix ${C}_{Y}=A\xb7{C}_{X}\xb7{A}^{T}$is diagonal. In the literature, the linear transformations

**A**is also known as Karhunen-Loève or Hotelling transforms [10].

The third step is to transform the PSI matrix **X** using the Hotelling transform to obtain the transformed PSI interferograms

**Y**is a principal component. The principal components are transformed interferograms with zero mean and diagonal covariance matrix (zero correlation). It can be demonstrated [13] that for a pure cosine signal like the one given in Eq. (2) the two first rows of

**Y**are the vectorized versions of${g}_{1}$and${g}_{2}.$Finally, the phase is calculated from the principal components asHowever, when dealing with actual interferograms, usually we do not know the correct sign of the global phase because the cosine and sine signals are arbitrarily assigned to the first and second principal components and, also, because the sign for every principal component is arbitrary. This issue can be solved if we have additional information about the phase steps. For example, if they increase monotonically, we could still recover the correct sign of the global phase [12]. For additional details about the PCA algorithm the reader can refer to [9,10] and [13].

#### 2.2 Isoclinics demodulation

In our experiments, we have used a motorized polariscope composed of a linear polarizer, two quarter-wave plates and a second linear polarizer that is denominated analyzer. In this system, the 2D photoelastic sample is placed between the quarter-wave plates and observed in transmission. For the isoclinics calculation we set the polariscope with the polarizers axis and quarter-wave plates fast axes parallel to the vertical Y direction, this configuration is called linear bright field configuration (LBF), and it would be equivalent to set the sample between two linear polarizers as shown in Fig. 1.

If we place a two-dimensional photoelastic sample between the quarter-wave plates, this is equivalent to introducing a retarder plate,${R}_{\theta}\left(\delta \right),$with a spatially variable fast axis direction$\theta \left(r\right)$and retardation$\delta \left(r\right),$the latter given as

where λ is the wavelength used to illuminate the sample,*C*the so-called photoelastic constant, ${\sigma}_{1,2}$the principal stresses and

*d*the sample thickness. The isoclinic angle

*θ*gives the direction of the 1st principal stress${\sigma}_{1}$with respect the X axes. The direction of the second principal stress${\sigma}_{2}$is at $90\xba$ with respect to that of${\sigma}_{1}.$

In the LBF configuration, the normalized output intensity is given by [1]

For the isoclinics, the phase sampling method consists in the rotation of the four elements of the polariscope by a set of *K* angles${\beta}_{k},$while keeping the sample fixed (See Fig. 1). One advantage of the PCA is that these angles can be unknown. For an angle ${\beta}_{k}$ the output intensity for the LBF configuration is given by [14]

*θ*from $W\left[4\theta \right],$we would not know whether

*θ*is the orientation of${\sigma}_{1}$or that of${\sigma}_{2}.$If we were only interested in the orientation of the local retardation axes, this uncertainty will not be problematic. The reason is that it can be interpreted as a global $\pi /2$ rotation of the reference system, meaning that the ${\sigma}_{1,2}$ orientations are interchanged, but the local orientation of both retardation axes is correct.

In addition, according to the former discussion about the PCA technique, there is a global sign uncertainty in the value of$W\left[4\theta \right]$returned by Eq. (12). In digital photoelasticity, the isoclinic parameter of interest is$W\left[2\theta \right]$because this magnitude determine the orientations of the two birefringence axes associated with${\sigma}_{1}$and${\sigma}_{2}$without ambiguity. The sign indetermination in $W\left[4\theta \right]$can be interpreted as a *π* phase shift and translates as a global$\pi /2$shift in W[2θ], that can transform $\mathrm{cos}\left(2\theta \right)=\mathrm{cos}\left(W\left[2\theta \right]\right)$ into $\mathrm{sin}\left(2\theta \right)=\mathrm{sin}\left(W\left[2\theta \right]\right)$and vice-versa.

#### 2.3 Isochromatics demodulation

Once the isoclinics have been measured, the next step is the demodulation of the isochromatic fringe pattern. To do so, we have to change the configuration of the linear bright field polariscope depicted in Fig. 1. For the isochromatics calculation, we are going to illuminate the sample with circularly polarized light. This will fix the angular positions of the Polarizer and first Quarter-wave plate to be at $45\xba$. On the other hand, the second quarter-wave plate and the analyzer can have any orientation. This configuration is shown in Fig. 2. Thus, if we designate ψ as the angle formed by the slow axis of the quarter wave plate with the X-axis and φ as the angle formed by the transmission axis of the analyzer and the X-axis, the normalized irradiance collected by the CCD camera can be written as [1]:

Here, direct application of PCA (or any PSA) will not yield the two principal components ${g}_{1}=\mathrm{cos}\delta $and ${g}_{2}=\mathrm{sin}\delta $. As we can see in Eq. (13), using circular illumination it is possible to isolate the $\mathrm{cos}\delta $ term (for example using $\left(\psi -\phi \right)=\pi /4$). However, the term $\mathrm{sin}\delta $will always be modulated by the isoclinics through the $\mathrm{sin}2\left(\varphi -\theta \right)$ term. That means that we need a special PSA for the isochromatic demodulation that could take this into account [4]

For the PCA, the solution we found is to use two different configurations of the circular polariscope depicted in Fig. 2. In one configuration we will set the second quarter wave-plate horizontal,$\phi =0$and in a second configuration we will lay the second quarter-wave plate at $45\xba$.

If we set the quarter-wave axis horizontal φ = 0 and we rotate the analyzer by a set of *K* values, ${\psi}_{k}{\text{\hspace{0.17em}}}_{k=1,2,\mathrm{...},K}$, the output intensity for each analyzer position will be

Next we orientate the second quarter-wave plate at $45\xba$, setting $\phi =\pi /4$ and rotate again the analyzer by a set of *K* angles, ${\psi}_{k}{\text{\hspace{0.17em}}}_{k=1,2,\mathrm{...},K}$. In this case, the output intensity for each analyzer position will be

Unwrapping of Eq. (19), together with the fact that, by definition,$\delta \ge 0,$will yield the isochromatic retardation *δ* at every location. Formally, Eq. (18) is the same result obtained in [4] but using an arbitrary number of (unknown) analyzer orientations, instead of the eight orientations used in [4]. The use of arbitrary and unknown orientations for the analyzer has two advantages: First is is not necessary to employ precise and expensive motorized components and, second, the possibility to reduce the SNR of the measured retardation by increasing the number of phase stepped images. This latter is a very useful feature when measuring samples with low order isochromatic fringes, as is the case of low stress, residual stresses or materials with a high photoelastic constant. In those cases, the sample is almost isotropic for all locations and the demodulation of the isoclinic and isochromatic phases can be seriously affected by noise.

In principle, Eqs. (12) and (19) solve the problem of the calculation of isoclinics and isochromatics from an arbitrary number of phase stepped images. However as we have mentioned, the way in which PCA operates introduces two uncertainties. For the isoclinics, we have a $\pi /2$ ambiguity in the $W\left[2\theta \right]$ phase, that can produce a swap between$\mathrm{cos}2\theta $and$\mathrm{sin}2\theta .$ For the isochromatics, we have a global sign uncertainty in the principal components given by Eqs. (15) and (17). These can be solved by imposing some a-priory knowledge on the phase steps, for example monotonicity [12]. In this work, we propose an alternative heuristic solution consisting in a thorough search of the best sign combination for the phasor described in Eq. (18). For this purpose, we calculate from the isoclinic phase $W\left[2\theta \right]$computed previously from $W\left[4\theta \right]$, Eq. (12), and from the isochromatic principal components of Eqs. (15) and (17) the following eight phasors

*m*gives the average amplitude of the phasor$z(r)$in all the positions inside the region of interest. In Eq. (21),$z*(r)$stands for the complex conjugate of $z(r).$ After computing

*m*for all the eight phasors, we will select the sign combination which yields the highest value of the average modulation.

The idea behind this procedure is that any sign error will greatly decrease the modulation of the phase map in those regions where the incorrect sign introduces errors. If${z}_{\mathrm{max}}$is the phasor with the maximum average modulation *m*, then the isochromatic phase will be calculated as

## 3. Experimental results

In order to test the proposed method, we have carried out a photoelastic experiment in which we have measured the isoclinics and isochromatics for two different objects: a planar disk in diametrical compression and a plastic ophthalmic lens which have been cut oversized so it is pressed by its metallic frame. We have computed the isochromatic and isoclinics maps for these samples using 20 images for determining the isoclinics and other 20 images to get the isochromatics. As an alternative, we have also demodulated the fringe patterns using the PSI algorithm proposed in reference [4] taking 4 images for retrieving the isoclinics and another 8 images for determining the isochromatics. All the measurements have been made using a motorized polariscope with a digital CCD camera for a fully automatic acquisition, and a dual white/monochromatic light source.

Regarding the first sample, the planar disk compressed along a diameter, we show in Fig. 3(a) the wrapped phase map corresponding to the spatial distribution of the angle 4θ (directly related to the isoclinics fringe pattern through the orientation angle θ) obtained after applying the proposed PCA algorithm while in Fig. 3(b) we have depicted the same phase map obtained through the application of a standard 4-step phase shifting algorithm [4]. As it can be seen in this figure, we have obtained a better estimation of the isoclinics orientation, particularly for the horizontal meridian of the disk which is the direction perpendicular to the direction in which the stress is applied to the sample. For the loaded disk, the isoclinic phase jump should be an horizontal straight line. So any deviation form straightness would indicate errors in the phase demodulation. As we can see in Fig. 3, for the PCA solution, these deviations in the phase jumps are less conspicuous that those presented in the PSA case.

Similarly, we have depicted in Fig. 4 the wrapped phase maps corresponding to the isochromatic fringe pattern for the diametrically compressed disk. In this case both methods arrive to a quite similar result due to the high fringe order of this particular pattern. However, some differences can be observed between the results provided by the two techniques. For example, as it can be seen in Fig. 4(b), a broken fringe can be noticed at the lower left of the disk close to the border, this broken fringe is not observed in the Fig. 4(a) corresponding to the PCA algorithm.

In order to show the ability of our algorithm to deal with low order photoelastic fringe patterns, we have demodulated the isochromatics and isoclinics fringe patterns corresponding to a plastic ophthalmic lens compressed by the metallic rim of the mount. Figure 5(a) shows an example of the photoelastic fringe pattern obtained using the polariscope in plane configuration according to the scheme of Fig. 1. The combination of a bright field planar polariscope and white light illumination allows for a better demodulation of the isoclinics by avoiding the null modulation zones that appear when the retardation δ becomes an integer multiple of 2π. In order to obtain the isochromatics it is necessary to change the configuration of the polariscope from the plane polariscope to the circular one represented in Fig. 2 changing, at the same time, the illumination from white to monochromatic (Na lamp) as it can be seen in the photoelastic pattern shown in Fig. 5(b).

In Fig. 6(a) and Fig. 6(b) we have represented the wrapped phase map W[4θ] corresponding to the isoclinics obtained with PCA and PSI, respectively. It is interesting to notice that as we are dealing with a low stress distribution, only low order fringes are present in both the isoclinics and isochromatics patterns. Regarding the latter ones, we have computed the wrapped phase maps W[δ] obtained with PSI, see Fig. 6(c) and using PCA, which has been represented in Fig. 6(b). As it can be seen in these figures, the result obtained with PCA is better than the one obtained using PSI. For example, a broken fringe can be appreciated near the isotropic zone of the PSI pattern. This can be better appreciated in Figs. 6(e) and 6(f) where we have depicted the maps corresponding to the sine of the retardation angle δ. We can see that the application of PCA results in a fringe pattern with better defined fringes than the ones obtained with the 8-step PSI algorithm. For checking the goodness of the isochromatic phase demodulation, the smoothness of the sine of the recovered retardation is a good indicator, because of the odd character of the sine which highlights any possible sign recovery errors in the isochromatics, as can be observed in Fig. 6(f).

The best recovery of the wrapped phase map corresponding to the retardation W[δ] obtained when applying the PCA is due to the fact that this is an intrinsically asynchronous technique which it is not limited in the number of images that it can process, so we have taken $N=20$ images for PCA and just 8 images for PSI. Of course, it would be possible to employ other synchronous PSI algorithms with a greater number of steps but, even in this case, PCA would be also advantageous because it would not require the calibration of the rotation angles of the polariscope while PSI techniques do require this procedure.

To perform the computations whose results are shown in Figs. 3 to 6 we have used the software Code 1 [17], written in MatLab programming language.

## 4. Conclusions

We have developed a new technique for processing photoelastic fringe patterns based on the use of PCA in phase shifting interferometry. The basis of the technique is the acquisition of several images using first a plane polariscope to obtain the isoclinics and, afterwards, a circular polariscope to get the isochromatics. In both cases, one or more elements of the polariscope are rotated randomly (for isoclinics both polariscopes and for isochromatics one quarter wave plate with two different positions for the analyzer) and an arbitrary number of images (unlimited for static samples) are taken in order to get the phase. After taking the images, the isoclinics are processed first in order to obtain the phase map 2π which is needed in order to compute the isochromatics. In this later case, due to the intrinsic sign uncertainty associated with the employment of a PCA algorithm, it is necessary to compute eight synthetic fringe patterns and their respective modulations. From this set of synthetic fringe pattern the one which presents greater modulation is selected as the resulting isochromatics wrapped phase map.

PCA algorithm can be easily implemented in an automatic polariscope in a fast and efficient way, being its main advantages its asynchronous nature which prevents the need of a calibration and its good noise rejection properties. The latter is due to its ability to work with any number of images, so the algorithm can be adjusted in such a way that the greater the noise, the greater the number of images taken to deal with this noise. We have proved these qualities by performing first an experiment in which we have compared the results obtained with our algorithm to those obtained with a standard eight step PSI algorithm for a known stress distribution (in this case a diametrically compressed disk). Similar results have been obtained in other experiment in which we have shown the ability of PCA to demodulate photoelastic fringe patterns associated with low stress distributions such the one produced by a metallic frame over an ill-fitted ophthalmic lens.

## Acknowledgments

The authors wish to thank the Spanish Ministerio de Economia y Competitividad, MINECO for its economic support to this work within the framework of project DPI2012-36103.

## References and links

**1. **K. Ramesh, *Digital Photoelasticity: Advanced Techniques and Applications, Volume 1*, (Springer-Verlag, 2000). Available at: https://books.google.com/books/about/Digital_Photoelasticity.html?id=f8hRAAAAMAAJ&pgis=1 [Accessed August 12, 2015]. [CrossRef]

**2. **M. Servín, J. A. Quiroga, and M. Padilla, *Fringe Pattern Analysis for Optical Metrology: Theory, Algorithms, and Applications*, (Wiley, 2014).

**3. **K. Ramesh, T. Kasimayan, and B. Neethi Simon, “Digital photoelasticiy - a comprehensive review,” J. Strain Analysis **46**(4), 245–266 (2011). [CrossRef]

**4. **J. A. Quiroga and A. González-Cano, “Phase measuring algorithm for extraction of isochromatics of photoelastic fringe patterns,” Appl. Opt. **36**(32), 8397–8402 (1997). [CrossRef] [PubMed]

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

**6. **J. C. Estrada, M. Servin, and J. A. Quiroga, “Easy and straightforward construction of wideband phase-shifting algorithms for interferometry,” Opt. Lett. **34**(4), 413–415 (2009). [CrossRef] [PubMed]

**7. **J. C. Estrada, M. Servin, and J. A. Quiroga, “A self-tuning phase-shifting algorithm for interferometry,” Opt. Express **18**(3), 2632–2638 (2010). [CrossRef] [PubMed]

**8. **P. Carré, “Installation et utilisation du comparateur photoelectrique et interferentiel du Bureau International des Poids et Mesures,” Metrologia **2**(1), 13–23 (1966). [CrossRef]

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

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

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

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

**14. **J. A. Quiroga and A. González-Cano, “Separation of isoclinics and isochromatics from photoelastic data with a regularized phase-tracking technique,” Appl. Opt. **39**(17), 2931–2940 (2000). [CrossRef] [PubMed]

**15. **J. A. Quiroga, E. Pascual, and J. Villa-Hernández, “Robust isoclinic calculation for automatic analysis of photoelastic fringe patterns,” Proc. SPIE **7155**, 715530 (2008). [CrossRef]

**16. **A. Asundi, L. Tong, and C. G. Boay, “Dynamic phase-shifting photoelasticity,” Appl. Opt. **40**(22), 3654–3658 (2001). [CrossRef] [PubMed]

**17. **J. A. Quiroga and J. A. Gomez-Pedrero, (2016): “Photoleastic fringe pattern demodulation through PCA”, figshare (2016) [Retreived February 18, 2016], https://dx.doi.org/10.6084/m9.figshare.2353396.v1