1. Introduction

Principal component analysis (PCA) is a multivariate technique for data analysis [1,2] that allows determining whether and how data from a given dataset are correlated. From a dataset table composed of N observations of M variables, PCA computes a linear transform of these variables into a set of $L \le N$ “principal components” which describes accurately the data variation. Therefore, PCA determines how the set of N observations are correlated and reduces the data complexity by keeping the most relevant linear combination of the original variables. PCA has many applications in fields such as finance [3], neuroscience [4,5], bioinformatics [6,7], machine learning [8], and image processing [9–15].

PCA has been successfully applied to the problem of fringe pattern demodulation using the phase-shifting technique [13–18]. In this application, the data corresponds to N images (fringe patterns or, simply, interferograms), each one converted into a column vector of M elements (pixels), assuming a constant phase jump between any pair of interferograms. The application of PCA to this dataset results in two significant principal components: the sine and the cosine of the interferogram phase. Therefore, the interferogram phase can be extracted through an arctangent operation.

Phase-shifting PCA presents some advantages: Firstly, it is fast; secondly, it does not require an initial guess of the phase-steps, which may even be randomly distributed. Thirdly, it shows good noise-rejection properties even for noisy or low contrast interferograms [19,20]. Recent improvements in the application of phase-shifting PCA are the application of VU factorization to deal with varying background illumination and uneven distribution of phase jumps [21], the use of Fourier analysis for improving the algorithm performance [22], and an enhanced phase-shifting PCA based on a robust weighted PCA technique [23].

One of the drawbacks of phase-shifting algorithms (including PCA) is that noisy or low contrast interferograms, usually require processing a significant number of images to get the desired accuracy in the phase retrieval. This leads to storing a considerable amount of data in the computer memory for large images. This issue is compounded by the fact that usually, it is not possible to know, a priori, the optimum number of images necessary to determine the phase map with the required accuracy, and thus, this number may be over or underestimated. Finally, current PCA phase-shifting algorithms are not well suited for “on the fly” phase demodulation, as they require taking a large number of images before processing them.

In this paper, we propose to implement an incremental phase-shifting PCA technique to address some of the issues listed before. Used in conjunction with a parameter that indicates the quality of the retrieved phase map, the incremental algorithm proposed here can stop once the requested phase retrieval accuracy is achieved, avoiding, in this way, processing more images than necessary. Therefore, the proposed algorithm allows for “on the fly” demodulation, which may be interesting for applications such as quality control or remote sensing. We checked the technique with both synthetic and experimental interferograms to show the feasibility of our proposed approach.

The rest of the paper is organized as follows. Next, we present the theoretical foundations of the technique. Afterward, we apply the method to both synthetic and experimental interferograms showing the goodness of the algorithm. Finally, conclusions are drawn at the end of the paper.

2. Theoretical foundations

For the sake of completeness, we first summarize the phase-shifting PCA algorithm [13–14] in the next subsection. Then, we introduce the proposed incremental phase-shifting PCA approach. In the following, matrices will be represented by bold-typeface uppercase letters, either Latin or Greek, vectors by bold-typeface lowercase letters, and scalars by letters (upper- or lowercase) with the standard typeface.

2.1 Standard phase-shifting PCA algorithm

The starting point of the phase-shifting PCA algorithm is a set of N interferograms with a given phase-shift between any pair of them. Mathematically, each interferogram is described by a matrix, ${{\boldsymbol I}_{\; k = 1,2, \ldots ,N}}$, being

The first step of the standard phase-shifting PCA algorithm is grouping the data in a matrix ${\mathbf X}$ defined as

Next step, we have the background suppression. Notice that we can subtract the mean computed column-wise, or the mean computed row-wise. The former one is equivalent to a spatial average along all the pixels of the image, while the latter one is equivalent to a temporal average [18]. In this work, we will use the temporal average as it is more effective for noise suppression so that we will rewrite matrix ${\mathbf X}$ as

The goal of the PCA algorithm is to determine a linear transformation that converts the correlation matrix ${\mathbf C} = {{\mathbf X}^\textrm{T}}{\mathbf X}$ (where $\textrm{T}$ stands for the transpose operator) into a non-zero positive-defined diagonal matrix, ${\mathbf \Lambda }$, with eigenvectors $\{{{\sigma_1},{\sigma_2}, \ldots {\sigma_N}} \}$ ordered in increasing order, so that, ${\sigma _1} > {\sigma _2} \ldots > {\sigma _N}$. In general, the first $L < N$ eigenvectors of matrix ${\mathbf \Lambda }$, are enough to describe the data variation accurately.

To compute ${\mathbf \Lambda }$, we will start with the singular value decomposition (SVD) of matrix ${\mathbf X}$

When PCA is applied to phase-shifted interferograms, it can be demonstrated [13,14] that there are only two statistically significant principal components so, in this case, $L = 2$, and $\hat{{\mathbf F}} = [{{{\mathbf f}_1},\; {{\mathbf f}_2}} ]$ where ${{\mathbf f}_{1,2}}$ are column vectors with M elements. Moreover, as it is shown in Refs. [13,14]., after arranging these first principal components as matrices, ${{\mathbf F}_{1,2}}$, with the same dimension of the interferograms, we have that [13,14]

2.2 Incremental phase-shifting PCA algorithm

The standard phase-shifting PCA algorithm presented above has some limitations. Firstly, to accurately demodulate noisy or low contrast fringe patterns, a relatively high number of images ($N \approx 20$) is required. According to Levy and Lindenbaum [24], the complexity of the standard PCA algorithm is $O({M{N^2}} )$ and it requires the storage of $O({MN} )$ memory units. With current image sizes in mind ($M \approx 5\cdot{10^6}\; \textrm{pixels}$), the standard phase-shifting PCA should process a large number of data ($MN \approx {10^8}$ bytes). In some cases, this may affect the performance of the algorithm, being even possible to get a computer memory overload. However, the capabilities of modern computers overcome this disadvantage to some extent. Secondly, it is impossible to know, a priori, the exact number of images required to demodulate some fringe pattern. To overcome these problems, we propose an incremental PCA algorithm coupled with a quantitative parameter to determine the goodness of the phase demodulation process.

In the literature, some incremental PCA algorithms have been proposed in the past [24–26]. The workflow of these algorithms is similar: they start with a matrix, ${\mathbf Y}$, formed by the first $R$ columns of matrix ${\mathbf X}$, Obviously, $R$ should be equal or greater than the number of statistically relevant principal components, so $R \ge L$. After applying PCA to matrix ${\mathbf Y}$ we obtain a first approximation to the principal components. This approximation is refined by sequentially processing the remaining columns of ${\mathbf X}$. Notice that not all of these incremental algorithms are suited for processing large amount of data, such as images. For example, the algorithm of Lippi et al. [25] requires the computation of $M \times M$ matrices. For a typical high resolution image of a 5 megapixel camera, this would computation of $2.5\cdot{10^{13}}$ bytes matrices, in uint8 precision, which is clearly not feasible.

In this work, we have selected the algorithm of Levy and Lindenbaum [24], which is suited for large images. However, there is a fundamental change in our implementation of Levy’s algorithm in how we perform background subtraction. In the work of Levy and Lindenbaum background suppression is performed individually for each image, this is equivalent to subtracting the spatial average, as commented before. However, we have found that, for the problem of fringe pattern processing, it is more effective, for noise suppression, to subtract the temporal average. In the standard PCA algorithm, given the data matrix ${\mathbf X} = [{{{\mathbf i}_1},{{\mathbf i}_2}, \ldots ,{{\mathbf i}_N}} ]$, spatial averaging is equivalent to subtracting the mean computed column-wise while temporal averaging implies subtracting the mean row-wise. Therefore, for the standard PCA algorithm, both spatial and temporal averaging can be implemented indistinctively. However, for the temporal PCA algorithm, to compute the temporal average we need to update sequentially the mean vector, so we have modified Levy and Lindenaum algorithm accordingly. Therefore, the algorithm proposed can be summarized in the following steps.

2.3 Determining the minimum number of interferograms

One of the advantages of the incremental phase-shifting PCA algorithm is that it demodulates the phase using the exact number of interferograms needed to get the desired accuracy. For this, it is necessary to use some quality parameter that determines the goodness of the phase demodulation. In Ref. [19], the following phase error indicator is proposed

In the next section, we will show the results obtained when we applied our algorithm to both simulated and experimental interferograms.

3. Experimental results

3.1 Numerical simulations

To check the feasibility of the proposed algorithm, we have carried out first a numerical simulation. We have generated a set of N = 20 synthetic fringe patterns with low modulation ($m = 0.15$) and high uniform noise ($SNR\; = \; 25$) as shown in Fig. 1. Notice that, in this work, we define the signal to noise ratio as $SNR = {b^2}/\sigma _\eta ^2$, where ${\sigma _\eta }$ is the noise standard deviation. The size of the fringe patterns is $512 \times 512$ pixels with 256 gray levels, so $M = {512^2} = 262144$ pixels. Each fringe pattern presented a random global phase shift, comprised within the interval $[{ - \pi ,\pi } ]$ rad. We have demodulated this set of fringe patterns using the standard PCA demodulation algorithm and the proposed incremental technique.

 

Fig. 1. First, a), and last, b), interferogram of the sequence employed to test the proposed incremental PCA algorithm. Both interferograms present low contrast, $m\; = \; 0.15$, and high uniform noise, $SNR\; = \; 5$, with a random phase shift between them.

Download Full Size | PDF

We will first present the results obtained with the standard PCA demodulation algorithm in Fig. 2. We have depicted in 2a) and 2b) the first and second principal components, ${{\boldsymbol F}_1}$ and ${{\boldsymbol F}_2}$, respectively, and the phase obtained after applying Eq. (9) in Fig. 2(c)). The actual phase (wrapped) is represented in Fig. 2(d)). As we can see, the PCA algorithm can demodulate properly even low contrast, noisy fringe patterns as the ones shown in Fig. 1.

 

Fig. 2. Results obtained after applying the standard PCA algorithm a) first and b) second principal component, c) demodulated wrapped phase and d) original wrapped phase. Notice that the algorithm has been able to demodulate correctly the phase. Still, the resulting phase map presents noise due to the noisy input fringe patterns, like those shown in Fig. 1.

Download Full Size | PDF

Next, we show the results of applying the proposed incremental algorithm. In Fig. 3, we present the algorithm's evolution by depicting the demodulated phase after processing 5, 10, 15, and, finally, 20 fringe patterns of the set using the incremental PCA algorithm. As expected, the result, after processing the whole set of fringe patterns with the incremental algorithm (Fig. 3(d)), is similar to that obtained using the standard PCA algorithm (Fig. 2(c)). However, the incremental PCA only needs to store at least two fringe patterns in the computer memory, while the standard algorithm should store the whole set.

 

Fig. 3. Phase demodulation with the incremental PCA algorithm after processing a) 5, b) 10, c) 15, and d) 20 fringe patterns of the set. Notice that the demodulated phase after processing the whole set of fringe patterns is similar to the phase depicted in Fig. 2(d)).

Download Full Size | PDF

Regarding performance, the incremental PCA carries out the demodulation of the entire set of $512 \times 512$ fringe patterns in 0.31 seconds, while the standard PCA algorithm employs a similar time of 0.32 seconds. In both cases, we use a computer with an Intel Core i7-8700 CPU with 32 Gb of RAM installed and a graphics card NVidia GeForce GTX 1050 with 1Gb of memory, running Matlab 2020a in Windows 10. The code employed can be found in Code 1, Ref. [27].

We have also compared our method to a robust PCA-based phase demodulation technique [28], especially suited for precise phase retrieval under harsh conditions. In Fig. 4 we show the results obtained after processing the whole set of 20 interferograms. We can see that the roust PCA method retrieves the phase with less error than incremental PCA but, as this method requires to compute an extended set of interferograms, the computation time rises to 38.8 seconds, a hundred-fold increment.

 

Fig. 4. a)Phase demodulated with the incremental PCA algorithm after processing 20 fringe patterns. b) Phase demodulated using the robust PCA-based algorithm described in Ref. [28].

Download Full Size | PDF

Besides the advantage of not storing at the same time the whole set of fringe patterns in the computer memory, the incremental PCA algorithm has the advantage of “on the fly” demodulation. This means, among other things, that we can stop the algorithm before processing the whole set of fringe patterns once we reach an acceptable demodulation result. For example, as we can see in Fig. 3, after processing 15 interferograms, the visual aspect of the demodulated phase (Fig. 3(c)) is quite similar to that after processing the whole set of 20 interferograms (Fig. 3(d)). Therefore, we need to define some parameter that gives us a quantitative indicator of the quality of the phase demodulation. In this way, we can process a set of fringe patterns until the quality parameter reaches a fixed value, avoiding, in this way, processing more fringe patterns than the necessary minimum to achieve the desired result.

As said before, in this work we will use the quality parameter proposed by Vargas et al. [19], and defined in Eq. (10). To show the evolution of this parameter, we have carried out another simulation in which we have demodulated three sets of synthetic fringe patterns corresponding to the same phase distribution depicted in Fig. 2(d)), but with different values of contrast and SNR. The first set is formed by noisy (SNR = 25) and low contrast (c = 0.15) fringe patterns, the same case of Figs. 1, 2, and 3. The second set consists of fringe patterns with medium contrast, c = 0.3, and a moderate SNR of 50. Finally, the last pattern set is formed by high-contrast (c = 0.6) fringes with negligible noise (SNR = 75). Each set contains 40 images, and we have demodulated them using the proposed incremental PCA algorithm. The demodulation results are shown in Fig. 5

 

Fig. 5. Demodulation results of three different sets of fringe patterns using the proposed incremental PCA algorithm. a) The first set formed by low-contrast and noisy fringe patterns, b) the second set formed by average-contrast and average-noise fringe patterns, and c) the third set formed by high-contrast and low noise fringe patterns. Notice that a global phase jump may exist between the recovered wrapped phase maps.

Download Full Size | PDF

As shown in Fig. 5, the algorithm successfully demodulated the three sets, obtaining a cleaner phase map for the high-contrast, low noise set, as expected. We have also computed the fringe pattern modulation computed as $m = \sqrt {{{\mathbf F}_1} + {{\mathbf F}_2}} $ where ${{\mathbf F}_1}$ and ${{\mathbf F}_2}$ are the significant principal components used to compute the phase in Eq. (9). The resulting modulation maps are plot in Fig. 6. Although these maps should show a uniform modulation, a double-frequency distortion appears. This distortion is linked to the fact that the sum $\mathop \sum \nolimits_n^N \cos {\delta _n}\sin {\delta _n}$ is not exactly zero, as assumed by the general phase-shifting PCA algorithm [14], which produces an error in the modulation that goes as $\sin ({2{\mathbf \Phi }} )$, hence the double-frequency distortion. The recovered phase also present this structure as we can see in Fig. 7, where we have represented the differences between the computed wrapped phase and the theoretical one used to generate the fringe patterns (also wrapped). Despite this structure, the RMS error diminishes as the SNR increases, as expected from the plots of Fig. 5, and the value of this magnitude is around 0.003 rads.

 

Fig. 6. Modulation of the three different sets of fringe patterns demodulated with the proposed incremental PCA algorithm for the three levels of SNR and contrast defined in Fig. 5. Notice that a double-frequency ripple distortion that appears due asynchronous nature of the phase-shifting PCA method

Download Full Size | PDF

 

Fig. 7. Difference between the wrapped phase demodulated with the proposed incremental PCA algorithm and the exact wrapped phase for the three levels of noise and contrast described in Fig. 5. Notice that a double-frequency ripple distortion that appears due asynchronous nature of the phase-shifting PCA method

Download Full Size | PDF

We have plotted in Fig. 8 the evolution of the phase-error indicator defined in Eq. (10), expressed as a percentage, against the number of iterations. We can see that, in general, the parameter tends to decrease with the number of images. However, while the evolution of the quality parameter is smooth for noisy fringe patterns, when the noise diminishes, the quality parameter presents oscillations, particularly for the low-noise, high-contrast interferograms depicted as green diamonds in Fig. 8. This behavior indicates that processing more images than necessary for such low-noise fringe patterns does not improve the quality of the phase recovered. For example, when SNR = 75 and c =0.6, the quality parameter reaches a minimum after processing ten images. We have checked that the phase map recovered after processing ten images is the same as the one depicted in Fig. 4(c)), obtained after processing 40 images. Therefore, a good strategy for implementing the quality parameter could be processing fringe patterns until the quality parameter reaches a local minimum value, avoiding, in this way, processing more images than necessary.

 

Fig. 8. Evolution of the phase-error indicator against the number of iterations for three sets of synthetic fringe patterns with the same phase but different contrast and noise. The triangles correspond to low-contrast (c = 0.35), high noise (SNR = 40), the squares to medium-contrast (c = 0.5) average noise (SNR = 80) and the diamonds to high-contrast (c = 0.7), low-noise (SNR = 240) patterns. Notice the decrement of the phase-error indicator value with the increasing number of iterations.

Download Full Size | PDF

3.2 Experimental results

Next, we will show the results obtained when we have applied our algorithm to demodulate experimental fringe patterns. In Fig. 9, we present the first and last image of a set of N = 19 interferograms whose size is $600 \times 800$ pixels, thus $M = 4.8\cdot{10^5}$ pixels. There is a fixed phase shift, not known in advance, between each pair of interferograms of this sequence. We can see that these interferograms present medium contrast with moderate noise. For this example, the interferograms are masked. However, this is not a problem for either the incremental or the standard PCA demodulation algorithms.

 

Fig. 9. First, a) and last b) fringe pattern of the sequence corresponding to an interferometric experiment. There is a global phase jump between these interferograms, otherwise identical.

Download Full Size | PDF

We have demodulated this set of interferograms with both the standard and incremental PCA algorithms. We have computed in both cases the evolution of the quality parameter, $\gamma $, or phase-error indicator, against the number of processed interferograms. In Fig. 10, we present this evolution plot showing that, in both cases, the quality parameter evolves similarly.

 

Fig. 10. Evolution of the quality parameter as a function of the number of images processed for the set of experimental interferograms depicted in Fig. 7. We can see a similar evolution for the incremental PCA (blue squares) and the standard (red triangles) algorithms.

Download Full Size | PDF

In Fig. 11(a)) we show the result after demodulating this fringe pattern set with the proposed incremental PCA algorithm, while the result obtained after applying the standard PCA algorithm is depicted in Fig. 11(b)). The result of Fig. 11(a)) has been obtained after processing 14 images, for which the Phase-error indicator of Fig. 10) reached a value lower than the threshold ${\gamma _0} = 10\; \% $. At the same time, the phase map depicted in Fig. 11(b)) has been obtained after processing the whole set of $N = 17$ fringe patterns. We can see a good agreement between the two methods, but it is noticeable that, with the proposed algorithm, we have to process only a fraction of the whole dataset. Notice that, to get the phase maps of Fig. 11, we have corrected a global phase jump between them.

 

Fig. 11. Wrapped phase obtained after demodulating the fringe pattern whose first and last image are presented in Fig. 9(a)) and b), respectively, using the proposed incremental algorithm. The result of the incremental PCA demodulation algorithm after processing 14 fringe patterns is shown in Fig. 11(a)) and the result of the standard PCA algorithm after processing the whole dataset of $N = 17$ fringe patterns is displayed in Fig. 11(b)). Notice the agreement between both phase maps.

Download Full Size | PDF

In Fig. 12, we show another pair of fringe patterns corresponding also with the first and last of a sequence of $N = 9$ fringe patterns obtained in an interferometry experiment. In this example, the size of the interferograms is $481 \times 641,$ so $M = 308321$ pixels, and there exists an unknown phase step between each interferogram.

 

Fig. 12. a) First and b) last fringe pattern of a sequence of interferograms. Compared to the fringe patterns depicted in Fig. 9, they present higher contrast and more noise, including a noticeable periodic noise in the lower left part of the images.

Download Full Size | PDF

The results of the demodulation of the fringe pattern sequence depicted in Fig. 12 are shown in Fig. 13. The fringe pattern showed in Fig. 13(a)) has been obtained after processing six fringe patterns, which was the point at which the phase error indicator reaches the threshold value ${\gamma _0} = 10\; \% $. Again, for comparison purposes, we present the result after demodulating the whole sequence of $N = 9$ fringe patterns with the standard PCA algorithm in Fig. 13(b)). A constant phase jump has been removed as we did before.

 

Fig. 13. Map corresponding to the wrapped phase obtained after demodulating the fringe pattern sequence corresponding to Fig. 12, with a) the proposed incremental PCA after processing six images and b) the standard PCA demodulation algorithm after processing the whole dataset.

Download Full Size | PDF

Finally, in Fig. 14, we show the evolution of the quality parameter. In this case, we can see a similar evolution for both the standard and incremental PCA algorithms again, with a local minimum of the quality parameter for six interferograms processed, suggesting that processing more than six for this set of interferograms fringe patterns could be redundant, as shown in Fig. 13. Both datasets and the code necessary for processing can be downloaded from Code 1, Ref. [27].

 

Fig. 14. Evolution of the quality parameter for the set of fringe patterns depicted in Fig. 9. Notice the slight increment of the quality parameter after processing four interferograms for both the standard and incremental PCA algorithms.

Download Full Size | PDF

4. Conclusions

We have devised an incremental principal component analysis algorithm for the fringe pattern demodulation problem adapting an existing incremental SDV algorithm. As a result, it is now possible to perform “on-the-fly” processing of fringe patterns without the need for storing a significant number of images in the computer memory. Furthermore, by using an error-phase quality parameter, it is possible to stop the algorithm once the desired accuracy is achieved, improving, in this way, the performance of the PCA algorithm as only the absolute minimum number of images are processed.

Numerical simulations have been performed to check the algorithm for different values of contrast and noise. The algorithm presents good noise rejection capabilities, being able to deal with low contrast and noisy fringe patterns with similar results to the standard PCA algorithm.

Finally, the results obtained for experimental fringe patterns show that the proposed incremental algorithm is feasible and gives similar results to the standard PCA algorithm. For these experimental fringe patterns, the quality parameter can be used to reduce the number of fringe patterns to be processed, incrementing, in this way, the efficiency of the standard PCA phase demodulation algorithm.