Accurate blind extraction of arbitrary unknown phase shifts by an improved quantum-behaved particle swarm optimization in generalized phase-shifting interferometry

Xuelian Yu; Renkang Song; Renkang Song; Xiufang Li; Pan Jiang; Kangwei Wang; Zhichang Zhang; Tao Shen

doi:10.1364/OSAC.2.003404

1. Introduction

Recently, digital phase-shifting interferometry (PSI) has been widely used as a powerful tool for the measurement of phase objects in studies of surface profiles of live cellular adhesion [1], surface measurement [2] and micro-structures characterization [3]. Standard PSI requires specific and constant phase shift ${{2\pi } \mathord{\left/ {\vphantom {{2\pi } N}} \right.} N}$, where $N \ge 3$ is an integer. This requirement, however, is often difficult to meet in realistic implementation. Therefore, additional efforts must be made to calibrate the phase shifter [4,5]. To circumvent this inconvenience, some phase-retrieval algorithms with unknown phase shifts have been developed [6,7]. Usually these methods have substantial computation loads, and some of them need dozens of (as many as 15) interferograms [6], or additional optical devices [7]. Generally, the phase-shifting algorithm (PSA) with more interferograms has high accuracy and low speed, and the PSA with less interferograms has high speed and low accuracy, it is difficult to obtain the high accuracy and speed simultaneously. To balance the computational time and accuracy, the research of non-iterative random PSA with less interferograms is essential. A series of wave-front reconstruction methods for two-frame phase-shifting interferometry are proposed [8–10], their method is fast and accurate, but the calculation method is complicated in the solution process. Cai et al. proposed the concept of generalized phase-shifting interference (GPSI) that can apply to arbitrary unknown phase shifts [11–12]. Several two-step GPSI methods [13–18] were introduced to achieve a higher precision of the extracted phase shift and a higher quality of the reconstructed image. Whether the method is based on the statistic property of the object phase distribution in the recording plane [11] or on the statistical nature of the diffraction field [12,13], it is necessary to introduce an error function. The unknown phase shifts are found by direct searching [12]. This iterative method heavily relies on the matching of the interferogram, and the exact phase shift is obtained after a dozen of iterations or even more. However, the step size choice of the phase shift iteration has a great influence on the calculation speed. If the step size is too small, more iterations may be required when increasing or decreasing the phase shifts. Whereas in the large step size case, it is easy to miss the minimum, so that one needs to reduce the step size and search in the opposite direction. Due to the problem of selecting the step size, the solution hovers around the minimum value, which is called the “sawtooth-solution phenomenon” in this paper. We manage to avoid this phenomenon by introducing the mutated adaptive quantum-behaved particle swarm optimization (MAQPSO) algorithm into the GPSI method.

In 2004, Sun et al. proposed a quantum-behaved particle swarm optimization (QPSO) algorithm that uses particle behavior in quantum space to simulate human behavior [19,20]. It has high randomness and can perform global search. Over the past decade, QPSO algorithm and its improved versions have been widely used to solve various complex problems in a variety of fields such as financial distribution [21], biomedicine [22], and fiber optics [23]. The AQPSO algorithm is developed from quantum-behaved particle swarm optimization (QPSO) which has been successfully used to numerical optimization [24], design FG plates [25], and engineering inverse problem [26]. The AQPSO algorithm has no velocity vector for the particles and the whole feasible solution space can be searched for the index modulation; only a parameter associated with the process of the iteration need to be adjusted. Moreover, the iterative equation of this algorithm is very different from that of PSO. The MAQPSO algorithm is an adaptive quantum-behaved particle swarm optimization (AQPSO) with mutation operation to improve the search performance of the AQPSO algorithm [27]. Thus, the MAQPSO algorithm is much simple and easy to implement. Here we introduce a scheme to solve the “sawtooth-solution phenomenon” caused by the direct search using MAQPSO algorithm [27]. It turns out that we can extract any unknown and unequal phase steps faster and accurately.

2. Principles

2.1. Principle of generalized phase-shifting interferometry algorithm

In many applications of PSI, the test object is located at certain distance from the recording (CCD) plane, and the CCD actually records the diffraction field from the object surface. In this case the statistical properties of the diffraction field may provide helpful hints for data processing. Using general method of PSI with arbitrary unknown phase shifts can enable extracting the phase steps and reconstructing the object field using only the measured constant reference wave intensity. This method can be valid for frame number $N \ge 3$.

Let us firstly describe the principle of the proposed method, by taking the four-step phase-shifting as an example. Assume that the diffraction field of an object in the recording plane ${\textrm{P}_{\textrm{H}}}$ is $O({x,y} )= {A_o}({x,y} )\exp [{i\varphi ({x,y} )} ]$, Where ${A_o}({x,y} )$ and $\varphi ({x,y} )$ are the real amplitude and phase distributions of the object wave on the recording plane ${\textrm{P}_{\textrm{H}}}$, respectively, and the reference wave is a plane wave with constant amplitude ${A_r}$; the intensity distributions of four interferograms can be written as

(1)$$\begin{array}{l} {I_0} = {A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos [\varphi ]\\ {I_1} = {A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos [{\varphi - {\alpha_1}} ]\\ {I_2} = {A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos [{\varphi - {\alpha_1} - {\alpha_2}} ]\\ {I_3} = {A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos [{\varphi - {\alpha_1} - {\alpha_2} - {\alpha_3}} ]\end{array}$$

where ${\alpha _1},{\alpha _2}$ and ${\alpha _3}$ ($0 \,< \,{a_j} \,< \,\pi ,j = 1,2,3$) are three constant phase steps. Note that the coordinates $({x,y} )$ are omitted for conciseness. It is possible to introduce three more quantities by the four measured intensities in Eq. (1)

(2)$$\begin{array}{l} p = \left\langle {|{{I_1} - {I_0}} |} \right\rangle = \left\langle {|{4{A_r}{A_o}\sin ({\varphi - {{{\alpha_1}} \mathord{\left/ {\vphantom {{{\alpha_1}} 2}} \right.} 2}} )} |} \right\rangle \sin ({{{{\alpha_1}} \mathord{\left/ {\vphantom {{{\alpha_1}} 2}} \right.} 2}} )\\ q = \left\langle {|{{I_2} - {I_1}} |} \right\rangle = \left\langle {|{4{A_r}{A_o}\sin ({\varphi - {\alpha_1} - {{{\alpha_2}} \mathord{\left/ {\vphantom {{{\alpha_2}} 2}} \right.} 2}} )} |} \right\rangle \sin ({{{{\alpha_2}} \mathord{\left/ {\vphantom {{{\alpha_2}} 2}} \right.} 2}} )\\ r = \left\langle {|{{I_3} - {I_2}} |} \right\rangle = \left\langle {|{4{A_r}{A_o}\sin ({\varphi - {\alpha_1} - {\alpha_2} - {{{\alpha_3}} \mathord{\left/ {\vphantom {{{\alpha_3}} 2}} \right.} 2}} )} |} \right\rangle \sin ({{{{\alpha_3}} \mathord{\left/ {\vphantom {{{\alpha_3}} 2}} \right.} 2}} )\end{array}$$

where $\left\langle {} \right\rangle$ represents the average over the whole frame. Suppose that ${A_o}$ and $\varphi$ are mutually independent and $\varphi$ is a spatially random distribution due to the Fresnel diffraction, one can reasonably expect all the terms $\left\langle {|{\ldots } |} \right\rangle$ on the right-hand side of Eq. (2) have approximately the same value $c$. Then we have

(3)$$p = c\sin ({{{{\alpha_1}} \mathord{\left/ {\vphantom {{{\alpha_1}} 2}} \right.} 2}} ),q = c\sin ({{{{\alpha_2}} \mathord{\left/ {\vphantom {{{\alpha_2}} 2}} \right.} 2}} ),r = c\sin ({{{{\alpha_3}} \mathord{\left/ {\vphantom {{{\alpha_3}} 2}} \right.} 2}} )$$

When the intensity of the object light and the reference light are recorded, the averaged contribution $c$ can be directly calculated to obtain an approximation, and a simple transformation can be established

(4)$${\alpha _1} = 2arc\sin ({{p \mathord{\left/ {\vphantom {p c}} \right.} c}} ),{\alpha _2} = 2arc\sin ({{q \mathord{\left/ {\vphantom {q c}} \right.} c}} ),{\alpha _3} = 2arc\sin ({{r \mathord{\left/ {\vphantom {r c}} \right.} c}} )$$

From Eq. (4), we immediately get an approximation of the phase shifts without adding more error functions. This result can be further improved by the improved mutated adaptive quantum-behaved particle swarm optimization. Firstly, we substitute ${\alpha _1},{\alpha _2}$ and ${\alpha _3}$ in Eq. (4) into Eq. (5) which was directly derived from the expressions of ${I_j}({j = 1,2,\ldots ,N} )$ [12],

(5)$$\begin{aligned}{A_o}\exp (i\varphi ) &= \frac{1}{{4{A_r}\sin [{{{({{\alpha_1} + {\alpha_3}} )} \mathord{\left/ {\vphantom {{({{\alpha_1} + {\alpha_3}} )} 2}} \right.} 2}} ]}}\\ &\times \left\{ {\frac{{\exp [{i{{({{\alpha_1} + {\alpha_2}} )} \mathord{\left/ {\vphantom {{({{\alpha_1} + {\alpha_2}} )} 2}} \right.} 2}} ]}}{{\sin [{{{({{\alpha_2} + {\alpha_3}} )} \mathord{\left/ {\vphantom {{({{\alpha_2} + {\alpha_3}} )} 2}} \right.} 2}} ]}}({{I_1} - {I_3}} )} \right.\left. { - \frac{{\exp [{i({{\alpha_1} + {{{\alpha_2}} \mathord{\left/ {\vphantom {{{\alpha_2}} 2}} \right.} 2} + {{{\alpha_3}} \mathord{\left/ {\vphantom {{{\alpha_3}} 2}} \right.} 2}} )} ]}}{{\sin [{{{({{\alpha_1} + {\alpha_2}} )} \mathord{\left/ {\vphantom {{({{\alpha_1} + {\alpha_2}} )} 2}} \right.} 2}} ]}}({{I_0} - {I_2}} )} \right\} \end{aligned}$$

Both the amplitude ${A_o}({x,y} )$ and the phase $\varphi ({x,y} )$ can be obtained by the measured constant intensity ${I_r} = {A_r}^2$, and we construct an error function

(6)$$\varepsilon = \sum {\left\{ \begin{array}{l} \textrm{|}{I_0} - ({A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos \varphi )\textrm{|}\\ + \textrm{|}{I_1} - [{A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos (\varphi - {\alpha_1})]\textrm{|}\\ + \textrm{|}{I_2} - [{A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos (\varphi - {\alpha_1} - {\alpha_2})]\textrm{|}\\ + \textrm{|}{I_3} - [{A_o}^2 + {A_r}^2 + 2{A_r}{A_o}\cos (\varphi - {\alpha_1} - {\alpha_2} - {\alpha_3})]\textrm{|} \end{array} \right\}}$$

where $\sum$ denotes summation over four frames. Theoretically $\varepsilon$ is a smooth function of $({{\alpha_1},{\alpha_2},{\alpha_3}} )$ and would approach to its minimum (zero) when ${\alpha _1},{\alpha _2}$ and ${\alpha _3}$ takes their exact values used to record the interferograms. Therefore, we can gradually change ${\alpha _1},{\alpha _2}$ and ${\alpha _3}$ by the MAQPSO algorithm to reduce $\varepsilon$ till its minimum.

2.2. Principle of the mutated adaptive quantum-behaved particle swarm optimization

The adaptive quantum-behaved particle swarm optimization (AQPSO) algorithm was proposed by Zhan et al. in 2009 [27]. Adaptive refers to dynamically updating the key parameters used in the algorithm and correctly performing policy calls, transformations, and settings based on the evolution state of the algorithm. The purpose of adaptation is to make the algorithm efficiently find the optimal solution or a satisfactory suboptimal solution. The AQPSO algorithm identifies the evolutionary state of a population in real time based on the distance distribution between each particle and all other particles. Then, the inertia weight and acceleration coefficient are adaptively updated according to the evolution state of the population to accelerate convergence. The MAQPSO algorithm uses Gaussian perturbation to apply appropriate momentum to the global optimal position to help the algorithm escape the stagnation/local extremum.

In order to increase the population diversity of the initial value of phase shifts and improve the convergence speed of the algorithm, we added the first mutation operator with Gaussian probability distribution before running the MAQPSO algorithm. The mutation operator reads

(7)$${p_{ed}}(m )= {p_{ed}}(m )\pm N({0,{\xi_1}\cdot \Delta {n_{range}}} )$$

where $N({0,{\xi_i}\cdot \Delta {n_{range}}} )$ is a normal distribution random number with a mean of 0 and a variance of ${\xi _i}\cdot \Delta {n_{range}}$. Here $\Delta {n_{range}}$ represents the modulation range of the initial ${\alpha _1},{\alpha _2},{\alpha _3}$, typically as $[{ - 0.005,0.005} ]$ (In this paper we set the iteration step precision to 0.0001, so the initial population setting is a higher precision range) and ${\xi _i}$ is the variation constant, taking the random constant inside $({0,1} )$.

Before we proceed, let us first outline the MAQPSO method. In the AQPSO algorithm, the particles move according to the following iterative equations [16,17]:

(8)$${p_{ed}}(m )= {\varphi _d}(m ){p_{ed}}(m )+ [{1 - {\varphi_d}(m )} ]{p_{gd}}(m )$$

(9)$${m_{best}}(m )= \frac{1}{M}\sum\limits_{e = 1}^M {{p_{gd}}(m )}$$

(10)$${x_{ed}}({m + 1} )= {p_{ed}}(m )\pm \beta (m )|{{m_{best}}(m )- {x_{ed}}(m )} |\ln [{{1 \mathord{\left/ {\vphantom {1 {{u_{ed}}(m )}}} \right.} {{u_{ed}}(m )}}} ]$$

where ${p_{ed}}$ is the attraction point of the particle swarm, ${p_{ed}}$ and ${p_{gd}}$ are the best position of the $e\textrm{ - th}$ particle and the position of the best particle among all the particles in a D-dimensional space, respectively; ${m_{best}}$ is defined as the mean value of all particles' ${p_{ed}}$, and ${u_{ed}}$ and ${\varphi _d}$ are random number distributed uniformly on $[{0,1} ]$, respectively; M is the total particle number, and $m$ denotes the iterative generation. The contraction-expansion coefficient (CEC) $\beta$ is the only adjustable parameter used to control the convergence speed of the algorithm.

In order to improve the search performance of the AQPSO algorithm, we introduce the second mutation operator with Gaussian probability distribution the ${p_{ed}}$ as follows:

(11)$${p_{ed}}(m )= {p_{ed}}(m )\pm N({0,{\xi_i}\cdot \Delta {n_{range}}} )$$

where $N({0,{\xi_i}\cdot \Delta {n_{range}}} )$ is among a set of normal random numbers with mean 0 and standard deviation ${\xi _i}\cdot \Delta {n_{range}}$, here $\Delta {n_{range}}$ denotes the range of the initial index modulation, and ${\xi _i}$ is a mutated positive constant less than 1.

The ways to update ${p_{ed}}$ and ${p_{gd}}$ are identical to the corresponding ones in a generic PSO [28], namely

(12)$${p_{ed}}({m + 1} )= \left\{ \begin{array}{l} {x_{ed}}({m + 1} ),({f[{{x_{ed}}({m + 1} )} ]< f[{{p_{ed}}(m )} ]} )\\ {p_{ed}}(m ),({f[{{x_{ed}}({m + 1} )} ]\ge f[{{p_{ed}}(m )} ]} )\end{array} \right.$$

(13)$${p_{gd}}({m + 1} )= \arg \mathop {\min }\limits_{1 \le e \le M} \{{f[{{p_{ed}}(m )} ]} \}$$

where $f$ is the fitness value function used to evaluate each particle in the swarm, as shown in Eq. (6).

To accelerate the proposed method, the CEC parameter $\beta$ is adaptive in the process of each iteration. To do that, the following error function is used to identify how close the particle is to the global best position

(14)$$F(m )= {{[{{f_e}(m )- {f_{gbest}}(m )} ]} \mathord{\left/ {\vphantom {{[{{f_e}(m )- {f_{gbest}}(m )} ]} {\min \{{abs[{{f_e}(m )} ],abs[{{f_{gbest}}(m )} ]} \}}}} \right.} {\min \{{abs[{{f_e}(m )} ],abs[{{f_{gbest}}(m )} ]} \}}}$$

where ${f_e}$ is the fitness value of the $e\textrm{ - th}$ particle and ${f_{gbest}}$ is the fitness value of the best particle, $abs[{{f_e}} ]$ and $abs[{{f_{gbest}}} ]$ denoting the absolute value of ${f_e}$ and ${f_{gbest}}$, respectively. Furthermore, $\min \{{abs[{{f_e}} ],abs[{{f_{gbest}}} ]} \}$ refers to the smaller one of $abs[{{f_e}} ]$ and $abs[{{f_{gbest}}} ]$. In this regard, the smaller the value of the error function $F(m )$ for a certain particle, the closer it to the best particle. The particles that are far away from the best particle should be given a smaller $\beta$, whereas those close to the best particle should be given larger value of $\beta$. In other words, $\beta$ is designed to be a nonlinear function $\beta ({\log F,e,m} )$. In our application, when changing log F from small value to large value, $\beta ({\log F,e,m} )$ are set from 0.3 to 0.01 which are different from that in Ref. [17]. This is because that the bigger αwill decrease the converging velocity of the MAQPSO algorithm.

In order to increase the population diversity of the initial value of phase shifts and improve the convergence speed of the algorithm, we add a mutation operator with Gaussian probability distribution in the MAQPSO algorithm. The mutation operator reads

(15)$${p_{ed}}(m )= {p_{ed}}(m )\pm N({0,{\xi_1}\cdot \Delta {n_{range}}} )$$

where $N({0,{\xi_i}\cdot \Delta {n_{range}}} )$ is a normal distribution random number with a mean of 0 and a variance of ${\xi _i}\cdot \Delta {n_{range}}$. Here $\Delta {n_{range}}$ represents the modulation range of the initial ${\alpha _1},{\alpha _2},{\alpha _3}$, typically as $[{ - 0.0001,0.0001} ]$ (In this paper, we set the iteration step precision to 0.0001) and ${\xi _i}$ is the variation constant, taking the random constant inside $({0,1} )$.

The improved MAQPSO algorithm is briefly summarized as follows:

1) The initial phase shifts ${\alpha _1},{\alpha _2},{\alpha _3}$ is obtained according to the Eq. (4), and initialize a population of particles with random positions ${x_{ed}}(0 )$ inside the phase shifts modulation range $\Delta {n_{range}}$. The position ${x_{ed}}(0 )$ of each particle represents the value of a set of phase shifts;
2) Evaluate the fitness value of each particle by substituting its position into Eq. (6) and assign the particle's best position to ${p_{ed}}(0 )$. Identify the best among ${p_{ed}}(0 )$ as ${p_{gd}}(0 )$. According to Eq. (8), ${p_{ed}}(0 )$ is calculated, and ${m_{best}}(0 )$ is also obtained by Eq. (9);
3) The CEC parameterα is chosen according to Eq. (14);
4) Updating the position of the particle population ${x_{ed}}(m )$ by Eq. (10);
5) Evaluate the fitness value of each particle by substituting its position into Eq. (6) and assign the particle's best position to ${p_{ed}}(m )$. Identify the best among ${p_{ed}}(m )$ as ${p_{gd}}(m )$;
6) The particle's best position ${p_{ed}}(m )$ is operated with the mutation probability and its fitness value is compared with the one before the mutation operation. If its fitness value is better, then update it and the fitness value;
7) Calculate ${p_{ed}}(m )$ and ${m_{best}}(m )$ by Eqs. (8) and (9), respectively;
8) For each particle, the fitness value of the particle's current position is calculated. If it is better than the previous position, then update it and the fitness value. Namely, if $f[{{x_{ed}}({m + 1} )} ]< f[{{p_{ed}}(m )} ]$, then ${p_{ed}}({m + 1} )= {x_{ed}}({m + 1} )$ and $f[{{p_{ed}}({m + 1} )} ]\textrm{ = }f[{{x_{ed}}({m + 1} )} ]$;
9) Identify the best among ${p_{ed}}(m )$ as the global best position of the current particle population ${p_{gd}}({m + 1} )$ and update the fitness value, namely ${p_{gd}}({m + 1} )= \arg {\min _{1 \le e \le M}}\{{f[{{p_{ed}}(m )} ]} \}$ and $f[{{p_{gd}}({m + 1} )} ]= f[{{p_{gd}}(m )} ]$;
10) Compare the global best position of the current particle population with the previous global best position, if it is better than $f[{{p_{gd}}({m - 1} )} ]$, update ${p_{gd}}(m )$ and its fitness value;
11) Repeat steps 3-10 until the assigned fitness value function is achieved.

Combining traditional method based on the statistical nature of the diffraction field and the MAQPSO algorithm, a new method is developed as shown in Fig. 1.

Fig. 1. Flow chart of the proposed method.

Download Full Size | PDF

3. Simulation

We carried out a series of simulations to verify the effectiveness of the proposed method. The schematic diagram of the device in the simulation is shown in Fig. 2. In these simulations, an artificially designed amplitude type object with a Chinese word was used as the test object. An axial plane wave of wavelength $\lambda = 632.8\textrm{nm}$ illuminates the surface and transmitted into a CMOS plane of $512 \times 512$ pixels with each pixel being 5.2µm × 5.2µm, the distance between the object plane and recording plane is $Z = 200\textrm{mm}$, and the diffractive surface (object) has a size of 5 mm (all the parameters were chosen according to sampling theory). Two-dimensional Fresnel diffraction of the complex amplitude of original object plane ${\textrm{P}_{\textrm{O}}}$ onto plane ${\textrm{P}_{\textrm{H}}}$ serves as object wave $O({x,y} )$. Different interferograms are computer generated with different nonstandard phase shifts.

Fig. 2. The device schematic diagram in the simulation experiment.

Download Full Size | PDF

For demonstration, here we present calculated results of the four-frame case in Fig. 3. In these simulations we set $\; {\alpha _1} = 1.39626({80^\circ } )$, ${\alpha _2} = 1.83260({105^\circ } )$ and ${\alpha _3} = 1.65806({95^\circ } )$. The initial values of ${\alpha _1},{\alpha _2}$ and ${\alpha _3}$ calculated directly with Eq. (4), without added the error function, which greatly reduces the computational complexity. It can be seen from Figs. 3(a)–(c) that the initial value calculated by our method have already reached quite high accuracy (the relative error is less than 0.7% in the worst case). Gaussian mutation operation is performed according to Eq. (7). Once the initial particle swarm is obtained, the accurate value of the phase shifts is calculated using the MAQPSO.

Fig. 3. Comparison of phase shifts extraction results of our method and the direct search for a four-frame PSI. (a)–(c) The iterative process of ${\alpha _1},{\alpha _2},{\alpha _3}$, (d) the number of repetitive experiments of our method and the direct search [12].

Download Full Size | PDF

Figures 3(a)–(c) clearly show that ${\alpha _1},{\alpha _2}$ and ${\alpha _3}$ can be retrieved almost exactly after 9 iterations. Note that we set the significant figure to 6 digits in the algorithm proposed in this paper, and the traditional direct search to 5 digits. Since the traditional direct search can only be iterated in steps of 0.1, 0.01, 0.001, and 0.0001 [12], the computation time of each additional significant figure is tremendously huge. On the contrary, the algorithm we proposed can quickly converge near the exact value, even with a large step size in the initial stage of the iteration. The step size is reduced in the iteration process to search for the exact value, avoiding the “sawtooth-solution phenomenon”. The calculation efficiency is not obviously affected because of the significant figure of the accurate value. In summary, compared with the traditional direct search method, our proposed method can improve the calculation accuracy without significantly reducing the computational efficiency, and greatly improve the practical application value of the algorithm. In order to test the stability of our method, we performed 20 repetitive simulation experiments. The results of the repetitive experiments are shown in Fig. 3(d). The processing time necessary using the proposed method is 1.4 s, while the direct search used in [12] lasts for more than 25 s (all results are obtained with a 3.2 GHz laptop and MATLAB 2018a).

Finally, we show the ability of this algorithm to retrieve the object phase distribution. The results are shown in Fig. 4, where the up and down rows correspond to two cases (the object phase and the difference). Figures 4(f) and 4(g) show the retrieved object phase by our method and the direct search approach [12]. We calculate the cross-correlation function between the object phases which are 0.9146, 0.8759 and 0.8615, respectively, as shown in Figs. 4(h)–(j). Figures 4(i) and 4(j) show that the retrieved object phase by our method has a higher correlation with the original object phase than the direct search.

Fig. 4. Simulation results of object phase reconstruction with the four-frame method with retrieved nonstandard phase shifts. (a) Original object phase, (b)–(e) four holograms with different phase shifts obtained on the image plane of the camera, namely ${I_1}$, ${I_2}$, ${I_3}$ and ${I_4}$, (f) retrieved object phase by our method, (g) retrieved object phase by direct search, (h) the difference between (f) and (g), (i) the difference between (a) and (f), (j) the difference between (a) and (g).

Download Full Size | PDF

4. Demonstration with experimental data

Some optical experiments of GPSI have also been carried out to investigate the performance of our method. The experimental setup is similar to that shown in Fig. 4 of Ref. [17] but use a solid-state laser. A solid-state laser of 532 nm is used as the light source. A transparency of United States Air Force (USAF) resolution target is used as the amplitude object and the reference beam is a normal plane wave. The central part of $1024 \times 1280$ pixels of a CMOS with pixel size of 5.2µm × 5.2µm is chosen as the frame size, and only its central part of $940 \times 940$ pixels is chosen as the frame size. The distance between the recording plane and object plane is $Z = 216\textrm{mm}$. With the preset phase shifts of $\; {\alpha _1} = 1.39626({80^\circ } )$, ${\alpha _2} = 1.83260({105^\circ } )$ and ${\alpha _3} = 1.65806({95^\circ } )$ we obtained four interferograms as shown in Figs. 5(a)–(d) (${I_1}$, ${I_2}$, ${I_3}$ and ${I_4}$, respectively), and Figs. 5(e) and 5(f) are the recorded intensity maps of ${I_o}$ and ${I_r}$. The extracted result of the phase shifts by our method are $\; {\alpha _1} = 1.39681$, ${\alpha _2} = 1.83325$ and ${\alpha _3} = 1.65843$, with a very small difference from its nominal value. The objective wave is retrieved with these value of $\; {\alpha _1} = 1.39681$, ${\alpha _2} = 1.83325$ and ${\alpha _3} = 1.65843$ by Eq. (5) in recording plane ${\textrm{P}_{\textrm{H}}}$, and then reconstructed in object plane ${\textrm{P}_{\textrm{O}}}$ by an inverse Fresnel transform. The direct search method uses the same process for calculation. The target image reconstructed with the direct search and the proposed method are shown in Figs. 5(g) and 5(h), respectively. Although there are some stains in the intensity maps of ${I_1}$, ${I_2}$, ${I_3}$ and ${I_4}$, their effects are considerably offset in this method and we can see that the reconstructed image has high quality with little noise.

In order to compare the spatial resolution of the images reconstructed by the direct search and the proposed method, we magnify the central regions of Figs. 5(g) and 6(h), respectively, as shown in Figs. 6(b) and 6(e). Then, the intensity distribution map corresponding to the red line in Figs. 6(b) and 6(e) is obtained, as shown in Figs. 6(c) and 6(f). It can be seen from Figs. 6(c) and 6(f) that the three peaks and valleys (the position of the ellipse) of the fifth pair of the fourth group can be clearly distinguished, so the spatial resolution of the digital holographic imaging system in this paper is 19.7 µm. Comparing the PV values (the difference between the highest point and the lowest point of a pair of peaks and valleys in the curve) of the three peaks and valleys of the fifth group of the fourth group in Figs. 6(c) and 6(f), the PV value of Fig. 6(f) is larger than the PV value in Fig. 6(c). Figure 6 demonstrates that phase extraction using the proposed method results in a more accurate phase shifts, which is advantageous for obtaining a higher quality reconstructed image.

Fig. 5. Optical experiment results. (a)–(d) the four interferograms ${I_1}$, ${I_2}$, ${I_3}$ and ${I_4}$ in plane ${\textrm{P}_{\textrm{H}}}$, (e) and (f) intensity ${I_o}$ of object wave and ${I_r}$ of reference wave in plane ${\textrm{P}_{\textrm{H}}}$, (g) and (h) the reconstructed images in plane ${\textrm{P}_{\textrm{O}}}$ by the direct search and the proposed method, respectively.

Download Full Size | PDF

Fig. 6. Comparison of reconstructed images by the direct search and the proposed method. (a) and (d) the reconstructed images in plane ${\textrm{P}_{\textrm{O}}}$ by the direct search and the proposed method, respectively, (b) and (e) enlarged views of the central part of (a) and (d), respectively, (c) and (f) the intensity distribution maps corresponding to the red lines at (b) and (e), respectively.

Download Full Size | PDF

5. Conclusion

In this paper, we have proposed a simplified but generalized algorithm for accurate and efficient retrieve of actual phase steps with arbitrary and unequal phase shifts in the PSI. The effectiveness of this algorithm is verified by a series of computer simulations and the optical experiments. Compared with existing methods, this algorithm can quickly converge to the exact solution, and can be used to avoid the “sawtooth-solution phenomenon” encountered in the traditional direct search approach. More importantly, this algorithm is fast (usually within 10 iterations in 1.408s, literally 18 times faster than the direct search), being able to yield results with a high degree of accuracy. Optical experiments have further yielded very satisfactory wave-front retrieval results. We believe that this proposed method can improve the practical applications of PSI for wave-front measurement and reconstruction.

Funding

National Natural Science Foundation of China (51677044); Natural Science Foundation of Heilongjiang Province (E2018047).

Acknowledgement

This work was supported by the Heilongjiang Provincial Key Laboratory of Quantum Manipulation & Control.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. W. M. Ash, L. Krzewina, and M. K. Kim, “Quantitative imaging of cellular adhesion by total internal reflection holographic microscopy,” Appl. Opt. 48(34), H144–152 (2009). [CrossRef]

2. V. Srivastava, T. Anna, and D. S. Mehta, “Full-field Hilbert phase microscopy using nearly common-path low coherence off-axis interferometry for quantitative imaging of biological cells,” J. Opt. 14(12), 125707 (2012). [CrossRef]

3. R. J. Edwards, S. W. Zhou, S. J. Hwang, K. Y. Mckeown, B. Wang, R. Bhaduri, P. J. Ganti, A. G. Yunker, J. A. Yodh, L. L. Rogers, G. Goddard, and Popescu, “Diffraction phase microscopy: monitoring nanoscale dynamics in materials science [invited],” Appl. Opt. 53(27), G33–43 (2014). [CrossRef]

4. N. A. Ochoa and J. M. Huntley, “Convenient method for calibrating nonlinear phase modulators for use in phase-shifting interferometry,” Opt. Eng. 37(9), 2501–2505 (1998). [CrossRef]

5. Y. J. Shen and J. M. Huntley, “Simple method to calibrate phase modulators for use in dynamic phase-shifting interferometry,” Opt. Eng. 43(12), 2998–3002 (2004). [CrossRef]

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

7. H. Kadono, Y. Bitoh, and S. Toyoka, “Statistical interferometry based on a fully developed speckle field: an experimental demonstration with noise analysis,” J. Opt. Soc. Am. A 18(6), 1267–1274 (2001). [CrossRef]

8. Z. T. Cheng and D. Liu, “Fast and accurate wavefront reconstruction in two-frame phase-shifting interferometry with unknown phase step,” Opt. Lett. 43(13), 3033–3036 (2018). [CrossRef]

9. Y. Zhang, X. B. Tian, and R. G. Liang, “Accurate and fast two-step phase shifting algorithm based on principle component analysis and Lissajous ellipse fitting with random phase shift and no pre-filtering,” Opt. Express 27(14), 20047–20063 (2019). [CrossRef]

10. Y. Zhang, X. B. Tian, and R. G. Liang, “Random two-step phase shifting interferometry based on Lissajous ellipse fitting and least squares technologies,” Opt. Express 26(12), 15059–15072 (2018). [CrossRef]

11. L. Z. Cai, Q. Liu, X. L. Yang, and Y. R. Wang, “Phase-shift extraction and wave-front reconstruction in phase-shifting interferometry with arbitrary phase steps,” Opt. Lett. 28(19), 1808–1810 (2003). [CrossRef]

12. L. Z. Cai, Q. Liu, and X. L. Yang, “Generalized phase-shifting interferometry with arbitrary unknown phase steps for diffraction objects,” Opt. Lett. 29(2), 183–185 (2004). [CrossRef]

13. X. F. Xu, L. Z. Cai, X. F. Meng, G. Y. Dong, and X. X. Shen, “Fast blind extraction of arbitrary unknown phase shifts by an iterative tangent approach in generalized phase-shifting interferometry,” Opt. Lett. 31(13), 1966–1968 (2006). [CrossRef]

14. X. F. Xu, L. Z. Cai, Y. R. Wang, X. L. Yang, X. F. Meng, G. Y. Dong, X. X. Shen, and H. Zhang, “Generalized phase-shifting interferometry with arbitrary unknown phase shifts: Direct wave-front reconstruction by blind phase shift extraction and its experimental verification,” Appl. Phys. Lett. 90(12), 121124 (2007). [CrossRef]

15. X. F. Xu, L. Z. Cai, Y. R. Wang, X. F. Meng, W. J. Sun, H. Zhang, X. C. Cheng, G. Y. Dong, and X. X. Shen, “Simple direct extraction of unknown phase shift and wavefront reconstruction in generalized phase-shifting interferometry: algorithm and experiments,” Opt. Lett. 33(8), 776–778 (2008). [CrossRef]

16. X. F. Xu, L. Z. Cai, Y. R. Wang, and D. L. Li, “Accurate phase shift extraction for generalized phase-shifting interferometry,” Chin. Phys. Lett. 27(2), 024215 (2010). [CrossRef]

17. X. F. Xu, L. Z. Cai, Y. R. Wang, and R. S. Yan, “Direct phase shift extraction and wavefront reconstruction in two-step generalized phase-shifting interferometry,” J. Opt. 12(1), 015301 (2010). [CrossRef]

18. X. F. Xu, L. Z. Cai, F. Gao, Y. L. Jia, and H. Zhang, “Detection and correction of wavefront errors caused by slight reference tilt in two-step phase-shifting digital holography,” Appl. Opt. 54(32), 9591 (2015). [CrossRef]

19. J. Sun, B. Feng, and W. B. Xu, “Particle swarm optimization with particles having quantum behavior,” Cong. Evolution. Comp.1, 325–331 (2004).

20. J. Sun, W. B. Xu, and B. Feng, “Adaptive parameter control for quantum-behaved particle swarm optimization on individual level,” IEEE Int. Conf. Systems, Man and Cybernetics4, 3049–3054 (2005).

21. Q. Niu, Z. Zhou, H. Y. Zhang, and J. Deng, “An Improved Quantum-Behaved Particle Swarm Optimization Method for Economic Dispatch Problems with Multiple Fuel Options and Valve-Points Effects,” Energies 5(9), 3655–3673 (2012). [CrossRef]

22. H. Mirvaziri and Z. S. Mobarakeh, “Improvement of EEG-based motor imagery classification using ring topology-based particle swarm optimization,” Biomed. Signal Proces. 32, 69–75 (2017). [CrossRef]

23. X. L. Yu, Y. Yao, J. J. Xiao, and J. J. Tian, “Optimal design of short fiber Bragg gratings with triangular spectrum,” Opt. Commun. 285(5), 631–637 (2012). [CrossRef]

24. J. He and H. Guo, “A Modified Particle Swarm Optimization Algorithm,” Int. J. Comput. Int. Sys. 11(10), 151–155 (2013).

25. C. Wang, T. T. Yu, J. L. Curiel-Sosa, N. G. Xie, and T. Q. Bui, “Adaptive chaotic particle swarm algorithm for isogeometric multi-objective size optimization of FG plates,” Struct Multidisc Optim 60(2), 757–778 (2019). [CrossRef]

26. S. Khan, S. Yang, and O. U. Rehman, “A dynamic particle swarm optimization method applied to global optimizations of engineering inverse problem,” COMPEL 37(1), 98–117 (2018). [CrossRef]

27. Z. H. Zhan, J. Zhang, and Y. Li, “Adaptive particle swarm optimization,” IEEE Trans. Syst., Man, Cybern. B 39(6), 1362–1381 (2009). [CrossRef]

28. J. Kennedy and R. C. Eberhart, “Particle swarm optimization,” IEEE Int. Conf. on Neural Networks, 1942–1948 (1995).

Accurate blind extraction of arbitrary unknown phase shifts by an improved quantum-behaved particle swarm optimization in generalized phase-shifting interferometry

Abstract

1. Introduction

2. Principles

2.1. Principle of generalized phase-shifting interferometry algorithm

2.2. Principle of the mutated adaptive quantum-behaved particle swarm optimization

3. Simulation

4. Demonstration with experimental data

5. Conclusion

Funding

Acknowledgement

Disclosures

References

Cited By

Figures (6)

Equations (15)

OSA Continuum