Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Demodulation of two-shot fringe patterns with random phase shifts by use of orthogonal polynomials and global optimization

Open Access Open Access

Abstract

We propose a simple and robust phase demodulation algorithm for two-shot fringe patterns with random phase shifts. Based on a smoothness assumption, the phase to be recovered is decomposed into a linear combination of finite terms of orthogonal polynomials, and the expansion coefficients and the phase shift are exhaustively searched through global optimization. The technique is insensitive to noise or defects, and is capable of retrieving phase from low fringe-number (less than one) or low-frequency interferograms. It can also cope with interferograms with very small phase shifts. The retrieved phase is continuous and no further phase unwrapping process is required. The method is expected to be promising to process interferograms with regular fringes, which are common in optical shop testing. Computer simulation and experimental results are presented to demonstrate the performance of the algorithm.

© 2016 Optical Society of America

1. Introduction

Phase-shifting interferometry (PSI) is a powerful technique in optical metrology to measure a wide variety of physical quantities, such as strain, temperature and surface deformation [1]. It generally collects several phase-shifted interferograms to recover encoded phase of the test quantity. Many algorithms have been developed to retrieve phase from multiple interferograms, such as three-frame, four- frame and five-frame algorithms [1]. However, these algorithms require at least three interterferograms with constant phase shifts having special and known values (e.g., π/3, π/2) to do the computation. These rigid requirements put critical demands on the stability of the employed phase shifter and environmental conditions.

To extend the versatility of PSI, a lot of interesting work trying to recover phase from interferograms with arbitrary phase shifts has been developed [2–7]. Chen et al. proposed a max-min algorithm to extract phase by comparing pixels within one interferogram and a set of interferograms [2]. Based on a least-square iterative procedure, Wang et al. used an advanced iteration algorithm (AIA) to retrieve phase from as few as three interferograms with unknown phase shifts [3]. Vargas et al. presented a principle component analysis-based method to compute two quadrature signals from three or more interferograms to estimate the phase [4]. They also developed an approach to obtain modulating phase from interfergrams affected by tilt aberrations [8]. All the algorithms mentioned above are able to cope with interferograms with unknown phase shifts, but still require at least three phase-shifted interferograms as input.

Generally, PSI suffers from several errors during interferograms acquisition, such as mechanical vibration, temperature fluctuation, and air turbulence. If all these errors are well controlled, it is practical for PSI to acquire sufficient number of interferograms to precisely reconstruct the encoded phase. Otherwise, fewer interferograms may help minimize induced errors. In principle, single-shot interferogram already contains all phase information to be measured. Several well-established techniques [9–24], such as Fourier transform [9], windowed Fourier transform (WFT) [10], wavelet transform [11], spiral phase transform (SPT) [12], Hilbert-Huang transform [18–20], regularized-phase tracking (RPT) [13, 14] and spatial carrier phase shifting [21–24], have already been successfully applied to single-shot fringe pattern demodulation. However, the problem of single-shot interferometry lies in that it cannot resolve sign ambiguity of the retrieved phase.

Considering all the facts above, two-shot phase-shifting interferometry serves as a good compromise between single-shot and multi-frame interferometry, and, therefore, has aroused a lot of interests [25–33]. Kreis et al. pioneered the work of phase retrieval from two-shot interferograms with unknown phase shifts and proposed a Fourier transform-based demodulation method [25]. The method works well when the phase shift ranges from π/3 to 2π/3, but may fail for small phase shifts or noisy fringe patterns. Vargas et al. proposed several practical algorithms for interferograms with arbitrary phase shifts, including the asynchronous self-tuning (AST) algorithm [26], the optical flow (OF) algorithm [27], and the Gram-Schmit (GS) orthogonalization algorithm [28]. Among them, the GS algorithm is simple and fast, and is able to give very accurate results. But it does not work well for interferograms with special phase shifts, i.e. 0 rad or π rad, in which cases the two fringe patterns become linearly dependent and orthonormalization is impossible. Ma et al. proposed a two-dimensional (2D) continuous wavelet transform (CWT) method for phase extraction from two-shot arbitrarily phase-shifted interferograms [29]. The algorithm is robust to noise and defects, but has the same problem as the GS algorithm, because of an involved interferogram orthonormalization process. Patorski and Trusiak also developed algorithms for amplitude demodulation from two-shot phase-shifted interferograms [34,35].

In this paper, we propose a two-shot phase demodulation algorithm for fringe patterns with arbitrary phase shifts. The algorithm is insensitive to noise or defects, and is able to extract very small phase shifts and retrieve phase from low fringe-number (less than one) interferograms. The principle, examples and discussions are presented below.

2. Principle

2.1 Demodulation of two-shot phase-shifted interferograms

Two-shot phase-shifted fringe patterns can be formulated as

I1(x,y)=a1(x,y)+b1(x,y)cos[ϕ(x,y)]+n1(x,y),
I2(x,y)=a2(x,y)+b2(x,y)cos[ϕ(x,y)+δ]+n2(x,y),
where x and y are the spatial coordinates; a1, a2 and b1, b2 are the background and the modulation terms, respectively; ϕ and δ are the phase map and the phase shift to be recovered, and n1, n2 are the additive noise. After normalization or pre-filtering [36–43], Eqs. (1) and (2) become
I1n(x,y)=cos[ϕ(x,y)],
I2n(x,y)=cos[ϕ(x,y)+δ],
where I1n and I2n are the normalized intensities.

Usually, it is reasonable to assume that the phase map ϕ is a smooth function, which is true for most cases, especially in optical shop testing. In this way, ϕ can be decomposed into a linear combination of finite terms of 2D orthogonal functions ψ, i.e.

ϕ(x,y)ϕf(x,y)=kckψk(x,y),
where ϕf is the fitted phase, k is the index, ck is the expansion coefficient of the orthogonal basis function ψk. Based on the least-square criteria, we construct a cost function to describe the closeness between the normalized fringe patterns and the estimates, i.e.
f(X)=(x,y)mask({I1n(x,y)cos[ϕf(x,y)]}2+{I2n(x,y)cos[ϕf(x,y)+δ]}2),
where mask is a 2D aperture that defines valid data region, X = [c1, c2, …, ck, δ]T is the solution vector, for which the cost function attains its minimum,
X=argminXf(X).
Practically, Eq. (7) can be evaluated through optimization and, whenever the solution X is estimated, it is substituted back into Eq. (5) to reconstruct the continuous phase map. Compared with Eq. (7) in [17], Eq. (6) excludes a correlation term, which was originally developed to improve the conditioning of the solution space of the cost function. This omission here is due to two facts that 1) one more fringe pattern is already added in the cost function to constrain the solution, and 2) simplifying the cost function speeds up the optimization.

Next, we will discuss the selection of orthogonal functions adopted to represent the phase map and the optimization algorithm used to seek the expansion coefficients.

2.2 Selection of 2D orthogonal functions

Theoretically, any function that is orthogonal on a predefined aperture can be employed as the basis. Many kinds of orthogonal functions are available to use, such as the Zernike polynomials [44], the Legendre polynomials [45,46] and the Chebyshev polynomials [45]. Since square aperture and circular aperture are two most common apertures in practice, we use the Legendre polynomials, which are orthogonal on a unit square, for interferograms with a square aperture, and the Zernike polynomials, which are orthogonal on a unit circle, for interferograms with a circular aperture (see Fig. 1). Both the two polynomials have corresponding relationship with classic aberrations [1,46]. If the shape of the data is neither circular nor square, the Gram-Schmidt technique [47] can be used to orthogonalize the Zernike or the Legendre polynomials to custom a set of orthogonal polynomials over the predefined aperture.

 figure: Fig. 1

Fig. 1 Definition of unit square and unit circle. Q(x, y) or Q(ρ, θ) is a point within the unit square and the unit circle.

Download Full Size | PDF

A. Legendre polynomials for interferograms with a square aperture

The one-dimensional (1D) non-orthonormal Legendre polynomials P are recursively defined as

Pm(x)=2m1mxPm1(x)m1mPm2(x)
where x ∈ [-1, 1] and the order m is a non-negative integer. The first two terms are P0(x) = 1, P1(x) = x. The 1D Legendre polynomials meet the following orthogonality condition
11Pm(x)Pm(x)dx=22m+1δmm
where δmm′ is the Kronecker delta.

The 2D non-orthonormal Legendre polynomials L of order n are defined by multiplying two 1D Legendre polynomials in the x and the y directions

Lj(x,y)=Pl(x)Pm(y)
where the order n = l + m, l and m are non-negative integers, and j = (n + 1)(n + 2)/2 is the index. The orthogonality of the 2D Legendre polynomials is
1111Lj(x,y)Lj(x,y)dxdy=22l+122m+1δllδmm
where Lj(x, y) = Pl(x)Pm(y), Lj′(x, y) = Pl′(x)Pm′(y). The 2D Legendre polynomials have the property of being limited to the range [-1, 1], i.e. Lj ∈ [-1, 1]. The Legendre polynomials up to 4th order and their relations with classic aberrations are given in Table 1 [46].

Tables Icon

Table 1. 2D Legendre polynomials up to 4th order

B. Zernike polynomials for interferograms with a circular aperture

The standard 2D non-orthonormal Zernike polynomials are defined as

Zj={Rnm(ρ),m=0,Rnm(ρ)cosmθ,m0andevenj,Rnm(ρ)sinmθ,m0andoddj,
where ρ ∈ [0, 1] and θ ∈ [0, 2π] are the normalized radial coordinate and the angular coordinate, respectively; j = (n + 1)(n + 2)/2 is the index, and the radial polynomials
Rnm(ρ)=s=0(nm)/2(1)s(ns)!s![(n+m)/2s]![(nm)/2s]!ρn2s,
where n and m are non-negative integers, which mean radial degree and azimuthal frequency, respectively, and satisfy nm = even. The orthogonality of the Zernike polynomials is
0102πZj(ρ,θ)Zj(ρ,θ)ρdρdθ=πδjj.
The non-orthonormal Zernike polynomials also have the property of being limited to the range [-1, 1], i.e. Zj ∈ [-1, 1]. The standard non-orthonormal Zernike polynomials up to 4th order and their relations with classic aberrations are given in Table 2 [1].

Tables Icon

Table 2. Standard non-orthonormal Zernike polynomials up to 4th order

2.3 Global optimization using the differential evolution algorithm

Usually, the constructed cost function [Eq. (6)] is highly multimodal and has numerous local minima. To find the optimal solution of the multidimensional problem, robust global optimization algorithms are essential. Here we use the differential evolution (DE) algorithm, which is a simple and efficient population-based stochastic search technique [48], to find the global optimum. It has been shown that, for most instances, the DE algorithm outperforms many other global optimization algorithms, such as genetic algorithm (GA) and simulated annealing (SA) in terms of required number of cost function evaluations (NFEs) necessary to locate the global minimum [48].

The DE algorithm works with parameter vectors having a population size NP (NP ≥ 3) and aims at evolving it towards the global minimum by repeated mutation, crossover and selection. The solution vector XG j (the jth individual in the population at the Gth generation) is denoted as

XjG=[c1,jG,c2,jG,,ck,jG,δjG]T=[x1,jG,x2,jG,,xD1,jG,xD,jG]T,
where D = k + 1 is the dimensionality, j = 1, 2, …, NP is the index of the individual in the population, NP is the population size and is normally in the range [5D, 20D]. A larger NP usually increases the probability of global convergence, but implies larger NFEs. The DE algorithm can be implemented through the following four steps.

Step 0: Initialization

The initial population should cover the entire search space as much as possible by uniformly randomizing the individuals. Denote the maximum and the minimum search bounds as Xmax = [x1max, x2max, …, ximax, …, xDmax]T and Xmin = [x1min, x2min, …, ximin, …, xDmin]T. The ith unknown in the jth individual at the 0th generation can be initialized as follows

xi,j0=ximin+rand(0,1)×(ximaxximin),
where i = 1, 2, 3, …, D, and rand(0, 1) is a random number uniformly distributed in [0, 1].

Step 1: Mutation

After initialization, a mutant vector is generated from the population using a mutation strategy DE/rand/1 (DE: differential evolution, rand: random selection of the vector to be mutated, 1: number of difference vectors used) [48], which can be expressed as

VjG=Xr1G+F(Xr2GXr3G),
where VG j = [vG 1,j, vG 2,j, …, vG D,j]T is the mutant vector; r1, r2 and r3 are three random integer indices uniformly selected in the range [1, NP] and satisfy r1r2r3j; F ∈ [0, 2] is a scaling factor, which controls the magnitude of the differential vector. In practice, a large F usually increases the probability of escaping a local minimum, but decreases the convergence speed. A dithering technique [49] that randomly selects F in a certain interval, for example, F = 0.3(1 + rand(0, 1)), for each generation is of great help to speed up the convergence.

Step 2: Crossover

Afte mutation, a trial vector can be generated by crossing over the target vector and the mutant vector using a binomial crossover operator, which is defined as

ui,jG={vi,jG,ifrandi(0,1)CRori=irand,xi,jG,otherwise,
where UG j = [uG 1,j, uG 2,j, …, uG D,j]T is the trial vector, randi(0, 1) is a random number in [0, 1] for the ith unknown, irand is a random integer index uniformly distributed in the range [1, D], and CR ∈ [0, 1] is the crossover rate. A large CR often helps the speed, but has the risk of premature convergence. The crossover operator copies the ith parameter of the mutant vector VG jto the same position of the trial vector UG j if randi(0, 1) ≤ CR or i = irand. Otherwise, it is copied from the target vector XG j. If the generated trial vector exceeds the search bounds Xmax and Xmin, it can be randomly and uniformly reinitialized according to Eq. (16). The crossover process is shown in Fig. 2.

 figure: Fig. 2

Fig. 2 Illustration of the crossover process. When randi(0, 1) ≤ CR or i = irand, copy the ith mutant element to the trial vector; Otherwise, copy the ith element in the target vector to the trial vector.

Download Full Size | PDF

Step 3: Selection

To determine whether the trial vector should become a new individual of the population at generation G + 1, a greedy criterion is used

XjG+1={UjG,iff(UjG)f(XjG)XjG,otherwise.
Whenever the cost function of the trial vector is no greater than that of the target vector, the trial vector replaces the target vector and enters the population of the next generation. Otherwise, it is discarded and the target vector remains.

Repeat step 1 to step 3 generation by generation until a preset convergence criterion are satisfied.

3. Simulation and experiments

To test the performance of the proposed algorithm, both numerical and real experiments were carried out. For all examples presented below, we used the following DE parameters: NP = 20D, CR = 1 and F = 0.3 [1 + rand(0, 1)].

3.1 Simulation

In the simulation, we tested the proposed algorithm from four different aspects, including, 1) high-order (n = 8) polynomials phase decomposition, 2) high-dimensional (D = 46) global optimization, 3) small phase shift (δ = 0.01π rad) extraction, 4) highly-noisy and defective fringe patterns demodulation.

Two-shot fringe patterns with a square aperture and 256 × 256 pixels [Figs. 3(a) and 3(b)] were simulated according to Eqs. (1) and (2) by setting a1 = a2 = 0 and b1 = b2 = 1. The interferograms were heavily corrupted with additive Gaussian white noise (AGWN) with a mean 0 and a standard deviation 1.2 and modulated by three diffraction rings. The theoretical phase [shown in Figs. 3(c) and 3(d) after wrapping] was modeled by the following expression

ϕ(x,y)=π×0.5{(x2+y2)10[cos(2πx)+cos(2πy)]},
where x, y ∈ [-0.9, 0.9]. The simulated phase shift δ between the two fringe patterns was a small value close to zero, i.e. 0.01π rad, which puts critical demands on the performance of the phase retrieval algorithm. Since δ is very small, the two-shot fringe patterns look very similar to each other and many two-shot algorithms will fail for such a case [25,26].

 figure: Fig. 3

Fig. 3 Simulated two-shot interferograms with a phase shift 0.01π. (a) and (b) two-shot interferograms with noise and defects; (c) and (d) 2D (wrapped) and 3D (unwrapped) true phase. Colorbar unit: rad.

Download Full Size | PDF

By using the Legendre polynomials up to 8th order and a search range ximax = -ximin = 10π, the cost function was evaluated 4.6 million times before the DE algorithm successfully located the global minimum. The retrieved 8th order Legendre coefficients are shown in Fig. 4(a). It is clear that, although the fringe patterns are very dense in the spatial domain, the expansion coefficients are quite sparse in the orthogonal domain and only a few terms play the dominant role. The finally reconstructed phase map and residual error map (with respect to the ground truth) are given in the first column of Fig. 5. The calculated quality index (Q index) of the reconstructed phase map and the root mean square (RMS) value of the error map are 0.908 and 0.139 rad, respectively, where the Q index [38,50] is defined as

Q=4σABμAμB(σA2+σB2)(μA2+μB2),
where A and B represent the ground truth and the estimated phase map; μA and μB are the mean values; 𝜎A and 𝜎B are the standard deviations; and 𝜎AB is the covariance of A and B. The dynamic range of the Q index is [-1, 1], where 1 means perfect match. The extracted phase shift in this case is 0.0098π rad, which is only 2% different from its true value 0.01π rad.

 figure: Fig. 4

Fig. 4 Extracted 8th order coefficients. (a) Legendre coefficients estimated by the proposed two-shot method, (b) Zernike coefficients estimated by the single-shot method in [17].

Download Full Size | PDF

 figure: Fig. 5

Fig. 5 Demodulation results of the two-shot interferograms in Fig. 3 using the proposed two-shot method, the single-shot method in [17], the Kreis method, the OF method and the GS method. First row: demodulated phase map, second row: corresponding error map.

Download Full Size | PDF

After intensity normalization [36] under the same conditions, the same interferograms were also evaluated using the single-shot method in [17], the Kreis method [25], the OF method [27] and the GS method [28]. The reconstructed phase and corresponding residual errors are shown in Fig. 5. Since the fringe patterns are very noisy, the solution space of the cost function is extremely complex and the single-shot method employing 8th order of Zernike polynomials failed to pinpoint the best solution and gave erroneous result. Figure 4(b) shows the estimated 45 Zernike coefficients retrieved by the single-shot method. The quantified Q indices of the reconstructed phase maps and the RMS values of the error maps by these methods are compared in Table 3. Among all the algorithms, the proposed method has the highest Q index and the smallest RMS value.

Tables Icon

Table 3. Summary of Q indices and RMS values of Figs. 5, 7 and 8

3.2 Experiments

Figure 6 shows two-shot phase extraction from two 256 × 256 pixels interferograms with open fringes by the proposed algorithm. The fringe patterns are simple in shape and 3rd order Zernike polynomials were used to represent the modulated phase. The optimal expansion coefficients estimated by the DE algorithm is shown in Fig. 6(c). The largest coefficients are the 2nd and the 3rd terms (Z2, Z3), which correspond to x-tilt and y-tilt (Table 2), respectively. The estimated phase shift in this case is −2.84 rad and the finally reconstructed phase map is shown in Fig. 6(d).

 figure: Fig. 6

Fig. 6 Phase extraction of two-shot interferograms with open fringes. (a) and (b) two experimentally acquired fringe patterns; (c) and (d) estimated 3rd order Zernike coefficients and reconstructed phase map (rewrapped for comparison).

Download Full Size | PDF

Figure 7 also shows phase extraction from two-shot 256 × 256 pixels interferograms but with closed fringes, which is very common in non-null interferometry, for example, in aspheric surfaces measurement [51]. Eleven-frame phase-shifted interferograms were first processed by the AIA method [3] and the recovered phase [Fig. 7(c)] was regarded as the reference. The proposed method employed up to 7th order Zernike polynomials to represent the phase. The evaluated expansion coefficients, the reconstructed phase map and the quantified error map (with respect to the AIA phase map) are shown in Figs. 7(d)-7(f), respectively. The Zernike coefficients mainly concentrate in the 4th and 11th terms, which corresponds to defocus and primary spherical aberrations, respectively. The quantified Q index and RMS value for this case are 0.865 and 0.236 rad, respectively. After the same normalization operation [36], the two-shot fringe patterns were also demodulated by the Kreis method [25] for comparison, and the extracted phase map and residual error are shown in Figs. 7(g) and 7(h), respectively. The quantified Q index and RMS value are 0.770 and 0.421 rad, respectively. Clearly, the phase retrieved by the proposed method is consistent with the AIA result and has a better quality than the result by the Kreis method.

 figure: Fig. 7

Fig. 7 Phase extraction of two-shot interferograms with closed fringes. (a) and (b) two experimental fringe patterns; (c) retrieved phase by the AIA method using 11 frame images serving as the reference; (d) – (f) estimated 7th order Zernike coefficients, reconstructed phase map and residual error map by the proposed method, respectively; (g) and (h) reconstructed phase map and residual error map by the Kreis method. Colorbar unit: rad.

Download Full Size | PDF

Finally, we demonstrate the low fringe-number interferograms demodulation capability of the proposed algorithm (see Fig. 8). Low fringe-number interferograms, which typically have less-than-one fringe (phase change smaller than 2π), are common in several cases, such as in null interfermetry to avoid retrace error [52], and in cat’s eye configuration for absolute testing [1]. Many algorithms in the literature fail to treat such kind of fringe patterns [53, 54]. Figures 8(a) and 8(b) show two experimentally acquired 256 × 256 pixels fringe patterns with less-than-one fringe. We also first evaluated eleven-frame phase-shifted interferograms using the AIA method [3] and treated the recovered phase [Fig. 8(c)] as the reference. Fifth-order Zernike polynomials and a search range ximax = -ximin = π were used for the DE optimization in the proposed algorithm. The retrieved Zernike coefficients, reconstructed continuous phase map and residual phase error (with respect to the AIA result) are shown in Figs. 8(d)-8(f), respectively. The quantified Q index and RMS value for this case are 0.930 and 0.171 rad, respectively. The phase extracted by the proposed method agrees with that of the AIA algorithm in both shape and magnitude. For comparison, after normalization under the same conditions, the interferograms were also evaluated using the GS method [28], and the unwrapped phase and residual error are shown in Figs. 8(g) and 8(h), respectively. The quantified Q index and RMS value are 0.963 and 0.365 rad. Although the GS method has slightly better Q index than the proposed method, the residual reconstruction error is bigger in magnitude. Since lower fringe-number interferograms usually have lower-spatial frequency and smaller phase magnitude, lower-order orthogonal polynomials and narrower search ranges could be used, which helps the optimization process. From this point of view, the proposed algorithm is well suited for this type of interferograms.

 figure: Fig. 8

Fig. 8 Phase retrieval of low fringe-number interferograms. (a) and (b) two experimental fringe patterns with less-than-one fringe; (c) retrieved phase by the AIA method using 11 frame images serving as the reference; (d) - (f) estimated 5th order Zernike coefficients, reconstructed phase map and residual error map by the proposed method, respectively; (g) and (h) reconstructed phase map and residual error map by the GS method. Colorbar unit: rad.

Download Full Size | PDF

Table 4 summarizes the parameters used in all four examples presented above.

Tables Icon

Table 4. Statistical results of all four experiments by the proposed methoda

4. Discussion and conclusion

Compared with the single-shot approach in [17], the proposed two-shot method has at least three advantages. First, the two-shot method is potentially capable of resolving the global sign ambiguity of the demodulated phase (i.e. ϕ or -ϕ). It is principally impossible to determine the global sign of the phase using the single-shot algorithm because of the even property of the cosine function, i.e. cos(ϕ) = cos(-ϕ). But in the two-shot algorithm, if we know the sign of the phase shift δ beforehand [i.e. δ ∈ (0, π) or δ ∈ (-π, 0)], the global sign ambiguity of the demodulated phase can be resolved because cos(ϕ + δ) ≠ cos(-ϕ + δ). All two-shot algorithms share the advantage over single-shot algorithms. Second, the two-shot algorithm is more robust and accurate. For both the two-shot algorithm and the single-shot method, they need to retrieve multiple unknown coefficients [(n + 1)(n + 2)/2 + 1 for the two-shot method or (n + 1)(n + 2)/2) for the single-shot method] through high-dimensional global optimization, which is not an easy task. The conditioning of the solution space of the cost function [Eq. (6) in this paper and Eq. (7) in [17]] is critical to the final convergence of the DE algorithm. Compared with that in the single-shot algorithm, the solution space conditioning in the two-shot method is improved because of tighter constraints by imposing one more phase-shifted image to the cost function. This enhancement makes the two-shot method more robust and accurate (see the comparison in Fig. 5), which is of particular importance for very noisy fringe patterns. Third, the two-shot method has the capability of demodulating fringe patterns with square apertures and small phase shifts, and registering the expansion coefficients with classic aberrations. This was also demonstrated in the numerical experiment (see Fig. 5).

We recommend using the Legendre polynomials for interferograms with square apertures and the Zernike polynomials for interferograms with circular apertures. In this way, the estimated values of the expansion coefficients are independent from the fitting order n and have direct relationships with classic aberrations [1, 46]. Alternatively, it is feasible to use the Legendre polynomials for circular apertures. But since they are not orthogonal over a unit disk, the expansion coefficients will not be independent any more. In other words, we may get totally different coefficients for different fitting order n. There is no problem if we do not care the physical meanings of these coefficients. But for some cases, like in optical shop testing, we may need to register the coefficients with classic aberrations (Tables 1 and 2). Inappropriate selection of the orthogonal polynomials may give erroneous results.

One fundamental assumption of the proposed algorithm is that the phase to be recovered is smooth and can be represented by finite terms of orthogonal functions with negligible truncation errors. This restricts the applications of the algorithm to interferograms with regular fringe shapes, which can be well modelled by low-order orthogonal polynomials (n ≤ 8 in this paper). For interferograms with highly-irregular (large phase gradient) or discontinuous fringes (step-like phase), e.g., in profilometry, it is possible to use a large number of orthogonal polynomials to approximate the phase function (just like Fourier approximation using the sine and cosine basis functions), but it may not be practical due to the high-dimensional global optimization process (for example, for n = 30, j = 496), in which the computation time goes up nonlinearly with the increase of the dimensionality j. Fortunately, most interferograms in PSI have regular and relatively low-frequency fringes, especially those in optical shop testing, where it is a standard practice to do wavefront fitting using 7th order Zernike or Legendre polynomials [1]. The proposed algorithm is expected to be promising for such cases.

For the determination of search bounds, as a general rule of thumb, it is a good start to set the values according to fringe numbers p in the interferogram, i.e. ximax = -ximin = pπ. Since both the Legendre and the Zernike polynomials limit themselves within the range [-1, 1], the preset search bounds guarantee the true solution is within the search space. In practice, some prior information, such as dominant aberrations or diminishing expansion coefficients with fitting order, could be of help to probably refine the search ranges to speed up the convergence.

The fitting order n discussed above determines how well the phase can be represented by finite terms of polynomials and the amount of the residual fitting error. The search bounds ximin and ximax determines whether the true solution is within the search space. Besides them, three parameters used in the DE algorithm, i.e. the population size NP, the scaling factor F and the crossover rate CR, which also need to be selected manually, determine whether the DE algorithm will converge to the optimal solution. Generally, NP ∈ [5D, 20D], F ∈ [0, 2] and CR ∈ [0, 1]. Larger NP and F usually increase the probability of global convergence, but implies larger NFEs, while a large CR helps the speed, but has the risk of premature convergence. We used a parameter combination NP = 20D, CR = 1 and F = 0.3 [1 + rand(0, 1)] throughout the paper and robustly obtained good results.

In conclusion, we presented a simple, straightforward two-shot phase demodulation algorithm based on phase representation by orthogonal functions and global optimization. The proposed phase demodulator is very robust to noise and local defects, and is particularly good at phase retrieval from low fringe-number or low-frequency interterferograms. It can also deal with two-shot interferograms with very small phase shifts. The technique is computationally expensive because of high-dimensional global optimization, and may not be suitable to process fringe patterns with highly-irregular or discontinuous fringes. The algorithm can be extended to retrieve phase from arbitrary-shot phase-shifted interferograms.

Acknowledgments

This work was supported by National Natural Science Foundation of China (NSFC) under Grant Nos. 60877043 and 61575061.

References and links

1. D. Malacara, Optical Shop Testing, 3rd ed. (Wiley, 2007).

2. X. Chen, M. Gramaglia, and J. A. Yeazell, “Phase-shifting interferometry with uncalibrated phase shifts,” Appl. Opt. 39(4), 585–591 (2000). [CrossRef]   [PubMed]  

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

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

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

6. J. Xu, Q. Xu, L. Chai, Y. Li, and H. Wang, “Direct phase extraction from interferograms with random phase shifts,” Opt. Express 18(20), 20620–20627 (2010). [CrossRef]   [PubMed]  

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

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

9. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]  

10. Q. Kemao, “Two-dimensional windowed Fourier transform for fringe pattern analysis: Principles, applications and implementations,” Opt. Lasers Eng. 45(2), 304–317 (2007). [CrossRef]  

11. J. Zhong and J. Weng, “Spatial carrier-fringe pattern analysis by means of wavelet transform: wavelet transform profilometry,” Appl. Opt. 43(26), 4993–4998 (2004). [CrossRef]   [PubMed]  

12. K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18(8), 1862–1870 (2001). [CrossRef]   [PubMed]  

13. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a two-dimensional regularized phase-tracking technique,” Appl. Opt. 36(19), 4540–4548 (1997). [CrossRef]   [PubMed]  

14. C. Tian, Y. Yang, D. Liu, Y. Luo, and Y. Zhuo, “Demodulation of a single complex fringe interferogram with a path-independent regularized phase-tracking technique,” Appl. Opt. 49(2), 170–179 (2010). [CrossRef]   [PubMed]  

15. R. Kulkarni and P. Rastogi, “Phase derivative estimation from a single interferogram using a Kalman smoothing algorithm,” Opt. Lett. 40(16), 3794–3797 (2015). [CrossRef]   [PubMed]  

16. C. Tian, Y. Yang, S. Zhang, D. Liu, Y. Luo, and Y. Zhuo, “Regularized frequency-stabilizing method for single closed-fringe interferogram demodulation,” Opt. Lett. 35(11), 1837–1839 (2010). [CrossRef]   [PubMed]  

17. C. Tian, Y. Yang, T. Wei, T. Ling, and Y. Zhuo, “Demodulation of a single-image interferogram using a Zernike-polynomial-based phase-fitting technique with a differential evolution algorithm,” Opt. Lett. 36(12), 2318–2320 (2011). [CrossRef]   [PubMed]  

18. M. Trusiak, K. Patorski, and K. Pokorski, “Hilbert-Huang processing for single-exposure two-dimensional grating interferometry,” Opt. Express 21(23), 28359–28379 (2013). [CrossRef]   [PubMed]  

19. M. Trusiak, K. Patorski, and M. Wielgus, “Hilbert-Huang processing and analysis of complex fringe patterns,” Proc. SPIE 9203, 92030K (2014).

20. K. Patorski, M. Trusiak, and K. Pokorski, “Single-shot two-channel Talbot interferometry using checker grating and Hilbert-Huang fringe pattern processing,” Proc. SPIE 9132, 91320Z (2014).

21. M. Kujawinska and J. Wojciak, “Spatial-carrier phase-shifting technique of fringe pattern analysis,” Proc. SPIE 61, 61–67 (1991). [CrossRef]  

22. M. Pirga and M. Kujawinska, “Two directional spatial-carrier phase-shifting method for analysis of crossed and closed fringe patterns,” Opt. Eng. 34(8), 2459–2466 (1995). [CrossRef]  

23. J. Xu, Q. Xu, and H. Peng, “Spatial carrier phase-shifting algorithm based on least-squares iteration,” Appl. Opt. 47(29), 5446–5453 (2008). [CrossRef]   [PubMed]  

24. Y. Du, G. Feng, H. Li, J. Vargas, and S. Zhou, “Spatial carrier phase-shifting algorithm based on principal component analysis method,” Opt. Express 20(15), 16471–16479 (2012). [CrossRef]  

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

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

27. J. Vargas, J. A. Quiroga, C. O. S. Sorzano, J. C. Estrada, and J. M. Carazo, “Two-step interferometry by a regularized optical flow algorithm,” Opt. Lett. 36(17), 3485–3487 (2011). [CrossRef]   [PubMed]  

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

29. J. Ma, Z. Wang, and T. Pan, “Two-dimensional continuous wavelet transform algorithm for phase extraction of two-step arbitrarily phase-shifted interferograms,” Opt. Lasers Eng. 55, 205–211 (2014). [CrossRef]  

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

31. 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]   [PubMed]  

32. F. Liu, Y. Wu, and F. Wu, “Phase shifting interferometry from two normalized interferograms with random tilt phase-shift,” Opt. Express 23(15), 19932–19946 (2015). [CrossRef]   [PubMed]  

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

34. K. Patorski and M. Trusiak, “Highly contrasted Bessel fringe minima visualization for time-averaged vibration profilometry using Hilbert transform two-frame processing,” Opt. Express 21(14), 16863–16881 (2013). [CrossRef]   [PubMed]  

35. K. Patorski, M. Trusiak, and T. Tkaczyk, “Optically-sectioned two-shot structured illumination microscopy with Hilbert-Huang processing,” Opt. Express 22(8), 9517–9527 (2014). [CrossRef]   [PubMed]  

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

37. M. Trusiak, M. Wielgus, and K. Patorski, “Advanced processing of optical fringe patterns by automated selective reconstruction and enhanced fast empirical mode decomposition,” Opt. Lasers Eng. 52, 230–240 (2014). [CrossRef]  

38. M. Trusiak, K. Patorski, and M. Wielgus, “Adaptive enhancement of optical fringe patterns by selective reconstruction using FABEMD algorithm and Hilbert spiral transform,” Opt. Express 20(21), 23463–23479 (2012). [CrossRef]   [PubMed]  

39. J. A. Quiroga, J. A. Gomez-Pedrero, and A. Garcia-Botella, “Algorithm for fringe pattern normalization,” Opt. Commun. 197(1-3), 43–51 (2001). [CrossRef]  

40. J. A. Guerrero, J. L. Marroquin, M. Rivera, and J. A. Quiroga, “Adaptive monogenic filtering and normalization of ESPI fringe patterns,” Opt. Lett. 30(22), 3018–3020 (2005). [CrossRef]   [PubMed]  

41. N. A. Ochoa and A. A. Silva-Moreno, “Normalization and noise-reduction algorithm for fringe patterns,” Opt. Commun. 270(2), 161–168 (2007). [CrossRef]  

42. M. B. Bernini, A. Federico, and G. H. Kaufmann, “Normalization of fringe patterns using the bidimensional empirical mode decomposition and the Hilbert transform,” Appl. Opt. 48(36), 6862–6869 (2009). [CrossRef]   [PubMed]  

43. K. Patorski and K. Pokorski, “Examination of singular scalar light fields using wavelet processing of fork fringes,” Appl. Opt. 50(5), 773–781 (2011). [CrossRef]   [PubMed]  

44. R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66(3), 207–211 (1976). [CrossRef]  

45. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables (Courier Corporation, 1964).

46. V. N. Mahajan, “Orthonormal aberration polynomials for anamorphic optical imaging systems with circular pupils,” Appl. Opt. 51(18), 4087–4091 (2012). [CrossRef]   [PubMed]  

47. D. Malacara-Hernandez, M. Carpio-Valadez, and J. J. Sanchez-Mondragon, “Wavefront fitting with discrete orthogonal polynomials in a unit radius circle,” Opt. Eng. 29(6), 672–675 (1990). [CrossRef]  

48. R. Storn and K. Price, “Differential Evolution – A Simple and Efficient Heuristic for global Optimization over Continuous Spaces,” J. Glob. Optim. 11(4), 341–359 (1997). [CrossRef]  

49. S. Das, A. Konar, and U. K. Chakraborty, “Two improved differential evolution schemes for faster global search,” in Proceedings of the 7th annual conference on Genetic and evolutionary computation, (ACM, 2005), pp. 991–998. [CrossRef]  

50. Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE Signal Process. Lett. 9(3), 81–84 (2002). [CrossRef]  

51. C. Tian, Y. Yang, and Y. Zhuo, “Generalized data reduction approach for aspheric testing in a non-null interferometer,” Appl. Opt. 51(10), 1598–1604 (2012). [CrossRef]   [PubMed]  

52. C. Tian, Y. Yang, T. Wei, and Y. Zhuo, “Nonnull interferometer simulation for aspheric testing based on ray tracing,” Appl. Opt. 50(20), 3559–3569 (2011). [CrossRef]   [PubMed]  

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

54. C. Luo, L. Zhong, P. Sun, H. Wang, J. Tian, and X. Lu, “Two-step demodulation algorithm based on the orthogonality of diamond diagonal vectors,” Appl. Phys. B 119(2), 387–391 (2015). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (8)

Fig. 1
Fig. 1 Definition of unit square and unit circle. Q(x, y) or Q(ρ, θ) is a point within the unit square and the unit circle.
Fig. 2
Fig. 2 Illustration of the crossover process. When randi(0, 1) ≤ CR or i = irand, copy the ith mutant element to the trial vector; Otherwise, copy the ith element in the target vector to the trial vector.
Fig. 3
Fig. 3 Simulated two-shot interferograms with a phase shift 0.01π. (a) and (b) two-shot interferograms with noise and defects; (c) and (d) 2D (wrapped) and 3D (unwrapped) true phase. Colorbar unit: rad.
Fig. 4
Fig. 4 Extracted 8th order coefficients. (a) Legendre coefficients estimated by the proposed two-shot method, (b) Zernike coefficients estimated by the single-shot method in [17].
Fig. 5
Fig. 5 Demodulation results of the two-shot interferograms in Fig. 3 using the proposed two-shot method, the single-shot method in [17], the Kreis method, the OF method and the GS method. First row: demodulated phase map, second row: corresponding error map.
Fig. 6
Fig. 6 Phase extraction of two-shot interferograms with open fringes. (a) and (b) two experimentally acquired fringe patterns; (c) and (d) estimated 3rd order Zernike coefficients and reconstructed phase map (rewrapped for comparison).
Fig. 7
Fig. 7 Phase extraction of two-shot interferograms with closed fringes. (a) and (b) two experimental fringe patterns; (c) retrieved phase by the AIA method using 11 frame images serving as the reference; (d) – (f) estimated 7th order Zernike coefficients, reconstructed phase map and residual error map by the proposed method, respectively; (g) and (h) reconstructed phase map and residual error map by the Kreis method. Colorbar unit: rad.
Fig. 8
Fig. 8 Phase retrieval of low fringe-number interferograms. (a) and (b) two experimental fringe patterns with less-than-one fringe; (c) retrieved phase by the AIA method using 11 frame images serving as the reference; (d) - (f) estimated 5th order Zernike coefficients, reconstructed phase map and residual error map by the proposed method, respectively; (g) and (h) reconstructed phase map and residual error map by the GS method. Colorbar unit: rad.

Tables (4)

Tables Icon

Table 1 2D Legendre polynomials up to 4th order

Tables Icon

Table 2 Standard non-orthonormal Zernike polynomials up to 4th order

Tables Icon

Table 3 Summary of Q indices and RMS values of Figs. 5, 7 and 8

Tables Icon

Table 4 Statistical results of all four experiments by the proposed methoda

Equations (21)

Equations on this page are rendered with MathJax. Learn more.

I 1 (x,y)= a 1 (x,y)+ b 1 (x,y)cos[ϕ(x,y)]+ n 1 (x,y),
I 2 (x,y)= a 2 (x,y)+ b 2 (x,y)cos[ϕ(x,y)+δ]+ n 2 (x,y),
I 1 n (x,y)=cos[ϕ(x,y)],
I 2 n (x,y)=cos[ϕ(x,y)+δ],
ϕ(x,y) ϕ f (x,y)= k c k ψ k (x,y) ,
f(X)= (x,y)mask ( { I 1 n (x,y)cos[ ϕ f (x,y)] } 2 + { I 2 n (x,y)cos[ ϕ f (x,y)+δ] } 2 ) ,
X= argmin X f(X).
P m (x)= 2m1 m x P m1 (x) m1 m P m2 (x)
1 1 P m (x) P m (x) dx= 2 2m+1 δ m m
L j (x,y)= P l (x) P m (y)
1 1 1 1 L j (x,y) L j (x,y)dxdy= 2 2l+1 2 2m+1 δ l l δ m m
Z j ={ R n m (ρ),m=0, R n m (ρ)cosmθ,m0andevenj, R n m (ρ)sinmθ,m0andoddj,
R n m (ρ)= s=0 (nm)/2 (1) s (ns)! s![(n+m)/2s]![(nm)/2s]! ρ n2s ,
0 1 0 2π Z j (ρ,θ) Z j (ρ,θ)ρdρdθ =π δ j j .
X j G = [ c 1,j G , c 2,j G ,, c k,j G , δ j G ] T = [ x 1,j G , x 2,j G ,, x D1,j G , x D,j G ] T ,
x i,j 0 = x imin +rand(0,1)×( x imax x imin ),
V j G = X r 1 G +F( X r 2 G X r 3 G ),
u i,j G ={ v i,j G ,if rand i (0,1)CRori= i rand , x i,j G ,otherwise,
X j G+1 ={ U j G ,iff( U j G )f( X j G ) X j G ,otherwise .
ϕ(x,y)=π×0.5{ ( x 2 + y 2 )10[ cos(2πx)+cos(2πy) ] },
Q= 4 σ AB μ A μ B ( σ A 2 + σ B 2 )( μ A 2 + μ B 2 ) ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.