An accurate phase shift extraction algorithm for phase shifting interferometry

Optics Communications 429 (2018) 144–151

Contents lists available at ScienceDirect

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

An accurate phase shift extraction algorithm for phase shifting interferometry Zhongsheng Zhai a,b, *, Li Zhou a , Yanhong Zhang a , Zhengqiong Dong a,b , Xuanze Wang a,b , Qinghua Lv c a b c

School of Mechanics and Engineering, Hubei University of Technology, China HuBei Key Lab of Manufacture Quality Engineering Wuhan, China Hubei Collaborative Innovation Center for High-efficient Utilization of Solar Energy, Hubei University of Technology, China

ARTICLE

INFO

ABSTRACT

Keywords: Interferometry Interference fringes Phase-shift extraction Ellipse fitting

The accuracy of phase shift extraction has a significant influence on measurement results in surface microtopography interferometry. Phase shifting errors are mainly caused by nonlinearity of the employed phase shifter, environmental turbulence, camera imperfection and so on. In this paper, a general algorithm based on Lissajous figures and ellipse fitting is proposed for extracting the phase distribution from a set of phaseshifting interferograms with random noise. Two sets of pixels with 𝜋∕2 phase difference in all the investigated interferograms are selected and used for ellipse fitting. Both numerical simulations and optical experiments have proven the validity, rapidity, and accuracy of the proposed method. Experiments show that the proposed method is a general phase extraction method, which can work for straight fringe patterns, circle fringes patterns and other anomalous features.

1. Introduction

In 1982, Morgan developed a least-squares iterative algorithm that estimates phases and their perturbation caused by linear time-dependent drifts [19]. In 1991, Okada et al. proposed a least-square-based iterative algorithm to solve a set of approximate linear equations iteratively, which allows the phase shift amounts and phase distributions to be determined simultaneously [20]. In 2007, Wang et al. proposed an advanced iterative algorithm (AIA) which can extract both initial phase distribution and phase shift amounts using three randomly shifted interferograms [21], overcoming the limitation found in conventional iterative algorithms that the number of frames must be at least four. In simulation, the phase extraction errors of the AIA algorithm with three frames are less than 0.0152 rad [15,21]. In 2010, Q. Kemao et al. proposed a windowed Fourier ridges and least squares fitting (WFRLSF) [22], and presented the phase shift errors of the algorithms: the AIA, the WFRLSF, the windowed Fourier transform (WFF) + AIA + WFF, and the WFF + WFRLSF + WFF. In simulation the phase shift errors are less than 0.06 rad, and as noise level is increased the phase shift error ranged from 0.026 to 0.29 rad. Recently, R. Zhu et al. extracted measurement phases from two phase-shifting fringe patterns using the spatial-temporal fringes method [23]. The surface error using this method is 1.8 × 10−3 𝜆. Compared with the iterative algorithms above, the approaches based on non-iterative solutions are intended to find the optimal results in less

High precision surface topography measurement has many applications in areas such as integrated circuits and MEMS [1–5]. Phase-shifting interferometry (PSI) is one of the most widely used techniques in surface topography measurement as it is non-contact, non-destructive and has high accuracy [6–10]. To conduct a complete phase-shift measurement, a piezoelectric ceramic transformer (PZT) is generally used as the phase shifter. However, due to the hysteresis nonlinearity of the PZT and the instability of the environment (including ambient vibration and air shock), random sampling errors are unavoidable in phase shift measurement. The accurate extraction of phase shifts is a challenging task and has significant influence on measurement results. To reduce the errors in the phase-shift extraction, many researchers have tried changing the configuration of interferometers or improving phase-shifting algorithms. For example, simultaneously phase shifting interferometry (SPSI) can effectively avoid environmental vibrations by collecting interferograms instantaneously [11,12], but the limitation of this method is the complicated hardware configuration. In fact, there has been more research focusing on modification of phase-shifting extraction algorithms. All the existing algorithms can be divided into two classes: iterative algorithms [13–15] and non-iterative algorithms [16–18]. *

Corresponding author at: School of Mechanics and Engineering, Hubei University of Technology, China. E-mail address: [email protected] (Z. Zhai).

https://doi.org/10.1016/j.optcom.2018.08.005 Received 26 March 2018; Received in revised form 26 July 2018; Accepted 4 August 2018 Available online 7 August 2018 0030-4018/© 2018 Elsevier B.V. All rights reserved.

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

time. Lissajous ellipse fitting, an example of a non-iterative approach, has proved to be an outstanding algorithm to extract phase shifts with high accuracy and efficiency [24]. In 1994, Farrell and Player used a pair of different pixels in fringe field (inter-pixel) to create Lissajous figures, from which they then calculated phase shift amounts, intensity bias and intensity modulation at each pixel using ellipse fitting, under the condition of both unequal and unknown phase steps [25]. The experimental result of wavefront reconstruction shows that the accuracy of the proposed algorithm is 0.032 rad. In addition, they mentioned that if the phase difference between the pixel pairs is close to ±𝑛𝜋, Bookstein’s algorithm would fail. In 2015, Fengwei Liu et al. proposed to correct the dynamic random phase shift errors by transforming the Lissajous ellipse to a unit circle (ETC) [26]. They deduced that the phase extraction error can be compensated with the unit circle and that the accuracy of correcting the phase extraction error is mainly dependent on the parameters of the ellipse. Experimental results show that the ETC method has similar precision comparable to AIA. Both Refs. [25] and [26] mentioned that when the phase difference between the pair of pixels equals 𝜋∕2 the extracted phase will be most accurate. However, a way of choosing the pixel pairs is lacking. Moreover, the accuracy and reliability of the phase extraction will be appreciably affected by random phase shifting noise if only a single pair of pixels is used for Lissajous ellipse generation. By using a series of pixel pairs with 𝜋∕2 phase difference, we introduce a high-precision phase extraction method based on a leastsquares algorithm and general ellipse fitting. The proposed method chooses the same region from all interferograms, and the average intensity of each selected region is calculated. Then two index numbers of interferograms sequence are obtained by seeking two points with 𝜋∕2 phase difference within the average intensity array. Following this, two groups of pixels with 𝜋∕2 phase difference in all the investigated interferograms are selected and used for ellipse fitting. The proposed method is not restricted in 3-interferograms, and it can effectively suppress random phase-shifting errors.

Fig. 1. The tendency of intensity curves with 35 fringes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. Select the pixels with the same gray-scale change tendency. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We seek two values in the array 𝐺i with 𝜋∕2 phase difference. According to the law that the phase difference between the maximum value and the middle value (amplitude = 0) is 𝑘𝜋 + 𝜋∕2, we can obtain the sequence numbers 𝑚 and 𝑛 from the following equation:

2. Theoretical analysis 2.1. Seeking two interferograms with 𝜋∕2 phase difference

⎧ ⎪𝑚 = 𝑖 ⎨ ⎪𝑛 = 𝑖 ⎩

The phase-shift fringe pattern generated by phase-shifting can be expressed as 𝐼𝑖 (𝑥, 𝑦) = 𝐴𝑖 (𝑥, 𝑦) + 𝐵𝑖 (𝑥, 𝑦) cos[𝜑(𝑥, 𝑦) + 𝜃𝑖 ] + 𝑁𝑖 (𝑥, 𝑦)

(1)

𝑎2 𝑏2 ∑ ∑ 1 𝐼 (𝑥 , 𝑦 ) (𝑎2 − 𝑎1 ) ∗ (𝑏2 − 𝑏1 ) 𝑘=𝑎 𝑙=𝑏 𝑖 𝑘 𝑙 1

(3)

The phase difference between the 𝑚th and the 𝑛th interferograms is close to 𝑘𝜋 + 𝜋∕2 + 𝜀, where 𝜀 is a deviation close to zero. The corresponding intensities 𝐼m and 𝐼n can be expressed as: { 𝐼𝑚 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜑(𝑥, 𝑦) + 𝜃𝑚 ] (4) 𝐼𝑛 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜑(𝑥, 𝑦) + 𝜃𝑚 + 𝑘𝜋 + 𝜋∕2 + 𝜀]

where 𝑖 denotes the sequence number of phase-shifting interferogram, (𝑥, 𝑦) represents the coordinate of an arbitrary pixel, 𝐼𝑖 (𝑥, 𝑦) is the intensity at pixel location (𝑥, 𝑦), 𝐴𝑖 and 𝐵𝑖 represent the background intensity and the modulation amplitude respectively, 𝜑(𝑥, 𝑦) is the initial phase, 𝑁𝑖 (𝑥, 𝑦) denotes the random noise, and 𝜃𝑖 describes the phase shift of the 𝑖th interferogram. 𝜃𝑖 is the main phase shift parameter that we need to obtain in each interference sequence. With ellipse fitting, the accuracy of extracted phase shift 𝜃𝑖 can be improved by using two sets of signals that are reliable with phase difference 𝜋∕2 [25]. The set of points can be selected in two interferograms that have a feature with 𝜋∕2 phase difference. Now we explain how to find these two interference patterns. Firstly, a center area (here, a rectangle with a length of 𝑎2 − 𝑎1 and a width of 𝑏2 − 𝑏1 as an example) of each interferogram is chosen to calculate the average intensity of all pixels in accordance with the expression. 𝐺𝑖 =

| 𝐺 + 𝐺min || 𝑖𝑓 ||𝐺𝑖 − max | = 𝑀𝑖𝑛 2 | | 𝑖𝑓 𝐺𝑖 = 𝐺max

where 𝜃𝑚 is the 𝑚th phase shift value. 2.2. Seeking the pixels with same phase change tendency in mth and nth interferograms Due to the unavoidable random noise, the Lissajous fitting accuracy will be unreliable if just a single pixel is used in each interferogram. Therefore, we propose to replace a single pixel by the average value of a group of pixels selected from the interferograms 𝐼𝑚 and 𝐼𝑛 in a certain way. Specifically, we take all pixels with the same grayscale change tendency within a certain range, such as (−𝜋 < phase < 0). For example, in Fig. 2, 𝑃1 , 𝑃2 and 𝑃3 represent the intensity distributions of three points ((𝑥p1 , 𝑦p1 ), (𝑥p2 , 𝑦p2 ), (𝑥p3 , 𝑦p3 )) in different interferograms. The red pixels ((𝑥p2 , 𝑦p2 ), (𝑥p3 , 𝑦p3 )) in the 𝑚th interferogram will be selected for Lissajous fitting because at each of these points the intensity is increasing, however the pixel (𝑥p1 , 𝑦p1 ) is not selected as its intensity is decreasing. Similarly for the 𝑛th interferogram, suitable pixels have been marked in green in Fig. 2.

(2)

1

where 𝐺i represents the average gray value of all the pixels in the selected region. The average value 𝐺i follows a sine curve, as show in Fig. 1, in which the maximum and minimum values of 𝐺i , i.e., 𝐺max and 𝐺min , are marked with red points. 145

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

Fig. 3. The pixels with same phase condition mapping to 𝑖th interferograms.

In particular, for the interferogram 𝐼𝑚 , we find the pixels which meet the phase condition (−𝜋 < phase < 0) by comparing the intensity values of each pixel between 𝐼𝑚 and 𝐼𝑚−1. If the intensity in location (𝑥, 𝑦) satisfies 𝐼𝑚 (𝑥, 𝑦) >𝐼𝑚−1 (𝑥, 𝑦), then the corresponding location coordinate will be saved into the array ArrayP1 , as shown in Eq. (5). Similarly, the pixel coordinates of interferogram 𝐼𝑛 are added to the array ArrayP2 when their intensities meet the condition 𝐼𝑛 (𝑥, 𝑦) >𝐼𝑛−1 (𝑥, 𝑦). { 𝐴𝑟𝑟𝑎𝑦𝑃1 [𝑝].𝑥 = 𝑥 𝑖𝑓 𝐼𝑚 (𝑥𝑝 , 𝑦𝑝 ) > 𝐼𝑚−1 (𝑥𝑝 , 𝑦𝑝 ) (5) 𝐴𝑟𝑟𝑎𝑦𝑃1 [𝑝].𝑦 = 𝑦 { 𝐴𝑟𝑟𝑎𝑦𝑃2 [𝑞].𝑥 = 𝑥 𝑖𝑓 𝐼𝑛 (𝑥𝑞 , 𝑦𝑞 ) > 𝐼𝑛−1 (𝑥𝑞 , 𝑦𝑞 ) (6) 𝐴𝑟𝑟𝑎𝑦𝑃2 [𝑞].𝑦 = 𝑦

2.3. A general ellipse fitting algorithm based on least-squares According to Eq. (8), we assume that: 𝜙𝑚 = 𝜑(𝑥1 , 𝑦1 ) + 𝜃𝑚 and 𝜑12 = 𝜑(𝑥1 , 𝑦1 ) − 𝜑(𝑥2 , 𝑦2 ) ≈ 𝑘𝜋 + 𝜋∕2, a general expression can be written as: ( ) 𝐼𝑚 (𝑥1 , 𝑦1 ) − 𝐴(𝑥1 , 𝑦1 ) cos 𝜙𝑚 = 𝐵(𝑥 𝑦 ) ( ) 1 ,( 1 ) 𝐼 (𝑥 , 𝑦2 )−𝐴(𝑥2 , 𝑦2 ) cos 𝜙𝑚 cos 𝜑12 − 𝑚 2𝐵(𝑥 2 , 𝑦2 ) sin(𝜙𝑚 ) = sin(𝜑12 )

where the phase value 𝜙𝑚 is wrapped in the range of (−𝜋, 𝜋), and the real phase distribution can be retrieved by using the arctangent function and unwrapping algorithm. On the basis of getting 𝜙𝑚 , the adjacent parameter is reduced to obtain the driven phase shift 𝛿𝑘 . From Eq. (9), a general ellipse form is found to be

The coordinates of the pixels saved in ArrayP1 have the same phase range (−𝜋 < phase < 0), as do the pixels saved in ArrayP2 . After finishing searching all the pixels in the interferograms 𝐼𝑚 and 𝐼𝑛 , the coordinates of pixels in the two arrays ArrayP1 and ArrayP2 will map into the same interferogram, as show in Fig. 3. As there is a 𝑘𝜋+𝜋/2 phase difference between 𝐼𝑚 and 𝐼𝑛 , the pixels between ArrayP1 and ArrayP2 also have a 𝑘𝜋+𝜋/2 phase difference when the pixel coordinates are mapped into the same interferogram. If the array size of the ArrayP1 and ArrayP2 are described by 𝑛1 and 𝑛2 respectively, it means that there are 𝑛1 pixels in the 𝐼𝑚 and 𝑛2 pixels in the 𝐼𝑛 which are used to calculate the average. For a general interferogram 𝐼𝑖 , we use the pixels which have the same coordinates in ArrayP1 as in ArrayP2 , to obtain the two intensity sequences shown in Eq. (7). ⎧ ⎪𝐼𝑖 (𝑥1 , 𝑦1 ) = ⎪ ⎨ ⎪𝐼 (𝑥 , 𝑦 ) = ⎪ 𝑖 2 2 ⎩

𝑛1 1 ∑ 𝐼 (𝐴𝑟𝑟𝑎𝑦𝑃1 [𝑝].𝑥, 𝐴𝑟𝑟𝑎𝑦𝑃1 [𝑝].𝑦) 𝑛1 𝑝=1 𝑖 𝑛2 1 ∑ 𝐼 (𝐴𝑟𝑟𝑎𝑦𝑃2 [𝑞].𝑥, 𝐴𝑟𝑟𝑎𝑦𝑃2 [𝑞].𝑦) 𝑛2 𝑞=1 𝑖

(9)

[𝐼𝑚 (𝑥1 , 𝑦1 ) − 𝐴(𝑥1 , 𝑦1 )]2 𝐵 2 (𝑥1 , 𝑦1 ) −2 cos(𝜑12 ) +

[𝐼𝑚 (𝑥1 , 𝑦1 ) − 𝐴(𝑥1 , 𝑦1 )] [𝐼𝑚 (𝑥2 , 𝑦2 ) − 𝐴(𝑥2 , 𝑦2 )] 𝐵(𝑥1 , 𝑦1 )

[𝐼𝑚 (𝑥2 , 𝑦2 ) − 𝐴(𝑥2 , 𝑦2 )]2 𝐵 2 (𝑥2 , 𝑦2 )

(10)

𝐵(𝑥2 , 𝑦2 )

= sin 2(𝜑12 )

and a conic curve can be written as [24,27,28] 𝑎1 𝐼m 2 (𝑥1 , 𝑦1 ) + 2𝑎2 𝐼m (𝑥1 , 𝑦1 )𝐼m (𝑥2 , 𝑦2 ) + 𝑎3 𝐼m 2 (𝑥2 , 𝑦2 ) +2𝑎4 𝐼m (𝑥1 , 𝑦1 ) + 2𝑎5 𝐼m (𝑥2 , 𝑦2 ) = 1

(7)

(11)

Comparing Eqs.(10) and (11), the coefficients in Eq.(11) can be obtained to be

In this way, 𝐼𝑖 (𝑥1 , 𝑦1 ) and 𝐼𝑖 (𝑥2 , 𝑦2 ) are the intensity sequences with smaller random errors in 𝑖th interference, which can be used in the fitting process instead of single pixel in 𝑖th interferogram. Taking the 𝑚th interferogram as an example (set 𝑖 = 𝑚), the equation is expressed as follows: { 𝐼𝑚 (𝑥1 , 𝑦1 ) = 𝐴(𝑥1 , 𝑦1 ) + 𝐵(𝑥1 , 𝑦1 ) cos[𝜑(𝑥1 , 𝑦1 ) + 𝜃𝑚 ] (8) 𝐼𝑚 (𝑥2 , 𝑦2 ) = 𝐴(𝑥2 , 𝑦2 ) + 𝐵(𝑥2 , 𝑦2 ) cos[𝜑(𝑥2 , 𝑦2 ) + 𝜃𝑚 ]

𝑎1 = 𝐵 2 (𝑥2 , 𝑦2 )∕𝑀 𝑎2 = −2𝐵(𝑥1 , 𝑦1 )𝐵(𝑥2 , 𝑦2 ) cos(𝜑12 )∕𝑀 𝑎3 = 𝐵 2 (𝑥1 , 𝑦1 )∕𝑀 𝑎4 = 2[𝐴(𝑥2 , 𝑦2 )𝐵(𝑥1 , 𝑦1 )𝐵(𝑥2 , 𝑦2 ) cos(𝜑12 ) − 𝐴(𝑥1 , 𝑦1 )𝐵 2 (𝑥2 , 𝑦2 )]∕𝑀 𝑎5 = 2[𝐴(𝑥1 , 𝑦1 )𝐵(𝑥1 , 𝑦1 )𝐵(𝑥2 , 𝑦2 ) cos(𝜑12 ) − 𝐴(𝑥2 , 𝑦2 )𝐵 2 (𝑥1 , 𝑦1 )]∕𝑀 𝑀 = 𝐵 2 (𝑥1 , 𝑦1 )𝐵 2 (𝑥2 , 𝑦2 ) sin2 (𝜑12 ) + 2𝐴(𝑥1 , 𝑦1 )𝐴(𝑥2 , 𝑦2 )𝐵(𝑥1 , 𝑦1 )𝐵(𝑥2 , 𝑦2 ) cos(𝜑12 ) −𝐴2 (𝑥1 , 𝑦1 )𝐵 2 (𝑥2 , 𝑦2 ) − 𝐴2 (𝑥2 , 𝑦2 )𝐵 2 (𝑥1 , 𝑦1 )

Here, the phase difference between 𝜑(𝑥1 , 𝑦1 ) and 𝜑(𝑥2 , 𝑦2 ) is approximately 𝜋∕2. If 𝐼𝑚 (𝑥1 , 𝑦1 ) and 𝐼𝑚 (𝑥2 , 𝑦2 ) are the transversal and longitudinal coordinates respectively, the theoretical trajectory of sequence points should be an ellipse. Parameters 𝐴 (𝑥1 , 𝑦1 ), 𝐴 (𝑥2 , 𝑦2 ), 𝐵 (𝑥1 , 𝑦1 ), 𝐵 (𝑥2 , 𝑦2 ) and initial phase difference (𝜑(𝑥1 , 𝑦1 ) − 𝜑(𝑥2 , 𝑦2 )) can be obtained through ellipse fitting.

Therefore, the parameters of in Eq.(8) are solved as 𝜑12 = arccos( 146

𝑎2 ) √ −2 𝑎1 𝑎3

(12)

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

It is easy to calculate the initial phase 𝜑(𝑥, 𝑦) of each pixel, and it is straightforward to solve the problem by using the linear three-parameter minimum sinusoidal method proposed by IEEE [29]. 3. Simulations and experiments Simulations were performed to verify the validity of our proposed method. The interference fringes shown in Fig. 4 were simulated with uniform and random phase shifts respectively. Random phase shifts were created by adding random numbers to the fixed interval 𝜃1 (0.1794 rad). The size of the fringe images simulated with MATLAB was 2000 × 3000 pixels. The extracted uniform and random phase shifts calculated by our proposed method are shown in Tables 1 and 2, respectively. From Fig. 5(a) and (b) it can be seen that their unwrapped phases are related linearly to the sequence of interferograms as the shift is a cumulative results. Tables 1 and 2 also display the error of the calculated phase shifts, equal to the difference between the extracted and actual values. We find that the extracted errors are less than ±0.0004 rad whether the phase shift is uniform or not, while the average phase-shift error of the AIA algorithm is less than 0.0143 rad [15]. Further, experiments on a plane mirror are conducted to demonstrate the effectiveness of the proposed method. The interferometer is carried out with a Mirau interference objective and the experimental setup as shown in Fig. 6. A PZT is employed as the phase-shifting inducer, used to drive the mirror to move in one direction with 18 steps (phase shift > 2𝜋). A white LED is used as the light source. The light passes through a bandpass filter (center wavelength 520 ± 2 nm, full width half max 10 ± 2 nm) to be used as the illumination source of the interferometer. A CCD (1292 × 964 pixels) camera is used to capture the interferograms.

Fig. 4. The fringes generated with (a) uniform phase-shift, and (b) random phase shift.

𝐴(𝑥1 , 𝑦1 ) =

𝐴(𝑥2 , 𝑦2 ) =

2𝑎3 𝑎4 − 𝑎2 𝑎5 𝑎2 2 − 4𝑎1 𝑎3 2𝑎1 𝑎5 − 𝑎2 𝑎4

𝑎2 2 − 4𝑎1 𝑎3 √ 𝑎3 𝑎 𝑎 2 + 𝑎1 𝑎5 2 − 𝑎2 𝑎4 𝑎5 (𝑎5 + 3 4 ) 𝐵(𝑥1 , 𝑦1 ) = 2 2 𝑎2 − 4𝑎1 𝑎3 𝑎2 2 − 4𝑎1 𝑎3 √ 𝑎 𝑎 2 + 𝑎1 𝑎5 2 − 𝑎2 𝑎4 𝑎5 𝑎1 𝐵(𝑥2 , 𝑦2 ) = 2 (𝑎5 + 3 4 ) 2 𝑎2 − 4𝑎1 𝑎3 𝑎2 2 − 4𝑎1 𝑎3

(13)

The parameters a1 , a2 , a3 , a4 and a5 in Eq. (11) can be calculated by least-squares fitting. 𝛿𝑘 is random interval of phase-shift and described in Eq. (14). Because the PZT movement starts from the second picture, the subscript k starts from 2 referring to the first interferogram. 𝜃𝑖 =

𝑖 ∑

𝛿𝑘

(14)

𝑖=2

Fig. 5. The extracted results of (a) uniform phase-shifts and (b) random phase-shifts.

Fig. 6. Experimental setup. 147

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

Fig. 7. Original interferogram and four different selected images from the original interferogram. Table 1 Extracted results from fringes with uniform phase step 0.1794 rad. No.

Real 𝜃𝑖

Extracted

Errors

No.

Real 𝜃𝑖

Extracted

Errors

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

0.1794 0.3588 0.5382 0.7176 0.8970 1.0764 1.2558 1.4352 1.6146 1.7940 1.9734 2.1528 2.3322 2.5116 2.6910 2.8704 3.0498

0.1793 0.3588 0.5381 0.7177 0.8969 1.0766 1.2555 1.4351 1.6150 1.7938 1.9734 2.1526 2.3320 2.5119 2.6910 2.8701 3.0498

−0.0001 0.0000 −0.0001 0.0001 −0.0001 0.0002 −0.0003 −0.0001 0.0004 −0.0002 0.0000 −0.0002 −0.0002 0.0003 0.0000 −0.0003 0.0000

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

3.2292 3.4086 3.5880 3.7674 3.9468 4.1262 4.3056 4.4850 4.6644 4.8438 5.0232 5.2026 5.3820 5.5614 5.7408 5.9202 6.0996

3.2294 3.4088 3.588 3.7672 3.9467 4.1265 4.3055 4.4848 4.6645 4.8436 5.0235 5.2028 5.379 5.5614 5.7407 5.9204 6.0995

0.0002 0.0002 0.0000 −0.0002 −0.0001 0.0003 −0.0001 −0.0002 0.0001 −0.0002 0.0003 0.0002 −0.0003 0.0000 −0.0001 0.0002 −0.0001

Table 2 Extracted results from fringes with random phase shift. No.

Real 𝜃𝑖

Extracted

Errors

No.

Real 𝜃𝑖

Extracted

Errors

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

0.1675 0.3305 0.5009 0.6561 0.8281 0.9962 1.1645 1.3376 1.5032 1.6695 1.8301 1.9888 2.1523 2.3207 2.4850 2.6453 2.8177

0.1676 0.3305 0.5009 0.6563 0.8279 0.9961 1.1641 1.3375 1.5029 1.6694 1.8299 1.9887 2.1520 2.3207 2.4849 2.6451 2.8178

0.0001 0.0000 0.0000 −0.0002 −0.0002 −0.0001 −0.0004 −0.0001 −0.0003 −0.0001 −0.0002 −0.0001 −0.0003 0.0000 −0.0001 −0.0002 0.0001

18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

2.9811 3.1480 3.3150 3.4808 3.6498 3.8155 3.9742 4.1361 4.3043 4.4650 4.6303 4.7926 4.9537 5.1156 5.2806 5.4485 5.6149

2.9811 3.1481 3.3150 3.4808 3.6499 3.8153 3.9741 4.1357 4.3041 4.4648 4.6302 4.7924 4.9537 5.1155 5.2805 5.4484 5.6148

0.0000 0.0001 0.0000 0.0000 0.0001 −0.0002 −0.0001 −0.0004 −0.0002 −0.0002 −0.0001 −0.0002 0.0000 −0.0001 −0.0001 −0.0001 −0.0001

Because the real phase shifts generated by the PZT (open-loop) is unknown, the accuracy of the proposed method only can be discussed in terms of multiple measurement results. According to Eq. (2), selecting the area is a key step to find the two points with 𝑘𝜋 + 𝜋∕2 phase difference. Four different areas (50*50, 100*100, 150*150, 200*200 pixels) represented by A50, A100, A150 and A200 are selected to perform the phase extraction. Fig. 7 shows the first interferogram and four different selected images from it. We see that the straight fringes are polluted by many small particles, which is due to dust particles adhering to the measured object (a mirror). We use these interferograms with high noise to test the effectiveness of the proposed method. Fig. 8 presents the distributions of gray mean value of the 18 interferograms calculated using Eq. (2) corresponding to the selected regions: A50, A100, A150 and A200. The maximum values and the middle values are obtained by Eq. (3) and marked in red in Fig. 8.

It can be observed that the values 𝑛 and 𝑚 vary as the selected region changes. The interference sequence gray mean values are obtained via Eq. (8), and Lissajous figures formed by these values are shown in Fig. 9. The phase difference 𝜑12 of different conditions were obtained as: −1.6943, −1.4628, −1.4628 and −1.6280 rad corresponding to select area A50, A100, A150 and A200. From the phase difference results, we can see that the selected two groups of sequences almost satisfy the orthogonal relation (i.e. phase difference is 𝜋∕2 or 1.5708 rad). The unwrapped phase shifts are calculated by Eq. (8), and Table 3 gives the phase shifts at every step generated by the PZT. In Table 3, A50 , A100 , A150 and A200 denote the phase shifts corresponding to the selected areas A50, A100, A150 and A200 respectively, and their average is defined as Aave . This is considered as the reference phase shift. The phase-shift errors for different selected areas are also presented in 148

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

Fig. 8. Mean intensity distributions of different selected regions in the interferogram sequence (a): A50, (b): A100, (c): A150 and (d): A200. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9. Lissajous figures for (a): A50, (b): A100, (c): A150 and (d): A200.

Table 3. The phase-shift errors of the proposed algorithm vary from 0.0 to 0.0208 rad. The average phase-shift errors eave are very close and are less than 0.0055 rad. If the wavelength is 520 nm, it means that the minimum phase step that can be recognized is 0.46 nm. Another experiment was carried out to further verify the performance of this method in circle interferograms. The reference mirror driven by the PZT has 10 phase steps (phase shift > 2𝜋), and the interferograms are captured by a CCD. Fig. 10 presents the first circle interferogram and the selected areas A50, A100, A150 and A200. With the same steps as those for the straight interferograms, Table 4 presents

the extraction results of the phase shifts. From Table 4, we find that the phase-shift errors vary from 0.0 to 0.0289 rad and the average phaseshift errors eave are less than 0.0083 rad. 4. Discussion When the phase difference between the pair of pixels used for Lissajous ellipse fitting is close to 𝜋∕2, the extracted phase will be most accurate. Random noise on the pixels will affect the accuracy of ellipse fitting, which will influence correctness of the phase-shift extraction. 149

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

Fig. 10. Circle interferogram and four different selected areas from it. Table 3 Phase shift errors (rad) for straight interferograms. No.i

A50

A100

A150

A200

Aave

A50 − Aave

A100 − Aave

A150 − Aave

A200 − Aave

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 eave

0.0000 0.3020 0.5591 0.9294 1.3815 1.8456 2.1755 2.4792 2.8908 3.2477 3.7015 4.2440 4.5032 4.8339 5.2410 5.5805 5.9562 6.3713

0.0000 0.3023 0.5508 0.9293 1.3611 1.8216 2.1771 2.4718 2.8847 3.2314 3.7014 4.2312 4.4913 4.8223 5.2276 5.5785 5.9277 6.3555

0.0000 0.3023 0.5508 0.9293 1.3611 1.8216 2.1771 2.4718 2.8847 3.2314 3.7014 4.2312 4.4913 4.8223 5.2276 5.5785 5.9277 6.3555

0.0000 0.3029 0.5477 0.9509 1.3718 1.8448 2.1755 2.4782 2.8896 3.2468 3.7008 4.2419 4.5016 4.8334 5.2290 5.6069 5.9432 6.3712

0.0000 0.3024 0.5521 0.9347 1.3689 1.8334 2.1763 2.4753 2.8875 3.2393 3.7013 4.2371 4.4969 4.8280 5.2313 5.5861 5.9387 6.3634

0.0000 −0.0004 0.0070 −0.0053 0.0126 0.0122 −0.0008 0.0040 0.0034 0.0084 0.0002 0.0069 0.0064 0.0059 0.0097 −0.0056 0.0175 0.0079 0.0053

0.0000 −0.0001 −0.0013 −0.0054 −0.0078 −0.0118 0.0008 −0.0034 −0.0027 −0.0079 0.0001 −0.0059 −0.0055 −0.0057 −0.0037 −0.0076 −0.0110 −0.0079 −0.0051

0.0000 −0.0001 −0.0013 −0.0054 −0.0078 −0.0118 0.0008 −0.0034 −0.0027 −0.0079 0.0001 −0.0059 −0.0055 −0.0057 −0.0037 −0.0076 −0.0110 −0.0079 −0.0051

0.0000 0.0005 −0.0044 0.0162 0.0029 0.0114 −0.0008 0.0030 0.0022 0.0075 −0.0005 0.0048 0.0048 0.0054 −0.0023 0.0208 0.0045 0.0078 0.0049

Table 4 Phase shift errors (rad) for circle interferogram. No.i

A50

A100

A150

A200

Aave

A50 − Aave

A100 − Aave

A150 − Aave

A200 − Aave

1 2 3 4 5 6 7 8 9 10 eave

0.0000 0.7046 1.5151 2.2953 2.9487 3.6909 4.5457 5.4221 6.1542 6.906

0.0000 0.7046 1.5151 2.2953 2.9487 3.6909 4.5457 5.4221 6.1542 6.906

0.0000 0.6651 1.4769 2.2683 2.9641 3.7023 4.5162 5.41 6.1787 6.9036

0.0000 0.7017 1.5062 2.302 2.9578 3.6996 4.5355 5.4265 6.1668 6.9262

0.0000 0.6940 1.5033 2.2902 2.9548 3.6959 4.5358 5.4202 6.1635 6.9105

0.0000 0.0106 0.0118 0.0051 −0.0061 −0.0050 0.0099 0.0019 −0.0093 −0.0044 0.0014

0.0000 0.0106 0.0118 0.0051 −0.0061 −0.0050 0.0099 0.0019 −0.0093 −0.0044 0.0014

0.0000 −0.0289 −0.0264 −0.0219 0.0093 0.0064 −0.0196 −0.0102 0.0152 −0.0068 −0.0083

0.0000 0.0077 0.0029 0.0118 0.0030 0.0037 −0.0003 0.0063 0.0033 0.0158 0.0054

Using the average intensity of multiple pixels instead of a single pixel can reduce the effect of random noise. From Figs. 7 and 10, it can be noted that the noise level in straight fringe pattern is higher than that in circle fringe pattern. However, the accuracy of phase-shift extraction is very similar. This shows that the proposed method is not sensitive to random noise or dust on the mirror. There is a clear difference between the phase shift errors in simulation and in experiments. In simulation, the illumination source in the whole field is assumed uniform, and the numeral interferograms did not add any noise. This means that using the proposed method we end up with a perfect ellipse, from which we can extract the phase shifts with high accuracy. However, in our experiment, the white LED and bandpass filter will result in the illumination being non-uniform. Since the pixels selected to form the Lissajous ellipse are not illuminated uniformly, ellipse fitting errors are inevitable. Therefore, the accuracy of phaseshift extraction is not high as that in simulations. The proposed method is a non-iterative one. The interferograms were processed on a MacBook Air (1.3 GHz Intel Core i5, 4 GB 1600

MHz DDR3). To obtain the results in Table 3 (18 interferograms at a resolution of 964 × 1292), the algorithm run time was 2.02 s, and for the results in Table 4 the time consumed was 0.39 s (10 interferograms at a resolution of 480 × 640).

5. Conclusions

In summary, we have proposed a general algorithm based on Lissajous figures and ellipse fitting to calculate the parameters of phaseshifted interferograms without using an iterative process. Using this method the initial phase, background intensity and modulation amplitude can be extracted. Numerical simulations and optical experiments have been performed, showing that the proposed method is noniterative, effective and accurate. In addition, the proposed algorithm gives a fast method to estimate the phase shifts. 150

Z. Zhai et al.

Optics Communications 429 (2018) 144–151

Acknowledgment

[9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

This work was supported by the National Natural Science Foundation of China (Nos. 51575164, 51405143), the Science and Technology Research Project of Department of Education of Hubei Province (No. D20161406), the China Postdoctoral Science Foundation (No. 2016M602269), and the Open Foundation of Hubei Collaborative Innovation Center for High-efficient Utilization of Solar Energy (No. HBSKFZD2014007). References [1] [2] [3] [4] [5] [6] [7] [8]

X. Jiang, D. Whitehouse, Phil. Trans. 370 (2012) 4154. W. Yang, X. Liu, W. Lu, X. Guo, Precis. Eng. 51 (2018) 348. L. Zhou, F. Zhou, D. Qu, X. Liu, W. Lu, Opt. Commun. 306 (2013) 174–178. W. Yang, X. Liu, W. Lu, C. Hu, Precis. Eng. 51 (2018) 40–47. W. Lu, X. Liu, L. Zhou, H. Heiyang, Measurement 45 (2012) 1633. B. Kemper, G. von Bally, Appl. Opt. 47 (2008) A52. X. Zhang, X. Zhang, M. Xu, Opt. Commun. 407 (2018) 367. Z. Lei, X. Liu, L. Zhao, L. Chen, Q. Li, T. Yuan, W. Lu, Opt. Commun. 359 (2016) 435.

151

L. Liu, M. Shan, Z. Zhong, B. Liu, G. Luan, M. Diao, Opt. Lett. 42 (2017) 4331. Z. Lei, X. Liu, L. Chen, W. Lu, S. Chang, Measurement 80 (2016) 1. R. Smythe, R. Moore, Opt. Eng. 23 (1984) 361. B. Hao, M. Shan, Z. Zhong, M. Diao, Y. Zhang, Opt. Lasers Eng. 51 (2013) 1278. Q. Kemao, H. Wang, W. Gao, Appl. Opt. 47 (2008) 5408. Z. Zhong, H. Bai, M. Shan, Y. Zhang, L. Guo, Opt. Lasers Eng. 97 (2017) 9. Z. Wang, B. Han, Opt. Lett. 29 (2004) 1671. W. Lu, M. Chang, P. Chen, W. Luo, J. Modern Opt. 61 (2014) S2. B. Kimbrough, Int. J. Optomechatronics 8 (2014) 304. G. Vladimir, I. Sergey, K. Roman, H. Dmitry, Opt. Eng. 52 (2013) 030501. C.J. Morgan, Opt. Lett. 7 (1982) 368. K. Okada, A. Sato, J. Tsujiuchi, Opt. Commun. 84 (1991) 118. Z. Wang, B. Han, Opt. Lasers Eng. 45 (2007) 274. Q. Kemao, H. Wang, W. Gao, L. Feng, S.H. Soon, Opt. Lasers Eng. 48 (2010) 684. R. Zhu, B. Li, R. Zhu, Y. He, J. Li, Opt. Express 24 (2016) 006814. C.T. Farrell, M.A. Player, Meas. Sci. Technol. 3 (1992) 953. C.T. Farrell, M.A. Player, Meas. Sci. Technol. 5 (1994) 648. F. Liu, Y. Wu, Opt. Express 23 (2015) 10794. F. Bookstein, Comput. Graphics Image Proc. 9 (1979) 56. K. Yatabe, K. Ishikawa, Y. Oikawa, Opt. Express 25 (2017) 29401. IEEE Std 1057-2007, IEEE Standard for Digitizing Waveform Recorders, http://doi. org/10.1109/IEEESTD.2008.8291822, 2008.