Random phase retrieval approach using Euclidean matrix norm of sum and difference map and fast least-squares algorithm

Optics Communications 460 (2020) 125174

Contents lists available at ScienceDirect

Optics Communications journal homepage: www.elsevier.com/locate/optcom

Random phase retrieval approach using Euclidean matrix norm of sum and difference map and fast least-squares algorithm Yu Zhang ∗ Institute of Materials Physics, College of Science, Northeast Electric Power University, Jilin, Jilin 132012, China State Key Laboratory of Applied Optics, Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences, Changchun, Jilin 130022, China

ARTICLE

INFO

Keywords: Phase shifting algorithm Euclidean matrix norm of sum and difference map Fast least-squares algorithm Pre-filtering

ABSTRACT In order to meet the requirements of the in-situ optical measurement, a random phase retrieval approach using the Euclidean matrix norm of sum and difference map and fast least-squares algorithm (EMNSD&FLSA) is proposed. It can obtain the high accuracy as the iterative algorithm and cost less computational time as the non-iterative algorithm simultaneously. No pre-filtering, iteration and relatively accurate initial phase lead to the high accuracy, fast iterative method and only two phase shifted interferograms lead to the high speed. It can obtain absolutely accurate result when the background intensity and modulation amplitude are perfect. Moreover, it is suitable for different levels of noises. If the high accuracy is required, it would be best to suppress the noise before using EMNSD&FLSA, and the phase shift would be best to far away from 0 rad and 𝜋 rad, furthermore, the numbers of patterns are best to be more than 2. Last but not least, it is robust for the circular, straight or complex fringes. The simulations and experiments verify the correctness and feasibility of EMNSD&FLSA.

1. Introduction With the development of optics manufacturing technology, in-situ measurement technology has also been widely developed, and the interferometry is an easy, high-speed, accurate optical testing tool which is always used in the in-situ measurement. Since the optical phase distribution can be easily extracted by several interferograms, the phase shifting interferometry (PSI) has been widely used in optical measurement [1–3]. In order to ensure the efficiency of in-situ measurement, both the accuracy and working time of PSI need to be considered. The performance of PSI mainly depends on the interferometer, environment and phase shifting algorithm (PSA), for the fixed interferometer and environment, outstanding PSA can suppress the different kinds of errors to improve the accuracy of PSI [4–6], and the PSA with random phase shift will be insensitive to the phase shift error due to the miscalibration of piezo-transducer (PZT), vibrational error, air turbulence, instability of the laser frequency, moreover, the PSA which is easy to implement and only needs less number of the phase shifted interferograms will cost less time. Conventional PSA needs a sequence of three, four or sometimes more interferograms with the fixed, known or unknown phase shifts [4– 18], in which the accuracy and computational time are affected by the more interferograms. In recent years, many PSAs with only two phase shifted interferograms have been reported, these methods have been found to be a good compromise between the phase retrieval

accuracy and computational time, moreover, the phase shift between the two frames can be random and unknown except 0 rad and 𝜋 rad, these features make the two-step PSA much more accurate and insensitive to the phase shift error, and compared to multi-step PSAs, two-step PSAs cost less time in both capturing the interferograms and computing the phase. In 1992, Farrell and Player [19] utilized Lissajous figures and ellipse fitting to calculate the phase difference between two interferograms, and in 2016, Liu et al. [20] proposed a PSA which can simultaneously extract the tested phase and phase shift from only two interferograms using Lissajous figure and ellipse fitting technology, but these two algorithms both need pre-filtering. In 2012, [21] presented a two-step demodulation based on the Gram– Schmidt orthonormalization method (GS), it requires subtracting the DC term by filtering before performing GS, in the same year, [22] presented a two-step demodulation algorithm based on extreme value of the interference (EVI), the DC component also needs to be filtered out by a high-pass filter before performing EVI. In 2015, [23] proposed an advanced two-step phase demodulation algorithm based on the orthogonality of diamond diagonal vectors (DDV), and in the same year, [24] proposed a two-step PSA based on the quotient of inner products of phase shifted interferograms (QIP), only the cosine of the phase shift can be obtained, if the phase shift is more than 𝜋, the accurate phase distribution cannot be obtained, the above two PSAs also need pre-filtering. In 2018, [25] introduced a fast and accurate

∗ Correspondence to: Institute of Materials Physics, College of Science, Northeast Electric Power University, Jilin, Jilin 132012, China. E-mail address: [email protected].

https://doi.org/10.1016/j.optcom.2019.125174 Received 29 October 2019; Received in revised form 16 December 2019; Accepted 22 December 2019 Available online 24 December 2019 0030-4018/© 2019 Elsevier B.V. All rights reserved.

Y. Zhang

Optics Communications 460 (2020) 125174

where ⟨⋅⟩ means the average over all the pixels, then the average of background intensity can be calculated as: ⟨ 𝑠𝑢𝑚 ⟩ 𝐼 . (5) 𝑎 ≈ ⟨𝑎𝑘 ⟩ = 𝑘 2 We can get the sum of 𝐼1𝑘 and 𝐼2𝑘 after filtering, ( ) ( ) 𝜃 𝜃 𝐼̃𝑘𝑠𝑢𝑚 ≈ 𝐼 𝑠𝑢𝑚 − 2𝑎 = 2𝑏𝑘 cos 𝜑𝑘 + cos . (6) 2 2

wavefront reconstruction method for two-frame PSI, it also uses the high-pass Gaussian filter firstly before perform the phase retrieval, the cosine value of the unknown phase shift between two interferograms is estimated directly by solving a quartic polynomial equation, and then the phase map is readily reconstructed. From the above literatures, we found a phenomenon that many two-step PSAs need pre-filtering which will cost more time and affect the accuracy before the phase extraction, these PSAs are not suitable for high precision in-situ measurement, in addition, some PSAs can only obtain the cosine value of phase shift, hence these PSAs are not real random PSAs since the range of phase shift is limited between 0 and 𝜋. Although the accuracy of iterative PSA is relatively high, general iterative PSAs need more than three interferograms, and most of them use a phase shift value which is set randomly as the initial value, if the initial value is far away from the real phase shift, the inaccurate phase may be obtained, and the computational time may be longer, even more serious is that the iterative PSA is out of work [7–9], few papers can be found to use the phase distribution as the initial value because it is hard to obtain or set. To meet the requirements of high accuracy and less working time for the in-situ measurement, a fast iterative PSA with two phase shifted interferograms and random phase shift is essential, the background intensity would be best to take part in the calculation, and it is best to set the relatively accurate phase distribution as the initial value. In this paper, we will discuss a two-step random and iterative PSA with no pre-filtering. Section 2 presents the principle and process of the proposed PSA based on the Euclidean matrix norm of sum and difference map and fast least-squares algorithm (EMNSD&FLSA). In Section 3 the simulation of EMNSD&FLSA is discussed, and the comparison with EMNSD and GS is performed. Section 4 evaluates the novel PSA with the experimental data. The conclusion is finally drawn in Section 5.

The Euclidean matrix norm of 𝐼̃𝑘𝑠𝑢𝑚 and 𝐼̃𝑘𝑑𝑖𝑓 can be written as } ) ( )]2 1∕2 ( 𝜃 𝜃 cos 2𝑏𝑘 cos 𝜑𝑘 + 2 2 𝑘=1 {𝐾 }1∕2 ( ) ∑ [ √ ( )] 𝜃 2 = 2 cos 𝑏𝑘 1 + cos 2𝜑𝑘 + 𝜃 . 2 𝑘=1 {

‖𝐼̃𝑠𝑢𝑚 ‖ = ‖ ‖

} 𝐾 [ ) ( )]2 1∕2 ( ∑ 𝜃 𝜃 sin 2𝑏𝑘 sin 𝜑𝑘 + 2 2 𝑘=1 {𝐾 }1∕2 ( ) ∑ [ √ ( )] 𝜃 2 = 2 sin 𝑏𝑘 1 − cos 2𝜑𝑘 + 𝜃 . 2 𝑘=1

‖ ̃𝑑𝑖𝑓 ‖ ‖𝐼 ‖ = ‖ ‖

𝐾 ∑ 𝑘=1

{

𝑏2𝑘 ≫

𝐾 ∑

(8)

( ) 𝑏2𝑘 cos 2𝜑𝑘 + 𝜃 .

(9)

𝑘=1

Then Eqs. (7) and (8) can be simplified as { 𝐾 }1∕2 ( ) ( ) ∑ √ 𝜃 ‖𝐼̃𝑠𝑢𝑚 ‖ ≈ 2 cos 𝜃 𝑏2𝑘 = 𝐶 cos . ‖ ‖ 2 2 𝑘=1

2.1. PSA based on Euclidean matrix norm of sum and difference map (EMNSD)

( ) 𝜃 ‖ ̃𝑑𝑖𝑓 ‖ √ ‖𝐼 ‖ ≈ 2 sin ‖ ‖ 2

In PSI, two phase shifted interferograms with total pixels of K can be described as ( ) 𝐼1𝑘 = 𝑎1𝑘 + 𝑏1𝑘 cos 𝜑𝑘 ( ). (1) 𝐼2𝑘 = 𝑎2𝑘 + 𝑏2𝑘 cos 𝜑𝑘 + 𝜃

{

𝐾 ∑

}1∕2 𝑏2𝑘

= 𝐶 sin

𝑘=1

( ) 𝜃 . 2

}1∕2 √ {∑𝐾 2 where 𝐶 = 2 . 𝑘=1 𝑏𝑘 Then, we can obtain the phase by Eqs. (3) and (6), ) ( 𝐼̃𝑘𝑑𝑖𝑓 𝜑𝑘 = arctan cot (𝜃∕2) ⋅ 𝑠𝑢𝑚 − 𝜃∕2. 𝐼̃

where 𝑘 = 1, 2, …, 𝐾 denotes the pixel position, 𝐼1𝑘 and 𝐼2𝑘 are the intensity of two phase shifted interferograms, 𝑎1𝑘 , 𝑎2𝑘 , 𝑏1𝑘 and 𝑏2𝑘 respectively represent the background intensity and modulation amplitude of the two phase shifted interferograms, 𝜑𝑘 is the tested phase, and 𝜃 is the phase shift value between two phase shifted interferograms. Generally for the background intensity and modulation amplitude distributions, both the fluctuation between different interferograms and the non-uniformity between different pixels exist, however, the subtraction can still filter most of the background intensity, for simplicity, we assume that the background intensity and amplitude modulation are irrelevant to the image index, only relevant to the pixel position in this sub-section, so 𝑎1𝑘 = 𝑎2𝑘 = 𝑎𝑘 , 𝑏1𝑘 = 𝑏2𝑘 = 𝑏𝑘 . Here, the sum and difference of 𝐼1𝑘 and 𝐼2𝑘 can be defined as: ( ) ( ) 𝜃 𝜃 𝑠𝑢𝑚 𝐼𝑘 = 𝐼1𝑘 + 𝐼2𝑘 = 2𝑎𝑘 + 2𝑏𝑘 cos 𝜑𝑘 + cos . (2) 2 2

(10)

(11)

(12)

𝑘

cot (𝜃∕2) can be calculated by Eqs. (10) and (11), ‖𝐼̃𝑠𝑢𝑚 ‖ ‖. cot (𝜃∕2) = ‖ ‖𝐼̃𝑑𝑖𝑓 ‖ ‖ ‖ and 𝜃∕2 can be easily calculated. Hence, the phase can be calculated by Eq. (14), ( ) ( ‖ ̃𝑠𝑢𝑚 ‖ ) ‖𝐼̃𝑠𝑢𝑚 ‖ 𝐼̃𝑑𝑖𝑓 𝐼 ‖ 𝑘 ‖ ‖ −1 𝜑𝑘 = tan ⋅ 𝑠𝑢𝑚 − tan−1 ‖ . ‖𝐼̃𝑑𝑖𝑓 ‖ 𝐼̃ ‖𝐼̃𝑑𝑖𝑓 ‖ ‖ ‖ ‖ ‖ 𝑘

(13)

(14)

2.2. Fast least-squares algorithm (FLSA) There are too more approximations in the above phase retrieval, the phase distribution extracted by Eq. (14) will be not accurate, hence, in the following, we will utilize an iterative method to obtain more accurate result. In Ref. [26], we have ever researched a method called fast iterative algorithm (FIA), it can obtain the relatively accurate phase by three interferograms, it costs relatively less time because only a limited number of pixels take part in the iterative process, but the initial value of the phase shift is set randomly which may affect the accuracy and speed of the iteration, the background intensity does not need to be calculated because we ignore the fluctuation of background intensity between different interferograms, however this is not the case, hence, we improve FIA, we design a method which has a good initial

( ) ( ) 𝜃 𝜃 𝐼̃𝑘𝑑𝑖𝑓 = 𝐼1𝑘 − 𝐼2𝑘 = 2𝑏𝑘 sin 𝜑𝑘 + sin . (3) 2 2 Then, we perform the filtering for Eq. (2), the filtering algorithm may introduce extra error and cost more time, hence we utilize a simple method to filter 𝐼𝑘𝑠𝑢𝑚 . ⟨ ( )⟩ 𝜃 2

(7)

where ‖⋅‖ represents the Euclidean matrix norm. If there is more than one fringe, the following approximation can be applied,

2. Principles

If there is more than one fringe, cos 𝜑𝑘 +

𝐾 [ ∑

≈ 0, hence

⟨ ( )⟩ ( ) ⟨ 𝑠𝑢𝑚 ⟩ 𝜃 𝜃 𝐼𝑘 = ⟨𝐼1𝑘 + 𝐼2𝑘 ⟩ = ⟨2𝑎𝑘 ⟩ + 2𝑏𝑘 cos 𝜑𝑘 + cos = ⟨2𝑎𝑘 ⟩ . 2 2 (4) 2

Y. Zhang

Optics Communications 460 (2020) 125174

Step 3: We also choose the same limited number of samples from Eq. (16), and we use the same method ( ) as Step 2 to calculate ( )𝑎2 and 𝑏2 , the only difference is 𝜂2 = 𝑏2 cos 𝜃2 , and 𝜉2 = −𝑏2 sin 𝜃2 , hence 𝜃2 can be calculated by Eq. (26) ( ) 𝜉 𝜃 = tan−1 − 2 . (26) 2 𝜂2

value, and it can accurately calculate the background intensity, last but not least, it can obtain accurate phase with less computational time by only two phase shifted interferograms, it is called fast least-squares algorithm (FLSA). In the following, we will introduce it in detail: Step 1: We can rewrite Eq. (1) as ) ( ( ) 𝜃 𝜃 𝐼1𝑘 = 𝑎1𝑘 + 𝑏1𝑘 cos 𝜑𝑘 = 𝑎1𝑘 + 𝑏1𝑘 cos 𝜑𝑘 + − 2 2 ( ) 𝜃 = 𝑎1𝑘 + 𝑏1𝑘 cos 𝛷𝑘 − . (15) 2

Because we can calculate 𝜃2 by steps 2 and 3, to ensure the accuracy of 𝜃 , we calculate 𝜃2 by averaging the results of Eqs. (25) and (26). 2 Step 4: In order to calculate the accurate phase, we rewrite Eqs. (2) and (3) for the chosen pixels as Eqs. (27) and (28), note that the non-uniformity of the background intensity and modulation amplitude are ignored, ( ) ( ) ( ) 𝜃 𝜃 𝐼𝑗𝑠𝑢𝑚 = 𝐼1𝑗 + 𝐼2𝑗 = 𝑎1 + 𝑎2 + 𝑏1 + 𝑏2 cos 𝜑𝑗 + cos 2 2 ( ) ( ) ( ) 𝜃 𝜃 + 𝑏1 − 𝑏2 sin 𝜑𝑗 + sin . (27) 2 2 ( ) ( ) ( ) 𝜃 𝜃 sin 𝐼̃𝑗𝑑𝑖𝑓 = 𝐼1𝑗 − 𝐼2𝑗 = 𝑎1 −𝑎2 + 𝑏1 + 𝑏2 sin 𝜑𝑗 + 2 2 ( ) ( ) ) ( 𝜃 𝜃 + 𝑏1 − 𝑏2 cos 𝜑𝑗 + cos . (28) 2 2 ( ) We 𝑏1 − 𝑏2 approximates ) that ( )calculate ( ) (zero, ) then 𝑏1 − 𝑏2 sin ( assume ( ) 𝜑𝑗 + 𝜃2 sin 𝜃2 and 𝑏1 − 𝑏2 cos 𝜑𝑗 + 𝜃2 cos 𝜃2 both close to zero,

) ( ( ) 𝜃 𝜃 𝐼2𝑘 = 𝑎2𝑘 + 𝑏2𝑘 cos 𝜑𝑘 + 𝜃 = 𝑎2𝑘 + 𝑏2𝑘 cos 𝜑𝑘 + + 2 2 ( ) 𝜃 = 𝑎2𝑘 + 𝑏2𝑘 cos 𝛷𝑘 + . (16) 2 Step 2: We select a limited number of samples at regular intervals from 𝛷k which can be calculated by EMNSD as the initial phase, ) ( ‖𝐼̃𝑠𝑢𝑚 ‖ 𝐼̃𝑑𝑖𝑓 𝜃 ‖⋅ 𝑘 . (17) 𝛷𝑘 = 𝜑𝑘 + = tan−1 ‖ ‖𝐼̃𝑑𝑖𝑓 ‖ 𝐼̃𝑠𝑢𝑚 2 ‖ ‖ 𝑘 and we select the same samples from Eq. (15), the expression of first phase shifted signal with the chosen samples is ( ) 𝜃 𝐼 𝑡 1𝑗 = 𝑎1𝑗 + 𝑏1𝑗 cos 𝛷𝑗 − . (18) 2 where 𝑗 = 1, 2, …, 𝑁 denotes the chosen pixel number with N the total number of chosen pixels. Provided that the background intensity 𝑎1𝑗 and modulation amplitude 𝑏1𝑗 are irrelevant to j, so 𝑎11 =(𝑎12) = ⋯ = 𝑎1𝑁 = 𝑎( 1 , 𝑏)11 = 𝑏12 =

⋯ = 𝑏1𝑁 = 𝑏1 . By setting 𝜂1 = 𝑏1 cos becomes

𝜃 2

, and 𝜉1 = 𝑏1 sin

𝜃 2

finally Eqs. (27) and (28) can be simplified as ( ) ( ) ( ) 𝜃 𝜃 𝐼𝑗𝑠𝑢𝑚 = 𝑎1 + 𝑎2 + 𝑏1 + 𝑏2 cos 𝜑𝑗 + cos . 2 2 ( ) ( ) ( ) 𝜃 𝜃 𝐼̃𝑗𝑑𝑖𝑓 = 𝑎1 −𝑎2 + 𝑏1 + 𝑏2 sin 𝜑𝑗 + sin . 2 2 Finally, the new phase can be calculated by ( ) ⎛ ( ) 𝐼̃𝑑𝑖𝑓 − 𝑎1 − 𝑎2 ⎞ 𝑗 𝜃 𝛷𝑗 = 𝜑𝑗 + 𝜃∕2 = tan−1 ⎜cot ⋅ 𝑠𝑢𝑚 ( )⎟. ⎜ 2 𝐼 𝑗 − 𝑎1 + 𝑎2 ⎟ ⎠ ⎝

, Eq. (18)

𝐼 𝑡 1𝑗 = 𝑎1 + 𝜂1 cos𝛷𝑗 + 𝜉1 sin𝛷𝑗 .

(19)

The sum of squared differences between the theoretical and actual value of the phase shifted signal can be expressed as 𝑆1 =

𝑁 𝑁 ∑ ( ( 𝑡 )2 )2 ∑ 𝑎1 + 𝜂1 cos𝛷𝑗 + 𝜉1 sin𝛷𝑗 − 𝐼1𝑗 . 𝐼 1𝑗 − 𝐼1𝑗 =

where 𝐼1𝑗 is actual value of the phase shifted signal with the chosen samples obtained by the experimental data. According to the least-squares theory [7–9], 𝑆1 should be minimum, for the known 𝛷𝑗 , 𝜕𝑆1 ∕𝜕𝑎1 = 0, 𝜕𝑆1 ∕𝜕𝜂1 = 0, 𝜕𝑆1 ∕𝜕𝜉1 = 0, so 𝑋1 = 𝑆1−1 𝑅1 .

(21)

⎡ ⎢ 𝑁 ⎢ ⎢𝑁 ( ) ⎢∑ 𝑆1 = ⎢ cos 𝛷𝑗 ⎢ 𝑗=1 ⎢∑ 𝑁 ( ) ⎢ sin 𝛷𝑗 ⎢ ⎣ 𝑗=1 [ 𝑋1 = 𝑎1

𝑅1 =

[𝑁 ∑ 𝑗=1

𝜂1

𝐼1𝑗

𝜉1

]𝑇

𝑁 ∑ 𝑗=1

𝑁 ∑ 𝑗=1 𝑁 ∑

( ) cos2 𝛷𝑗

𝑗=1 𝑁 ∑

( ) ( ) sin 𝛷𝑗 cos 𝛷𝑗

𝑗=1

⎤ ⎥ ⎥ 𝑗=1 ⎥ 𝑁 ∑ ( ) ( )⎥ sin 𝛷𝑗 cos 𝛷𝑗 ⎥ . ⎥ 𝑗=1 ⎥ 𝑁 ∑ ( ) ⎥ sin2 𝛷𝑗 ⎥ 𝑗=1 ⎦ 𝑁 ∑

( ) cos 𝛷𝑗

( ) sin 𝛷𝑗

( ) 𝐼1𝑗 cos 𝛷𝑗

𝑁 ∑

]𝑇 ( ) 𝐼1𝑗 sin 𝛷𝑗 .

(24)

𝑗=1

𝜂1 and 𝜉1 can be obtained by Eq. (21), and ( ) 𝜉1 𝜃 = tan−1 . 2 𝜂1

𝜃 2

(31)

Note that, the unwrapped process is needed in the final step, and this is the only unwrapped process in the proposed method. Although EMNSD&FLSA is an iterative algorithm, it costs less time since it only needs two phase shifted interferograms, and only a limited number of samples are chosen to take part in the iterative process, moreover, only one time of unwrapped process is used in the whole calculation, this timesaving method is often used in the iterative algorithm, in addition, this algorithm does not need ‘‘real’’ filtering which costs more time. Except for the timesaving performance, EMNSD&FLSA is a accurate algorithm due to the iterative process and no pre-filtering. In the practical situation, generally for the background intensity and modulation amplitude distributions, both the fluctuation between different interferograms and the non-uniformity between different pixels exist, in Section 2.1, the fluctuation between different interferograms is ignored, but the phase extracted by EMNSD is only set as the initial value of the iteration, hence it will not affect the accuracy, then, in Section 2.2, the non-uniformity between different pixels is ignored, it will affect the accuracy, but all two-step PSAs cannot calculate the background intensity and modulation amplitude of every pixel because

(22)

(23)

.

(30)

From Eq. (31), we find that 𝑏1 and 𝑏2 are useless under the assumption that the background intensity and modulation amplitude are uniform for the different pixels, hence in the following calculation, we only need to calculate 𝑎1 , 𝑎2 and 𝜃2 . ( ) Step 5: Repeat steps 2 to 4 until 𝑅𝑀𝑆 𝛷𝑙 − 𝛷𝑙−1 < 𝜀, the iteration 𝜃 terminates, and the accurate 2 can be obtained, where 𝜀 is the predefined converging threshold of iteration, i.e., 10−5 rad, and l represents the number of iterations. Step 6: Perform step 4 using the whole samples of Eqs. (27) and (28), then the accurate phase distribution can be obtained by ( ( )) ( ) 𝐼̃𝑑𝑖𝑓 − 𝑎1 − 𝑎2 𝜃 𝜑𝑘 = tan−1 cot ⋅ 𝑘𝑠𝑢𝑚 ( (32) ) − 𝜃∕2. 2 𝐼𝑘 − 𝑎1 + 𝑎2

(20)

𝑗=1

𝑗=1

(29)

can be calculated by (25)

We can also extract √ the background intensity 𝑎1 , and the modulation amplitude 𝑏1 = 𝜂 2 1 + 𝜉 2 1 from Eq. (21). 3

Y. Zhang

Optics Communications 460 (2020) 125174

the number of unknowns is more than that of equations, that is why the accuracy of two-step PSAs is not very high, for FLSA, the background intensity and modulation amplitude of different interferograms can be calculated, it can improve the accuracy of general two-step PSAs. In the following simulation, we will discuss that how the non-uniform and variable background intensity and modulation amplitude affect the accuracy.

Table 1 RMS phase errors of EMNSD, EMNSD&FLSA, GS and EVI in three different situations (rad).

3. Simulation

Table 2 The computational time of EMNSD, EMNSD&FLSA, GS and EVI in three different situations (s).

To evaluate the performance of the proposed method, we perform several simulations, different background intensity and modulation amplitude distributions, different phase shifts, different levels of noises, different fringe numbers, different kinds of fringes and different chosen samples will be studied in the following, moreover, we compare EMNSD, EMNSD&FLSA with the well-evaluated two-step PSAs based on Gram–Schmidt orthonormalization method (GS) and extreme values of the interference (EVI), all computations are performed with the CPU of Intel(R) Core(TM) i7-6700 and 8 GB memory, and we use the Matlab software for coding. Firstly, we perform EMNSD, EMNSD&FLSA, GS and EVI to process two phase shifted interferograms with the circular fringes in three ( ) different situations, the tested phase is set as 𝜑 = 4𝜋 𝑥2 + 𝑦2 . In situation 1, the background intensity and modulation amplitude distributions of the two interferograms are uniform and invariable, 𝑎1 = 𝑎2 = 1, 𝑏1 = 𝑏2 = 1. In situation 2, the background intensity and [ ( )] modulation amplitude are set as 𝑎𝑖 (𝑥, 𝑦) = 𝑁𝑎 exp −0.02 𝑥2 + 𝑦2 [ ( 2 )] and 𝑏𝑖 (𝑥, 𝑦) = 𝑁𝑏 exp −0.02 𝑥 + 𝑦2 respectively, 𝑁𝑎 of the 1st and 2nd interferograms are set as 1 and 0.95, 𝑁𝑏 of the 1st and 2nd interferograms are set as 0.95 and 0.9. In situation 3, we add noise with SNR of 20 dB which is generated by the function awgn in Matlab to situation 2, it is the most complex situation which is more similar to the actual case, the effectiveness of the proposed method will be verified in this situation. With the above parameters setting, two simulated phase shifted interferograms with 401 × 401 pixels and the phase shift of 2 rad in situation 3 are generated, as shown in Fig. 1(a) and (b), the phase distribution is shown in Fig. 1(c). Moreover, for EMNSD&FLSA, only 81 × 81 pixels are uniformly selected to take part in the iterative process, it will highly save time, and the predefined converging threshold of iteration is 10−4 rad. Note that, in the following, the PV value of the phase distribution means the difference between the maximum and minimum phase, and the RMS value means the RootMean-Square of phase distribution, moreover, the RMS phase error means the Root-Mean-Square of the phase error. Fig. 2 shows the phase distributions and phase error distributions of four different methods in three different situations, N represents the index of the situation, the 1st to 3rd lines show the phase distributions calculated by EMNSD, EMNSD&FLSA, GS and EVI in three different situations, and the phase error distributions are displayed in the 4th to 6th lines. We can see that the phase distributions are similar to Fig. 1(c), that is to say, these four methods are effective in the different situations. However, the phase error distributions are different for the different situations and methods, from the phase error distributions in Fig. 2 and Table 1 which is the RMS phase errors of EMNSD, EMNSD&FLSA, GS and EVI in three different situations, we can see that, for the different situations, the phase errors of GS and EVI are larger than that of EMNSD and EMNSD&FLSA, and they are similar for the different situations because the filtering error is the main error for GS and EVI, the phase error of EMNSD&FLSA is less than that of EMNSD because of the iterative process of FLSA, also, the original phase shifted interferograms with the background intensity take part in the calculation, moreover, for situation 1, the phase error of EMNSD&FLSA is close to zero, that is to say, it can obtain absolutely accurate result when the experimental condition is perfect, but EMNSD cannot do this since there are some approximations, such as Eqs. (4) and (9), for EMNSD and EMNSD&FLSA in situation 3, the main error is noise,

N

EMNSD

EMNSD&FLSA

GS

EVI

1 2 3

0.0261 0.0426 0.1298

6.5926 × 10−5 0.0255 0.1251

0.1704 0.1672 0.1818

0.1853 0.1673 0.1822

N

EMNSD

EMNSD&FLSA

GS

EVI

1 2 3

0.65 0.64 0.65

0.74 0.75 0.75

2.93 2.93 2.9

3.04 3.03 2.83

so the phase error distributions are similar to the noise distribution, and we can conclude that, for EMNSD and EMNSD&FLSA, the more complex the situation, the larger the phase error is. We also compare the computational time of these four different methods, as shown in Table 2, we can conclude that, for each PSA, the computational time in the different situations are similar because of the same pixels, and EMNSD&FLSA costs a little longer time than EMNSD due to the iterative process, however, the iterative process costs less time which is less than 0.1 s because only a limited of pixels are chosen to take part in the iteration, and the computational time of EMNSD and EMNSD&FLSA is further less than that of GS and EVI because the filtering process costs more time. Furthermore, it is important to note that the consuming time of the phase wrapping has been added to the computational time of every PSA, for the single phase wrapping process, it costs 0.53 s for 401 × 401 pixels, so we can know that, the EMNSD&FLSA which excludes the phase wrapping process only costs less than 0.3 s, it is very fast. To analyze the effects of different phase shifts for four different methods, we calculate the RMS phase errors of EMNSD, EMNSD&FLSA, GS and EVI with different phase shifts in situation 3 because situation 3 is most complex and close to the practical situation, and to study the effects of different levels of noises to the range of phase shift, we analyze two different levels of noises with SNR of 20 dB and 30 dB, the results are shown in Fig. 3. For GS and EVI, the RMS phase errors are irrelevant to the different phase shifts since the main error is the filtering error, the effect of the filtering error is more larger than that of the different phase shifts, and we find a strange phenomenon that the range of phase shift with 20 dB noise which is between 0.2 rad and 3.0 rad is wider than that with 30 dB noise which is between 0.5 rad and 2.6 rad for GS and EVI because the filtering process not only filters the background intensity, but also filters the noise, and the larger the noise, the better the filtering effect is, and the RMS phase errors for different levels of noises for GS and EVI are similar because most of the noise has been filtered, the main part of phase error is the filtering error. However, the RMS phase errors of EMNSD and EMNSD&FLSA are relevant to the different phase shifts and levels of noises, we found that the farther away the phase shift from 0 rad and 𝜋 rad, the smaller the RMS phase error is for different levels of noises, and the range of phase shift with 20 dB noise which is between 0.5 rad and 2.7 rad is smaller than that with 30 dB noise which is 0.2 rad and 3.0 rad, that is to say, the less complex the situation, the larger the working range of phase shift is for EMNSD and EMNSD&FLSA, and the RMS phase error of EMNSD and EMNSD&FLSA with 20 dB noise is larger than that with 30 dB noise since they use the original phase shifted interferograms with noise to obtain the phase, and they cannot suppress the noise. Finally, we find that the accuracy of EMNSD&FLSA is a little higher than that of EMNSD for different levels of noises and phase shifts, for 20 dB noise, when the phase shift is between 0.7 rad and 2.3 rad, the RMS phase error of EMNSD and EMNSD&FLSA is less than that of GS and EVI, the phase error of EMNSD and EMNSD&FLSA is larger than that of GS and EVI for other range of phase shift, and for 30 dB, the 4

Y. Zhang

Optics Communications 460 (2020) 125174

Fig. 1. Simulated phase shifted interferograms with the circular fringes and phase distribution. (a) and (b) the first and second interferograms, (c) the theoretical phase distribution (PV = 25.1327 rad, RMS = 5.3249 rad).

Fig. 2. The phase distributions and phase error distributions of EMNSD, EMNSD&FLSA, GS and EVI in three different situations.

RMS phase error of EMNSD and EMNSD&FLSA is less than that of GS and EVI for all phase shifts. Hence, if the high accuracy is required, it would be best to suppress the noise before using EMNSD&FLSA, and the phase shift would be best to far away from 0 rad and 𝜋 rad. We know that two-step PSA is easily influenced by the noise, hence we estimate the noise effect to four different methods in the following. The SNR of noise is set from 20 dB to 80 dB, other parameters are same as the circular fringes in situation 3. We plot the RMS phase errors of the different levels of noises for four different methods, as shown in Fig. 4. Since the large filtering error of GS and EVI, the RMS phase errors of GS and EVI are larger than that of EMNSD and EMNSD&FLSA for any level of noise, and the RMS phase errors of GS and EVI are similar, moreover, when the SNR of noise is less than 50 dB, the larger the noise, the larger the RMS phase error of GS and EVI is, when the SNR of noise is more than 50 dB, the RMS phase error

is ruleless. For EMNSD and EMNSD&FLSA, the larger the noise, the larger the RMS phase error is for any level of noise, also, the phase errors are relatively small and stable when the SNR of noise is more than 50 dB, in this situation, the main phase error is caused by the non-uniform and variable background intensity, modulation amplitude and the intrinsic error of the algorithm, in addition, EMNSD&FLSA has more outstanding performance in regard to the different levels of noises than EMNSD. In Section 2.1, in order to obtain the accurate phase distribution, we assume that there is more than one fringe in the interferogram, in the following, we vary the number of the phase shifted fringe patterns to obtain the range of the fringe numbers using EMNSD&FLSA, provided ( ) that the tested phase distribution 𝜑 = 𝑘𝜋 𝑥2 + 𝑦2 , which k represents the number of the phase shifted fringe patterns in one interferogram, other parameters are same as the circular fringes in situation 3. As can 5

Y. Zhang

Optics Communications 460 (2020) 125174

Fig. 3. The RMS phase errors of EMNSD, EMNSD&FLSA, GS and EVI with different phase shifts in situation 3 with 20 and 30 dB noise.

error is decreasing with the increase of fringe numbers, that is to say, the approximation error from Eqs. (4) and (9) is decreasing with the increase of fringe numbers. When the numbers of patterns are more than 2, the RMS phase errors are relatively small and similar, in this case, the approximation error is nearly stable, hence, we can conclude that the numbers of patterns are best to be more than 2 if high accuracy is requested. To verify the robustness of the proposed method, we also simulate the straight and complex fringes, the comparisons of EMNSD, EMNSD&FLSA, GS and EVI are also performed in the following. For the straight fringes, the theoretical phase is set as 𝜑 = 4𝜋𝑥, and for the complex fringes, the phase is set as 𝜑 = 4𝜋𝑥 + 4𝑝𝑒𝑎𝑘𝑠 (401), other parameters are same as the circular fringes in situation 3. Fig. 5 shows one of the simulated interferograms with the straight fringes and phase distribution, and Fig. 6 represents the simulated results of the straight fringes using EMNSD, EMNSD&FLSA, GS and EVI. The interferogram with the complex fringes which are asymmetrical is shown in Fig. 7(a), the complex phase distribution is drawn in Fig. 7(b), and the simulated results are shown in Fig. 8. For the straight fringes, the RMS phase errors of EMNSD, EMNSD&FLSA, GS and EVI are respectively 0.1259 rad, 0.1238 rad, 0.1832 rad and 0.1901 rad, the computational time of them are respectively 0.65 s, 0.75 s and 2.91, 2.81 s. For the complex fringes, the RMS phase errors of EMNSD, EMNSD&FLSA, GS and EVI are 0.1277 rad, 0.1242 rad, 0.2066 rad and 0.2292 rad, the computational time of them are respectively 0.65 s, 0.75 s, 2.91 s and 2.81 s. From the above simulations, we can get the same conclusion as the above circular fringes, hence we can conclude that EMNSD, EMNSD&FLSA, GS and EVI are all effective for the circular, straight and complex fringes, and EMNSD&FLSA can obtain relatively accurate phase distribution and cost relatively less time simultaneously. Finally, we perform the proposed method with different chosen samples to compare the accuracy and computational time, note that, the samples are chosen with the regular intervals, the computational time and RMS phase errors of the circular, straight and complex fringes with different chosen samples are shown in Table 4, where T and P represent computational time and RMS phase error, and Cir, Str

Fig. 4. RMS phase errors of different methods with different levels of noises.

Fig. 5. Simulated interferogram with the straight fringes and phase distribution. (a) One of the simulated interferograms, (b) the theoretical phase distribution (PV = 25.1327 rad, RMS = 7.2733 rad).

be seen from Table 3, when the fringe number is less than 0.8, the RMS phase error is relative larger, and the ratio of RMS phase to RMS phase error is less than 10 (In general, the ratio of RMS phase to RMS phase error is more than 10 in the accurate measurement). For the range of fringe numbers between 0.9 and 1.2, the RMS phase errors are unstable. When the range of fringe numbers are between 1.3 and 2.0, the ratio is increasing with the increase of fringe numbers, and the RMS phase

Fig. 6. Simulated results of the straight fringes using EMNSD, EMNSD&FLSA, GS and EVI.

6

Y. Zhang

Optics Communications 460 (2020) 125174 Table 3 RMS phase and RMS phase errors with different numbers of the phase shifted fringe patterns using EMNSD&FLSA. Patterns

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

RMS phase (rad) RMS phase error (rad)

0.6656 0.3551

0.7987 0.2322

0.9319 0.1440

1.0650 0.1322

1.1981 0.1287

1.3312 0.1318

1.4643 0.1290

1.5975 0.1305

Patterns

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

RMS phase (rad) RMS phase error (rad)

1.7306 0.1291

1.8637 0.1289

1.9968 0.1284

2.1300 0.1272

2.2631 0.1268

2.3962 0.1267

2.5293 0.1265

2.6625 0.1263

Patterns

3.0

4.0

5.0

15

25

35

45

RMS phase (rad) RMS phase error (rad)

3.9937 0.1255

5.3249 0.1251

6.6561 0.1247

19.9684 0.1248

33.2807 0.1247

46.5929 0.1247

59.9052 0.1248

Table 4 Computational time and RMS phase errors of different fringes with different chosen samples. Samples

6 × 6

11 × 11

21 × 21

41 × 41

81 × 81

101 × 101

201 × 201

401 × 401

Cir

T (s) P (rad)

0.66 0.1545

0.66 0.1257

0.68 0.1254

0.68 0.1257

0.75 0.1251

0.78 0.1250

1.15 0.1249

3.82 0.1249

Str

T (s) P (rad)

0.68 0.1263

0.68 0.1244

0.68 0.1238

0.68 0.1239

0.75 0.1238

0.76 0.1238

1.12 0.1238

3.44 0.1238

Com

T (s) P (rad)

0.67 0.1259

0.67 0.1248

0.68 0.1242

0.69 0.1243

0.76 0.1242

0.80 0.1242

1.16 0.1242

3.01 0.1242

samples with 401 × 401, the best chosen samples are 81 × 81, the conclusion is different from Ref. [26] because they are two different methods, EMNSD&FLSA is a two-step PSA, FIA is a two-step PSA. Moreover, when the predefined converging threshold of iteration is 10−4 rad, the number of iterations is 9, and the iterative process only costs less than 0.1 s with 81 × 81 chosen pixels, hence, EMNSD&FLSA can obtain relatively high accuracy with less time even though it is an iterative algorithm. 4. Experiments Fig. 7. Simulated interferogram with the complex fringes and phase distribution. (a) One of the simulated interferograms, (b) the theoretical phase distribution (PV = 57.8814 rad, RMS = 11.0403 rad).

In the following, three groups of experiments are performed to do the phase retrieval by EMNSD, EMNSD&FLSA, GS and EVI. We respectively capture the circular, straight and complex fringes. Standard 4-step PSA is a highly accurate and high-speed PSA, it needs four phase-shifted interferograms with the phase shifts 0, 𝜋/2, 𝜋 and 3𝜋/2, and its accuracy is easily affected by the phase shift error, in our experiment, the phase shift error will be very small because four phaseshifted interferograms with the phase shifts 0, 𝜋/2, 𝜋 and 3𝜋/2 are extracted from a single image snapshotted by the polarization camera, hence the highly accurate phase extracted by standard 4-step PSA can be set as the reference phase. For the 1st experiment, the phase shifted interferograms with the circular fringes (401 × 401 pixels) are collected, moreover, for EMNSD&FLSA, only 81 × 81 pixels are uniformly selected to take part in the iterative process as the simulation, one of the interferograms is shown in Fig. 9(a), the reference phase

and Com represent the circular, straight and complex fringes. For the different fringes, the computational time is increasing with the increase of chosen samples, and when the chosen samples are less than 101 × 101, the computational time is increasing slowly with the increase of chosen samples, but when the chosen samples are more than 101 × 101, the computational time is increasing quickly with the increase of chosen samples. Moreover, when the chosen samples are 6 × 6, the RMS phase error is largest for all kinds of fringes because the chosen samples are too few, and the RMS phase errors are relatively stable when the chosen pixels are more than 81 × 81. Hence, to obtain high accuracy and cost less computational time simultaneously, for the

Fig. 8. Simulated results of the complex fringes using EMNSD, EMNSD&FLSA, GS and EVI.

7

Y. Zhang

Optics Communications 460 (2020) 125174

Fig. 9. Experimental results of the circular fringes. (a) One of the phase shifted interferograms, (b) the reference phase distribution extracted by 4-step PSA (PV = 43.4713 rad, RMS = 8.9328 rad), (c), (d), (e) and (f) the phase distributions extracted by EMNSD (PV = 43.8433 rad, RMS = 8.9343 rad), EMNSD&FLSA (PV = 43.8187 rad, RMS = 8.9328 rad), GS (PV = 43.9180 rad, RMS = 8.8639 rad) and EVI (PV = 43.9071 rad, RMS = 8.8641 rad), (g), (h), (i) and (j) the differences between the reference and phase distributions extracted by EMNSD, EMNSD&FLSA, GS and EVI.

Fig. 10. Experimental results of the straight fringes. (a) One of the phase shifted interferograms, (b) the reference phase distribution extracted by 4-step PSA (PV = 22.2026 rad, RMS = 6.0503 rad), (c), (d), (e) and (f) the phase distributions extracted by EMNSD (PV = 22.3561 rad, RMS = 6.0422 rad), EMNSD & FLSA (PV = 22.2407 rad, RMS = 6.0458 rad), GS (PV = 24.0806 rad, RMS = 6.0492 rad) and EVI (PV = 24.1965 rad, RMS = 6.0473 rad), (g), (h), (i) and (j) the differences between the reference and phase distributions extracted by EMNSD, EMNSD&FLSA, GS and EVI.

8

Y. Zhang

Optics Communications 460 (2020) 125174

Fig. 11. Experimental results of the complex fringes. (a) One of the phase shifted interferograms, (b) the reference phase distribution extracted by 4-step PSA (PV = 37.9868 rad, RMS = 8.0988 rad), (c), (d), (e) and (f) the phase distributions extracted by EMNSD (PV = 38.6386 rad, RMS = 8.0858 rad), EMNSD&FLSA (PV = 38.6187 rad, RMS = 8.0850 rad), GS (PV = 38.0239 rad, RMS = 8.0315 rad) and EVI (PV = 37.9909 rad, RMS = 8.0325 rad), (g), (h), (i) and (j) the differences between the reference and phase distributions extracted by EMNSD, EMNSD&FLSA, GS and EVI.

phase distributions extracted by EMNSD, EMNSD&FLSA, GS, EVI and the reference phase distribution are 0.0516 rad, 0.0470 rad, 0.1433 rad and 0.1475 rad, the computational time are 0.45 s, 0.53 s, 1.23 s, and 1.36 s. From the above results, we can see that, all four methods are effective for the different fringes in the practical experiments, and EMNSD&FLSA has the best performance among the four methods for all kinds of fringes. Through the above experiments, we can get the same conclusions as the simulations, we verify that, for the circular, straight and complex fringes, the proposed EMNSD&FLSA without pre-filtering can obtain relatively high measurement accuracy, and costs relatively less time by only two randomly phase shifted interferograms.

distribution is plotted in Figs. 9(b), and 8(c), (d), (e) and (f) show the extracted phase distribution using EMNSD, EMNSD&FLSA, GS and EVI which are similar to Fig. 9(b), it means that EMNSD, EMNSD&FLSA, GS and EVI can get the accurate results, and Fig. 9(g), (h), (i) and (j) present the differences between the reference and phase distributions extracted by EMNSD, EMNSD&FLSA, GS and EVI, the RMS values are respectively 0.0861 rad, 0.0806 rad, 0.1452 rad and 0.1584 rad, further indicating that the accuracy of EMNSD&FLSA is higher than that of EMNSD, GS and EVI have the lowest accuracy because of the filtering process, moreover, the computational time of EMNSD, EMNSD&FLSA, GS and EVI are 1.08 s, 1.34 s, 3.56 s and 3.76 s, the computational time of the experimental data is longer than that of the simulated data, although the pixels of the interferograms are both 401 × 401, the experimental condition may worse than the simulated condition, hence the computational time will be affected, furthermore, as the simulated results, EMNSD costs least time, and the computational time of EMNSD&FLSA is only a little longer than that of EMNSD, FLSA costs less time which is less than 0.1 s. Then, the 2nd and 3rd experiments with the straight and complex fringes are performed, and we use the asymmetrical fringes as the complex fringes, the pixels of the interferograms with the straight fringes are also 401 × 401, and the pixels of the interferograms with the complex fringes are 201 × 201, other conditions are same as the above circular fringes. Figs. 10 and 11 show the results of the straight and complex fringes, for the straight fringes, the RMS values of the difference between the phase distributions extracted by EMNSD, EMNSD&FLSA, GS, EVI and the reference phase distribution are respectively 0.0697 rad, 0.0687 rad, 0.2326 rad and 0.2570 rad, the computational time are 0.73 s, 0.88 s, 3.42 s, and 3.66 s, and for the complex fringes, the RMS values of the difference between the

5. Conclusion In this paper, a two-step PSA based on the Euclidean matrix norm of sum and difference map and fast least-squares algorithm is presented. Firstly, the sum and difference maps are obtained, and the background intensity is filtered by subtracting the mean of sum map, after calculating the Euclidean matrix norm of sum and difference maps, the initial phase are extracted, then the original interferograms with the background intensity take part in the calculation of FLSA to obtain the more accurate phase, although FLSA is an iterative algorithm, it also costs less time as the non-iterative algorithm since only a limited of chosen pixels take part in the iterative process. We have compared EMNSD, EMNSD&FLSA and well-evaluated GS, EVI by the simulated and experimental data. The proposed method can achieve the high accuracy as the iterative method and cost less time as the non-iterative method with no-filtering and only two randomly phase shifted interferograms, it removes the restriction that two-step PSA is difficult to 9

Y. Zhang

Optics Communications 460 (2020) 125174

obtain the high accuracy because most of the two-step PSAs filter the background intensity firstly before performing the phase retrieval, and the filtering process may introduce large error, in addition, it can obtain absolutely accurate phase distribution when the background intensity and modulation amplitude are perfect, and it is suitable for different levels of noises, furthermore, it is roughness for the circular, straight and complex fringes, lastly, if the higher accuracy is requested, the noise is best to be suppressed, and it is best to choose a phase shift which is far away from 0 rad and 𝜋 rad, the numbers of patterns are best to be more than 2. The simulations and experiments demonstrate the effectiveness of the proposed method. In summary, this proposed method is a power tool for the phase retrieval with random phase shift.

[10] J. Vargas, J. Quiroga, T. Belenguer, Phase-shifting interferometry based on principal component analysis, Opt. Lett. 36 (8) (2011) 1326–1328, http://dx. doi.org/10.1364/OL.36.001326. [11] J. Vargas, J. Quiroga, T. Belenguer, Analysis of the principle component algorithm in phase-shifting interferometry, Opt. Lett. 36 (12) (2011) 2215–2217, http://dx.doi.org/10.1364/OL.36.002215. [12] J. Deng, K. Wang, D. Wu, X. Lv, C. Li, J. Hao, J. Qin, W. Chen, Advanced principal component analysis method for phase reconstruction, Opt. Express 23 (9) (2015) 12222–12231, http://dx.doi.org/10.1364/OE.23.012222. [13] J. Xu, W. Jin, L. Chai, Q. Xu, Phase extraction from randomly phase-shifted interferograms by combining principal component analysis and least squares method, Opt. Express 19 (21) (2011) 20483–20492, http://dx.doi.org/10.1364/ OE.19.020483. [14] K. Yatabe, K. Ishikawa, Y. Oikawa, Improving principal component analysis based phase extraction method for phase-shifting interferometry by integrating spatial information, Opt. Express 24 (20) (2016) 22881–22891, http://dx.doi.org/10. 1364/OE.24.022881. [15] K. Yatabe, K. Ishikawa, Y. Oikawa, Simple, flexible and accurate phase retrieval method for generalized phase-shifting interferometry, J. Opt. Soc. Amer. A 34 (1) (2017) 87–96, http://dx.doi.org/10.1364/JOSAA.34.000087. [16] K. Yatabe, K. Ishikawa, Y. Oikawa, Hyper ellipse fitting in subspace method for phase-shifting interferometry: Practical implementation with automatic pixel selection, Opt. Express 25 (23) (2017) 29401–29416, http://dx.doi.org/10.1364/ OE.25.029401. [17] H. Wang, C. Luo, L. Zhong, S. Ma, X. Lu, Phase retrieval approach based on the normalized difference maps induced by three interferograms with unknown phase shifts, Opt. Express 22 (5) (2014) 5147–5154, http://dx.doi.org/10.1364/ OE.22.005147. [18] Y. Xu, Y. Wang, Y. Ji, H. Han, W. Jin, Three-frame generalized phase-shifting interferometry by a Euclidean matrix norm algorithm, Opt. Lasers Eng. 84 (2016) 89–95, http://dx.doi.org/10.1016/j.optlaseng.2016.04.011. [19] C.T. Farrell, M.A. Player, Phase step measurement and variable step algorithms in phase-shifting interferometry, Meas. Sci. Technol. 3 (1992) 953–958, http: //dx.doi.org/10.1088/0957-0233/3/10/003. [20] F. Liu, J. Wang, Y. Wu, F. Wu, M. Trusiak, K. Patorski, Y. Wan, Q. Chen, X. Hou, Simultaneous extraction of phase and phase shift from two interferograms using lissajous figure and ellipse fitting technology with Hilbert-Huang prefiltering, J. Opt.-UK 18 (2016) 105604, http://dx.doi.org/10.1088/2040-8978/18/10/ 105604. [21] J. Vargas, J. Quiroga, C. Sorzano, J. Estrada, J. Carazo, Two-step demodulation based on the Gram–Schmidt orthonormalization method, Opt. Lett. 37 (3) (2012) 443–445, http://dx.doi.org/10.1364/OL.37.000443. [22] J. Deng, H. Wang, F. Zhang, D. Zhang, L. Zhong, X. Lu, Two-step phase demodulation algorithm based on the extreme value of interference, Opt. Lett. 37 (22) (2012) 4669–4671, http://dx.doi.org/10.1364/ol.37.004669. [23] C. Luo, L. Zhong, P. Sun, H. Wang, J. Tian, X. Lu, Two-step demodulation algorithm based on the orthogonality of diamond diagonal vectors, Appl. Phys. B 119 (2015) 387–391, http://dx.doi.org/10.1007/s00340-015-6087-z. [24] W. Niu, L. Zhong, P. Sun, W. Zhang, X. Lu, Two-step phase retrieval algorithm based on the quotient of inner products of phase-shifting interferograms, J. Opt. 17 (2015) 085703, http://dx.doi.org/10.1088/2040-8978/17/8/085703. [25] Z. Cheng, D. Liu, Fast and accurate wavefront reconstruction in two-frame phaseshifting interferometry with unknown phase step, Opt. Lett. 43 (13) (2018) 3033–3036, http://dx.doi.org/10.1364/ol.43.003033. [26] Y. Zhang, X. Tian, R. Liang, Timesaving phase retrieval approach based on difffference map normalization and fast iterative algorithm, Opt. Lasers Eng. 121 (2019) 18–28, http://dx.doi.org/10.1016/j.optlaseng.2019.03.010.

CRediT authorship contribution statement Yu Zhang: Conceptualization, Funding acquisition, Methodology, Writing - original draft, Writing - review & editing. Acknowledgments This work was supported by the National Natural Science Foundation of China (NSFC) (61905039); Jilin Scientific and Technological Development Program (20190701018GH); Education Department of Jilin Province (JJKH20190691KJ) and State Key Laboratory of Applied Optics. References [1] D. Malacara, Optical Shop Testing, third ed., John Wiley & Sons, Inc, 2007, Chap 1-7. [2] J.H. Bruning, D.R. Herriott, J.E. Gallagher, D.P. Rosenfeld, A.D. White, D.J. Brangaccio, Digital wavefront measuring interferometer for testing optical surfaces and lenses, Appl. Opt. 13 (11) (1974) 2693–2703, http://dx.doi.org/10.1364/ AO.13.002693. [3] C. Tian, S. Liu, Two-frame phase-shifting interferometry for testing optical surfaces, Opt. Express 24 (16) (2016) 18695–18708, http://dx.doi.org/10.1364/ OE.24.018695. [4] D. Malacara, Optical Shop Testing, third ed., John Wiley & Sons, Inc, 2007, Chap 14. [5] P.J. de Groot, Vibration in phase-shifting interferometry, J. Opt. Soc. Amer. A 12 (2) (1995) 354–365, http://dx.doi.org/10.1364/JOSAA.12.000354. [6] L.L. Deck, Suppressing phase errors from vibration in phase-shifting interferometry, Appl. Opt. 48 (20) (2009) 3948–3960, http://dx.doi.org/10.1364/AO.48. 003948. [7] Z. Wang, Advanced iterative algorithm for phase extraction of randomly phaseshifted interferograms, Opt. Lett. 29 (14) (2004) 1671–1673, http://dx.doi.org/ 10.1364/OL.29.001671. [8] J. Xu, Q. Xu, L. Chai, Iterative algorithm for phase extraction from interferograms with random and spatially nonuniform phase shifts, Appl. Opt. 47 (3) (2008) 480–485, http://dx.doi.org/10.1364/AO.47.000480. [9] Y. Chen, P. Lin, C. Lee, C. Liang, Iterative phase-shifting algorithm immune to random phase shifts and tilts, Appl. Opt. 52 (14) (2013) 3381–3386, http: //dx.doi.org/10.1364/AO.52.003381.

10