## Abstract

Phase extraction methods based on the principal component analysis (PCA) can extract objective phase from phase-shifted fringes without any prior knowledge about their shift steps. Although it is fast and easy to implement, many fringe images are needed for extracting the phase accurately from noisy fringes. In this paper, a simple extension of the PCA method for reducing extraction error is proposed. It can effectively reduce influence from random noise, while most of the advantages of the PCA method is inherited because it only modifies the construction process of the data matrix from fringes. Although it takes more time because size of the data matrix to be decomposed is larger, computational time of the proposed method is shown to be reasonably fast by using the iterative singular value decomposition algorithm. Numerical experiments confirmed that the proposed method can reduce extraction error even when the number of interferograms is small.

© 2016 Optical Society of America

## 1. Introduction

Phase-shifting interferometry (PSI) is the well-known and widely used technique for measuring objective phase from interferograms [1–6]. In PSI, *N* fringe images are acquired by shifting the phase, and the objective phase modulating the fringes is extracted from the images. For the extraction, some traditional algorithms require information of the shifting steps. However, it is not easy to obtain precise information about them because of practical implementation error of measuring instruments. Therefore, there has been many attempts to estimate the unknown phase-shifting steps and/or objective phase from observed data [7–18].

A relatively new approach for the phase estimation is based on the principal component analysis (PCA) [19–24]. It treats each pixel of *N* fringe images as an element of an *N*-dimensional vector space. Then, the data vectors corresponding to all pixels within the image are lying on a specific linear subspace of the vector space because of the sinusoidal nature of fringes. PCA is a data analysis method that can automatically identify a basis of such linear subspace. Therefore, PCA can decompose the cosine and sine components, which is necessary for the phase extraction, from the data vectors. The PCA method has several advantages including its low computational cost, non-iterative nature which does not require any initial guess, easiness of treating spatially varying background illumination and contrast, and easiness of implementation. Although the PCA method has such attractive properties, many fringe images are required in order to obtain an accurate result. This requirement prevents its applicability to some situations where obtaining a lot of interferograms is difficult. Therefore, it is better to consider an extension of the PCA method so that more accuracy can be attained with less interferograms.

In this paper, a simple extension of the PCA method is proposed. The basic idea of the proposed method is to integrate spatial information into the analysis process so that phase extraction error related to random noise is reduced. It is quite simple because the only modification is in the construction of a data matrix from fringe images. Therefore, the proposed method inherits most of the advantages of the PCA method. Although its computational cost becomes higher, use of the iterative singular value decomposition (SVD) algorithm allows reasonably fast computation. The most important aspect of the proposed method is that combining it to existing PCA-based algorithms is easy. That is, the extension is independent from the experimental set-up and applications. Thus, it has possibility of improving every PCA-based phase extraction algorithm.

This paper is organized as follows. In Section 2, the PCA method is explained in terms of low-rank matrix decomposition using SVD. Then, the proposed method is introduced in Section 3, where its connection with spatial derivative is also discussed in the same section. Results for numerical experiments are shown in Section 4, and the paper is concluded in Section 5.

## 2. PCA-based phase extraction

Let *I*^{[n]}(*x*) be an *n*-th image of phase-shifted fringes:

*x*∈ ℝ

^{2}denotes the position vector,

*ϕ*is the objective phase,

*B*is background illumination,

*A*is amplitude of the fringe, and

*δ*

^{[n]}is the unknown phase shift. Its discrete sample at (

*i*,

*j*)-th pixel will be denoted by ${I}_{i,j}^{[n]}$. The aim of phase extraction is to reconstruct

*ϕ*, which contains the desired information, from

*N*shifted fringe images.

#### 2.1. Data matrix constructed from fringes

In PCA-based phase extraction, shifted fringes at each pixel are concatenated to form an *N*-dimensional column vector,

*x*

_{1},

*x*

_{2}, …,

*x*] is the horizontal concatenation of {

_{N}*x*}, and

_{n}*x*is the transpose of

^{T}*x*. This vector can be regarded as a data vector obtained from a single pixel. Then, a data matrix is constructed by concatenating

*d*for every pixels horizontally,

_{i,j}*N*×

*N*

_{v}

*N*

_{h}, where

*N*

_{v}and

*N*

_{h}are respectively numbers of vertical and horizontal pixels. Note that one may construct

*D*in a different manner because the order of concatenation is not important here.

Since Eq. (1) can be rewritten as

*d*is sum of three components: where

_{i,j}*o*= [1, 1, …, 1]

*is the*

^{T}*N*-dimensional vector whose elements are all one,

*B*,

_{i,j}*C*=

_{i,j}*A*cos(

_{i,j}*ϕ*

*) and*

_{i,j}*S*= −

_{i,j}*A*sin(

_{i,j}*ϕ*

*) are scalers independent from the phase shift, and*

_{i,j}*c*and

*s*are vectors defined only by

*δ*

^{[}

^{n}^{]}as

Therefore, rank of the matrix *D* is three, i.e., *D* is sum of the three rank-one matrices defined by horizontal concatenation of each column vector, *B _{i,j} o*,

*C*, and

_{i,j}c*S*, as in Eq. (3).

_{i,j}s#### 2.2. Phase extraction by decomposing data matrix into rank-one matrices

SVD is a method of decomposing a matrix *D* into three matrices as

*U*and

*V*are orthonormal matrices whose sizes are respectively

*N × N*and

*N × N*

_{v}

*N*

_{h}, and Σ is a diagonal matrix consisting of singular values as Σ = diag(

*σ*

_{1},

*σ*

_{2}, …,

*σ*) with

_{N}*σ*

_{1}≥

*σ*

_{2}≥ ⋯

*σ*. This decomposition can be interpreted as where

_{N}*u*and

_{n}*v*are respectively

_{n}*n*-th column vectors of

*U*and

*V*. That is, a matrix can be decomposed into rank-one matrices $\{{u}_{n}{v}_{n}^{T}\}$ through SVD.

As the data matrix *D* in Eq. (3) consists of the three rank-one matrices, SVD is applied to decompose *D* into them. If the decomposition is correct, *C _{i,j}* and

*S*can be separated from the matrices. Then, the objective phase is extracted by

_{i,j}*z*] denotes the principal value of complex argument of

*z*.

In the above method, accuracy of the decomposition directly affect the accuracy of phase extraction. To improve the accuracy, mean values of each column vector of *D* is usually subtracted in advance. Then, the first two right singular vectors, *v*_{1} and *v*_{2} corresponding to the first and second singular values, contain *C _{i,j}* and

*S*as their elements. Thus,

_{i,j}*ϕ*

*can be reconstructed by using SVD up to its global sign.*

_{i,j}We should note that most of the articles about PCA-based methods adopted the eigenvalue decomposition for explanation instead of SVD. Since

*U* obtained by SVD of *D* coincides with the one obtained by the eigenvalue decomposition of *DD*^{T} (by assuming that there is no eigenvalue degeneracy). Therefore, the above SVD-based explanation is consistent with the original PCA-based method. Moreover, most of the articles states that the background illumination term *B _{i,j} o* in Eq. (5) must be removed before the decomposition. For many cases, PCA-based methods still works without subtracting it by selecting appropriate vectors corresponding to the cosine and sine terms. However, it is definitely better to remove it in advance because the rank-one matrices likely be mixed up in the resultant decomposition that causes extraction error.

## 3. Reducing random noise related error

In this section, we propose a simple modification of the above PCA-based phase extraction method to improve its accuracy. The basic idea is to consider the spatial relation of adjacent pixels within the decomposition.

Since fringes have specific structure based on their sinusoidal nature, adjacent pixels of a fringe image also have specific properties. To utilize such relation automatically, the data matrix of the PCA method is extended here. The extension allows better separation of random noise because noise does not follow the specific structure of fringes.

In the following, spatial derivative of a fringe image is considered first because it can be regarded as relationship of the neighborhood. Then, an extended data matrix is introduced, whose property is discussed in connection with the derivative. Finally, the proposed method is explained as a PCA method applied to the extended data matrix, where the ordinary PCA method is included as a special case.

#### 3.1. Data matrix constructed from spatial derivative of fringes

Since fringes are generated through the cosine function, their gradient should also have similar properties to them. Derivative of a fringe [Eq. (1)] is

*f*′ denotes spatial derivative of

*f*. Spatial dependency was omitted for simplicity. Then, a data vector,

Based on the above observation, it may be possible to improve accuracy of estimation of *c* and *s* by including *d*′, which can be approximated by difference between adjacent pixels, into the ordinary data matrix *D* in Eq. (3). However, the present paper will not consider such gradient-based approach because estimating derivative from real data tends to decrease the signal-to-noise ratio (SNR) that can be another source of error. What we want to emphasize here is that a data matrix constructed from relation between adjacent pixels, or derivative, also possesses the low-rank structure.

#### 3.2. Extending data matrix by combining adjacent pixels

Here, a data matrix constructed from extended data vectors is considered. Let
$\tilde{d}$ be a 2*N*-dimensional column vector defined by

*N*×

*N*

_{v}

*N*

_{h}data matrix,

*i*,

*j*)-th pixel but also information of its adjacent. Therefore, it should be natural to think about decomposing $\tilde{D}$ so that spatial information of data is taken into account.

The extended data vector can be interpreted as

*0*= [0, 0, …, 0]

*is the*

^{T}*N*-dimensional vector whose elements are all zero. This equation indicates that $\tilde{D}$ is a low-rank matrix because each vector in the right hand side of Eq. (17) ends up with low-rank data matrices as discussed in the previous sections. Moreover, singular values of $\tilde{D}$ should concentrate in the few leading singular values because

*d*and ${{d}^{\prime}}_{i,j}$ share similar properties as in Eqs. (5) and (14).

_{i,j}On the other hand, random noise does not have such adjacent relationship that holds at every pixel within the fringe images. Therefore, components related to noise should spread over all singular vectors after the decomposition. Since only the first two singular vectors are needed for phase extraction, their SNR increases as number of singular vectors increases. That is, the energy of noise in each singular value decreases because of its spreading over other singular values, while the energy of fringe components concentrates in the first two singular values. Thus, extracted phase has less random-noise-related error as more data vectors are concatenated into the extended data matrix.

#### 3.3. Proposed method

Based upon the above discussion, we propose an extension of the PCA-based phase extraction method by combining an extended data matrix. The proposed method is defined as the following procedure:

- Import
*N*fringe images $\{{I}_{i,j}^{[1]},{I}_{i,j}^{[2]},\dots ,{I}_{i,j}^{[N]}\}$. - Select
*M*adjacent pixels centered at (*i*,*j*)-th pixel, including (*i*,*j*). - Subtract mean values from each ${\tilde{d}}_{i,j}$ to obtain zero-mean extended data vectors ${\stackrel{\u2323}{d}}_{i,j}={\tilde{d}}_{i,j}-{\mu}_{i,j}$, where
*μ*is the mean of ${\tilde{d}}_{i,j}$ calculated from its_{i,j}*MN*elements. - Construct an
*MN*×*N*_{v}*N*_{h}data matrix*Ď*by concatenating the zero-mean extended data vectors $\{{\tilde{d}}_{i,j}\}$ horizontally as in Eq. (16). - Perform SVD to decompose
*Ď*into rank-one matrices:*Ď*=*U*Σ*V*^{T}. - Extract objective wrapped phase from the first two right singular vectors of
*Ď*as $\widehat{\varphi}=\text{Arg}[{v}_{1}-\iota {v}_{2}]$.

In the last step, instead of using *V* directly, one may project data by Σ^{−1}*U ^{T}* to obtain

*V*as in the original PCA method. This way of obtaining

*V*is useful when the images contain an untrusted region where no fringe pattern is available. In that case,

*U*and Σ are decomposed from data only in a trusted region, and then whole images including the untrusted region are obtained by multiplying Σ

^{−1}

*U*from left.

^{T}As discussed in the previous subsection, the proposed method should eliminate more noise when more adjacent pixels are integrated in the extended data matrix. On the other hand, the proposed method requires more computational resources because *Ď* contains *M* times more elements than the ordinary data matrix, and the typical computational cost of SVD is *O*(*M*^{2}*N*^{2}*N*_{v}*N*_{h}). Fortunately, only the first two singular vectors of *Ď* are necessary for extracting the phase. Therefore, an iterative algorithm can be applied to compute only the necessary vectors with less computational time comparing to the full SVD. Note that such iterative SVD algorithm should be easily found in common linear algebra libraries (for example in MATLAB, it is implemented as the
`svds` command by default).

In the proposed method, the number of adjacent pixels *M* is the parameter balancing degree of noise reduction and computational time. In addition, geometry of the adjacent pixels (choice of the *M* pixels which are considered as adjacent) also has impact on the extraction result. For concrete examples with different *M*, five
${\tilde{d}}_{i,j}$ used in the numerical experiments in the next section are shown here:

**Case (b):***M*= 9,$$\begin{array}{l}{\tilde{d}}_{i,j}=[{d}_{i,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-1}^{T},\\ \phantom{\rule{1.8em}{0ex}}{d}_{i+1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j-1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j-1}^{T}{]}^{T},\end{array}$$**Case (c):***M*= 13,$$\begin{array}{l}{\tilde{d}}_{i,j}=[{d}_{i,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j-1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+2,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-2}^{T}{]}^{T},\end{array}$$**Case (d):***M*= 21,$$\begin{array}{l}{\tilde{d}}_{i,j}=[{d}_{i,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j-1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+2,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-2}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+2,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+2,j-1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+1,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j-2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j-2}^{T}{]}^{T},\end{array}$$**Case (e):***M*= 25,$$\begin{array}{l}{\tilde{d}}_{i,j}=[{d}_{i,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j-1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+2,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i,j-2}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+2,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+2,j-1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j+1}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j-1}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+1,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-1,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+1,j-2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j-2}^{T},\\ \phantom{\rule{2em}{0ex}}{d}_{i+2,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i+2,j-2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j+2}^{T},\phantom{\rule{0.2em}{0ex}}{d}_{i-2,j-2}^{T}{]}^{T},\end{array}$$which are depicted in Fig. 1.

Note that boundary of the fringe images must be treated properly so that *M* adjacent pixels are reasonably defined at the boundary (for example,
${d}_{{N}_{\mathrm{v}}+1,{N}_{\mathrm{h}}}$ is outside the image, which cannot be defined in a trivial manner). Since the performance of the proposed method depends on the choice of a boundary condition including periodic, symmetric, replicative, and constant boundary conditions, choosing a better condition is necessary to obtain a good result. However, the best boundary condition should depend on the application. As this paper is not particular about a specific application, we will simply ignore the boundary in the numerical experiments in the next section for avoiding this ambiguity. That is, we treat an *N*_{v} *× N*_{h} image as the (*N*_{v} − 4) × (*N*_{h} − 4) image so that adjacent pixels are available for all pixels within that region.

## 4. Simulation

In this section, the effectiveness of the proposed method is examined by two numerical simulations: fringes obtained from equispaced and random shifting steps. The reason of considering equispaced shifting steps here is that the proposed method can be compared with the standard *N*-step phase-shifting algorithm which is only applicable to the fringes with equispaced shifts.

#### 4.1. Fringes obtained from equispaced shifting steps

For the first examination, the proposed method was applied to *N* fringes obtained from known equispaced phase shifts:

*n*∈ {1, 2, …,

*N*},

*ε*is the Gaussian random noise whose standard deviation is one, and

*σ*≥ 0 is the standard deviation. This simple setting was chosen because a reference for measuring the performance can be obtained by the standard

*N*-step phase-shifting algorithm [1],

The proposed method was compared with this reference to confirm its effectiveness. For the parameters, *B* = 2, *A* = 1,
$\varphi =4\pi ({x}_{1}^{2}+{x}_{2}^{2})$ was chosen, where *x*_{1}, *x*_{2} ∈ [−1, 1], and the image size was 300 × 300. The standard deviation of noise *σ* was chosen from {0.1, 0.2, 0.3, 0.4, 0.5}.

Figure 2 shows an example of obtained result for *N* = 3 and *σ* = 0.5, where the first three images in the first row represent a fringe image and phase maps, and other images illustrate error maps. Error maps correspond to the phase-shifting algorithm in Eq. (24), the ordinary PCA method described in Section 2, and the proposed method. There are five results for the proposed method corresponding to Eqs. (18)–(22) which are depicted in Fig. 1. Remind that the ordinary PCA method is a special case of the proposed method:
${\tilde{d}}_{i,j}={d}_{i,j}(M=1)$.

From the figure, it is seen that the *N*-step shifting algorithm and the PCA method obtained the similar results. This is consistent with the analysis of the ordinary PCA method in [20] which introduced an interesting link between the two methods. On the other hand, every result obtained by the proposed method has less error than them. Therefore, its ability of noise separation mentioned in the previous section is confirmed. The root-mean-square error (RMSE) of the obtained phase is shown in Fig. 3, which is calculated from unwrapped and mean subtracted results. It also confirmed the above observations.

Computational time of the proposed method, which is implemented by MATLAB using the
`svds` command, is shown in Table 1. It was measured by using five fringes (*N* = 5), whose image size was 300 × 300, with an Intel core i7 3.0 GHz processor. The construction of each extended data matrix was also included. Since *M* = 1 corresponds to the ordinary PCA method, *M* > 1 results in a larger extended data matrix, for example, the matrix for *M* = 25 consists of 25 times more elements than the ordinary PCA method. Therefore, the computational cost is obviously higher for larger *M*. In contrast, the extracted results were better for larger *M* as in Figs. 2 and 3. Thus, *M* should be adjusted to balance the performance and computational time. Although the proposed method needs more time for computation than the ordinary PCA method, it should be fast enough in practice.

#### 4.2. Fringes obtained from random shifting steps

The proposed method was also applied to randomly phase-shifted fringes:

where {*δ*

^{[}

^{n}^{]}} was obtained randomly from the uniform distribution within [0, 2

*π*]. Since the shift steps are not aligned, the

*N*-step phase-shifting algorithm cannot be applied any more. Therefore, the standard self-calibration method, the advanced iterative algorithm (AIA) [10], was used as a reference instead.

*A*,

*B*and

*ϕ*were the same as in the previous subsection.

An example of extracted results for *N* = 3 and *σ* = 0.5 is shown in Fig. 4, where {*δ*^{[n]}} = {2.3592, 2.9314, 4.3212} was obtained randomly. Although the noise level *σ* was same as in the previous example in Fig. 2, the extraction error was higher than that owing to ambiguity on the phase shift.

Since the fringe images were generated by a totally random manner, mean values of RMSE were calculated from 100 trials. However, some extreme trials had huge error comparing to the average (typically 10–100 times greater than the average) that make interpretation of the results quite difficult. Therefore, outliers whose error values were greater than a standard deviation from the mean values were excluded. The resultant mean RMSE is shown in Fig. 5. It confirms that the proposed method also works well for randomly shifted fringes.

We should note that, the proposed method with *M* ∈ {13, 21, 25} was less stable than that with *M* ∈ {5, 9} when the number of fringe images *N* was small (say, *N* ≤ 5). Therefore, we suggest that the case (b) in Eq. (19) (*M* = 9) seems a good choice in terms of performance, computational time, and stability.

## 5. Conclusions

For improving phase extraction from randomly phase-shifted fringes, a simple modification of the PCA method was proposed. The data matrix was extended to contain adjacent pixels so that spatial relation of a fringe is taken into account. Some discussions on the extended data matrix in connection with spatial derivative were also presented in order to show the basic principle of the proposed method.

The numerical simulations indicated that the proposed method can extract objective phase with less error than the PCA method. Therefore, by using the proposed method, it is possible to obtain a similar result from less interferograms, or a more accurate result from the same number of interferograms. Although the computational time becomes longer by the modification, it can be calculated reasonably fast thanks to the iterative SVD algorithm. Since the proposed method is a quite simple modification, it should be easy to combine with the existing PCA-based methods. Thus, the proposed method can be considered as a promising extension of the PCA method which should be effective for many applications of PSI. One possible application of the proposed method is the single-shot PSI whose implementation error on the phase shifters may not be easy to be reduced [4–6].

For future works, the boundary condition, which was ignored in this paper for simplifying the discussion, must be investigated. The dependency of the proposed method on geometry of the adjacent pixels should be examined as well. Moreover, combining the proposed method with other PCA-based methods [21–24], or optimization methods as in [25,26], can be the next step to improve the overall performance.

## Funding

Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for JSPS Fellows (15J08043, 16J06772).

## References and links

**1. **D. Malacara, M. Servín, and Z. Malacara, *Interferogram Analysis For Optical Testing*, 2nd ed. (CRC, 2005). [CrossRef]

**2. **H. Schreiber and J. H. Bruning, “Phase shifting interferometry,” in *Optical Shop Testing*, 3rd ed., D. Malacara, ed. (John Wiley & Sons, Inc., 2007). [CrossRef]

**3. **I. Yamaguchi and T. Zhang, “Phase-shifting digital holography,” Opt. Lett. **22**(16), 1268–1270 (1997). [CrossRef] [PubMed]

**4. **Y. Awatsuji, M. Sasada, and T. Kubota, “Parallel quasi-phase-shifting digital holography,” Appl. Phys. Lett. **85**(6), 1069–1071 (2004). [CrossRef]

**5. **D. I. Serrano-García, A. Martínez-García, N.-I. Toto-Arellano, and Y. Otani, “Dynamic temperature field measurements using a polarization phase-shifting technique,” Opt. Eng. **53**(11), 112202 (2014). [CrossRef]

**6. **K. Ishikawa, K. Yatabe, N. Chitanont, Y. Ikeda, Y. Oikawa, T. Onuma, H. Niwa, and M. Yoshii, “High-speed imaging of sound using parallel phase-shifting interferometry,” Opt. Express **24**(12), 12922–12932 (2016). [CrossRef] [PubMed]

**7. **K. Okada, A. Sato, and J. Tsujiuchi, “Simultaneous calculation of phase distribution and scanning phase shift in phase shifting interferometry,” Opt. Commun. **84**(3), 118–124 (1991). [CrossRef]

**8. **C. T. Farrell and M. A. Player, “Phase step measurement and variable step algorithms in phase-shifting interferometry,” Meas. Sci. Technol. **3**(10), 953–958 (1992). [CrossRef]

**9. **I.-B. Kong and S.-W. Kim, “General algorithm of phase-shifting interferometry by iterative least-squares fitting,” Opt. Eng. **34**(1), 183–188 (1995). [CrossRef]

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

**11. **A. Patil and P. Rastogi, “Approaches in generalized phase shifting interferometry,” Opt. Lasers Eng. **43**(3–5)475–490 (2005). [CrossRef]

**12. **H. Guo, Y. Yu, and M. Chen, “Blind phase shift estimation in phase shifting interferometry,” J. Opt. Soc. Am. A. **24**, (1)25–33 (2007). [CrossRef]

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

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

**15. **F. Liu, Y. Wu, and F. Wu, “Correction of phase extraction error in phase shifting interferometry based on Lissajous figure and ellipse fitting technology,” Opt. Express **23**(8), 10794–10807 (2015). [CrossRef] [PubMed]

**16. **C. Meneses-Fabian and F. A. Lara-Cortes, “Phase retrieval by Euclidean distance in self-calibrating generalized phase-shifting interferometry of three steps,” Opt. Express **23**(10), 13589–13604 (2015). [CrossRef] [PubMed]

**17. **Y. Xu, Y. Wang, Y. Ji, H. Han, and W. Jin, “Three-frame generalized phase-shifting interferometry by a Euclidean matrix norm algorithm,” Opt. Lasers Eng. **84**, 89–95 (2016). [CrossRef]

**18. **X. Xu, X. Lu, J. Tian, J. Shou, D. Zheng, and L. Zhong, “Random phase shifting interferometry based on independent component analysis,” Opt. Commun. **370**, 75–80 (2016). [CrossRef]

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

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

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

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

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

**24. **J. Deng, K. Wang, D. Wu, X. Lv, C. Li, J. Hao, J. Qin, and W. Chen, “Advanced principal component analysis method for phase reconstruction,” Opt. Express **23**(9), 12222–12231 (2015). [CrossRef] [PubMed]

**25. **K. Yatabe and Y. Oikawa, “Convex optimization based windowed Fourier filtering with multiple windows for wrapped phase denoising,” Appl. Opt. **55**(17), 4632–4641 (2016). [CrossRef] [PubMed]

**26. **K. Yatabe, K. Ishikawa, and Y. Oikawa, “Compensation of fringe distortion for phase-shifting three-dimensional shape measurement by inverse map estimation,” Appl. Opt. **55**(22), 6017–6024 (2016). [CrossRef] [PubMed]