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

Phase extraction from interferograms with unknown tilt phase shifts based on a regularized optical flow method

Open Access Open Access

Abstract

A novel method is presented to extract phase distribution from phase-shifted interferograms with unknown tilt phase shifts. The proposed method can estimate the tilt phase shift between two temporal phase-shifted interferograms with high accuracy, by extending the regularized optical flow method with the spatial image processing and frequency estimation technology. With all the estimated tilt phase shifts, the phase component encoded in the interferograms can be extracted by the least-squares method. Both simulation and experimental results have fully proved the feasibility of the proposed method. Particularly, a flat-based diffractive optical element with quasi-continuous surface is tested by the proposed method with introduction of considerably large tilt phase shift amounts (i.e., the highest estimated tilt phase shift amount between two consecutive frame reaches 6.18λ). The phase extraction result is in good agreement with that of Zygo’s MetroPro software under steady-state testing conditions, and the residual difference between them is discussed. In comparison with the previous methods, the proposed method not only has relatively little restrictions on the amounts or orientations of the tilt phase shifts, but also works well with interferograms including open and closed fringes in any combination.

©2013 Optical Society of America

1. Introduction

The temporal phase-shifting interferometry (PSI) has been generally accepted as the most accurate technique for wave-front reconstruction based on automated interferogram analysis [1]. PSI electronically records a series of interferograms while the reference phase of the interferometer is changed, so as to extract the phase information encoded in the variations in the intensity pattern of the recorded interferograms [2]. However, the accuracy of PSI is limited by the phase-shifting uncertainties resulting from the imperfect conditions in the real testing environments, such as the miscalibration, non-linear responses, aging of the phase shifter [3], and the mechanical vibration [1].

As a result, much effort has been made to improve the PSI methods, to alleviate the influence of phase-shifting uncertainties on phase extraction accuracy [1, 315]. The papers [310] assume the phase steps are constant for all pixels in one interferogram but typically different between interferograms. However, sometimes such an assumption may be invalid, for instance, due to the unbalanced piezoelectric effect in the phase shifter or instability of the optical platform, the phase shifter may probably introduces non-negligible orientation errors during the shift [1, 1115]; or as the extreme situations discussed in the paper [1], the strong environmental vibration will induce quite large tilt phase shift due to the relative motions between the test surface and the interferometer. In all those cases, the phase steps related to a specific interferogram are not longer constant for all pixels but vary with a tilt function across the field.

Several methods have been proposed to compensate the tilt-shift errors [1115]. The method given in the paper [11] is based on a first-order Taylor series expansion of the phase shift errors, which iteratively update the parameters related to the tilt phase shifts until numerical convergence and the required accuracy are achieved. However, as the ratio of the tilt-shift errors to the unknown phase steps increases the first-order approximation will gradually lose the accuracy, then the compensation performance of the tilt-shift errors will get worse. In the papers [12, 14] the interferograms are divided into a set of small blocks, and the phase steps within an individual block are viewed as an unknown constant. In the paper [12] these constants are solved by using calculated contrast maps, while in the paper [14] they are optimized by the least-squares method [4]. The size of the block in the papers [12, 14] should be reasonably defined beforehand. The size of the block would be very small if the tilt-shift amount is large. In that case, the estimation of tilt-shift parameters will be much coarse and even untruthful. The paper [13] proposes a method to extract the tilt-shift parameters from the set of zero-crossing points of the difference between the fringe intensities. However, the amounts and orientations of the tilt phase shifts are crucial to the performance of that method (for instance, if the amount of the tilt phase shift is less than 1λ, it is impossible to estimate the related tilt-shift parameters by that method), which may not be carefully controlled in the real testing environments. In the papers [1, 15] each interferogram is processed with the Fourier transform method (spatial carriers are assumed in the interferograms) so as to solve the related tilt phase shift parameters, and then the modulation phase can be extracted with the least-squares method. The paper [1] even has implemented experiments using a Mach-Zehnder interferometer, where the phase-shifts are introduced by induced vibrations without the need of any phase-shifter device. As the paper [1] includes affine registration preprocessing of the interferograms regarding for the camera movement due to strong vibration, it is expected to be applicable in a very hostile testing environment. However, if the interferograms are composed of closed fringes or the carrier frequency is inadequate, the methods given in the papers [1, 15] will be erroneous.

In this paper a novel method is proposed to cope with unknown tilt phase shifts. Firstly, a set of interferograms are captured, and the coarse wrapped phase shift distributions across the whole field between the selected interferograms are estimated in order, based on the regularized optical flow method [10]. Secondly, those phase shift distributions are subsequently filtered and unwrapped, the parameters related to the tilt phase-shifting planes are determined with the help of an efficient 2-D frequency estimation method [16], and the sign ambiguity of the tilt phase shifts among interferograms is eliminated by comparing the coarse phase demodulation results obtained with the regularized optical flow method. Finally, accurate demodulation phase can be extracted from all the interferograms by the least-squares method, with all the estimated tilt phase shifts. The proposed method has some more flexibility than the previous methods [1, 1115], i.e., it not only has relatively little restrictions on the amounts or orientations of the tilt phase shifts, but also works well with interferograms including open and closed fringes in any combination.

It should be pointed out that the extracted phase by the proposed method would have indetermination in the global sign if no prior knowledge about the phase shifts is available. However, such indetermination is not a particular problem of our method, because it is common to any asynchronous approach, if no information is given about the phase-shifts [8].

The rest of the paper is organized as follows: the principle of the proposed method is provided in detail in section 2, simulation and experimental results are reported in section 3 and section 4, respectively. Finally, conclusions are drawn in Section 5.

2. Principle

As for the set of interferograms with unknown tilt phase shifts, the intensity located at pixel(x,y)can be represented as [11, 13]:

In(x,y)=A(x,y)+B(x,y)cos[φ(x,y)+δn(x,y)]=A(x,y)+B(x,y)cos[φ(x,y)+(kxnx+kyny+dn)]
where A(x,y) is the background, B(x,y) is the modulation of the fringe pattern, φ(x,y)is the phase to be determined, and δn(x,y) represents the phase shift related to thenth interferogram. δn(x,y) is located on a so-called phase-shifting plane, which can be described with the tilt-shift parameterskxn,kyn,anddn. Without loss of generality, we can defineδ1(x,y)0. In principle, once the phase shiftsδn(x,y) are known, all the unknowns includingA(x,y),B(x,y), φ(x,y)can be uniquely determined with a minimum of three interferograms [13]. For instance, with the interferogramsI1,I2,I3, the wrapped phase φ can be solved as [13]:
φ=arctanI3I2+(I1I3)cosδ2+(I2I1)cosδ3(I1I3)sinδ2+(I2I1)sinδ3
Here the pixel coordinates have been omitted for simplicity.

In practice, if possible, it is always preferable to take more than three interferograms for the phase extraction. The reasons for it can be easily understood. On one hand, for pixels where one or both of the phase shifts are equal or close to integer numbers of 2π, Eq. (2) will be invalid and the corresponding phase extraction result will be untruthful, as the equivalent number of interferograms is reduced for those pixels; On the other hand, as a rule of thumb, the PSI methods’ resistance ability to noise can be improved as more interferograms are involved.

Figure 1 shows the flow chart of the proposed method in this paper.

 figure: Fig. 1

Fig. 1 The flow chart of the proposed method.

Download Full Size | PDF

Obviously, the accurate estimation of all the phase-shifting plane parameters {kxn,kyn,dn;n=2,3,...N}[see Eq. (1)] would be crucial to the performance of the phase extraction procedure. As we will see below, the proposed method would provide an accurate and fast solution for it.

Next, we will first briefly review the regularized optical flow method for the two-step interferometry, and then further introduce the procedures of the proposed method step by step. To be simple and intuitive, the descriptions in sections 2.1-2.4 are based on the case of two interferograms.

2.1 The regularized optical flow method for the two-step interferometry

Supposing there are two interferogramsI1,I2, and their high-pass filtered background-suppressed versions are denoted as I˜1 and I˜2. Then we have:I˜1B(x,y)cos[φ(x,y)], I˜2=B(x,y)cos[φ(x,y)+δ(x,y)], where δ(x,y)represents the phase shift distribution between the I1 and I2, which is located at the phase-shifting plane: δ(x,y)=kxx+kyy+d.

The paper [10] has put forward a regularized optical flow method to extract the wrapped demodulation phase from two interferograms with an arbitrary unknown but constant phase shift inside the range(0,2π), with a single exception of δ=π, i.e. kx=ky=0 and dπ are assumed. First, the fringe direction is obtained with a regularized optical flow method in an iterative way

uk+1=u¯kIx[Ixu¯k+Iyv¯k+It]/(ρ2+Ix2+Iy2)vk+1=v¯kIy[Ixu¯k+Iyv¯k+It]/(ρ2+Ix2+Iy2)
η=arctan(ν/u)
where Ixand Iy represent the derivatives of I˜1with respect to x,y, while Itrepresents the difference between the two background-suppressed interferograms, i.e. It=I˜2I˜1; uk+1 and vk+1are the velocity components obtained in the iteration k+1, whileu¯kandv¯k correspond to the mean value of uand v in a defined neighborhood; ρis the regularizing parameter that weighs the smoothness of uand v;ηrepresents the fringe direction map. Then, the modulation phase can be extracted using the spiral phase transform
SPT{}=FT1{(ωx+iωyωx2+ωy2)FT{}}
φest=arctan(iexp(iη)SPT{I˜1}I˜1)
where SPT{}denotes the spiral phase transform; ωxand ωyare the coordinates in the spectral domain; φest is the extracted phase by the regularized optical flow method.

2.2 Computation of the wrapped phase shift distribution between two interferograms

Actually, the phase extraction method given in the paper [10] can provide us with more information than the authors revealed.

On one hand, the phase shift between the two interferograms can be easily obtained with minor extension to that method, if its distribution is uniform across the field (i.e.kx=ky=0). We find that the shape of the direction map η is almost invariant as the phase shift changes, but the global sign of it is dependent on the range value of the phase shiftδ. Specifically, if we denote the resultant direction map as η1(δ(0,π)), and η2 (δ(π,2π)), respectively, then we will haveη2η1. As a result, we further have

φest(x,y){mod(φ(x,y),2π),0<δ<πmod(φ(x,y),2π),π<δ<2π
Here, if we suppose0<δ<π, and exchange the order of the two interferograms, i.e. letting
I˜1'(x,y)=I˜2(x,y)B(x,y)cos[φ(x,y)]I˜2'(x,y)=I˜1(x,y)cos[φ(x,y)δ]=cos[φ(x,y)+δ]
whereφ(x,y)=φ(x,y)+δ, andδ=2πδ(π,2π); mod()is a function related to the modulo operation. Then we can obtain another estimation result φest(x,y) of the demodulation phase with Eqs. (3)-(6). Meanwhile, from previous analysis, it is easy to deduce that: φest(x,y)mod(φ(x,y)+δ,2π). Thus we can further get the following equation
δest(x,y)=mod(φest(x,y)+φest(x,y),2π)2πδ
Similarly, if we supposeπ<δ<2π, we will have
δest(x,y)=mod(φest(x,y)+φest(x,y),2π)δ
From Eqs. (9) and (10), we can conclude that by adding the two extracted phase results obtained with the regularized optical flow method, an estimation of the phase shift can be obtained which is specially wrapped inside the approximate range (π,2π). The accuracy of the phase shift estimation can be improved by averaging it across the whole field, i.e., lettingδest=mean[δest(x,y)], where the functionmean()is related to the average operation. Additionally, at some pixels the estimated phase shiftφest(x,y)will fall outside the range(π,2π), mainly due to the noise in the interferograms and the truncation effect of the modulo operation. To compensate for it, we will add 2πto these pixels where φest(x,y)(0,π/2).

On the other hand, the phase shift distribution can also be analyzed as described above even if it is non-uniform across the whole field, since the regularized optical flow method is implemented in a local way. However, if the phase shift is equal to kπ(kZ)somewhere in the field, the computed wrapped phase shift distribution δest(x,y)would include some flip lines. In this case, it should be correctly unwrapped inside the approximate range(0,2π), which can be implemented as follows.

2.3 Unwrap the phase shift distribution

First, the map δest(x,y)is filtered with a mean filter (the size of filter adopted in this paper is 10×10pixels) to generate a map denoted asδ¯est(x,y), then a binary map denoted asBM1can be further produced as

BM1(x,y)={1,|δ¯est(x,y)π|π/20,|δ¯est(x,y)π|>π/2
Subsequently, the connected-component labeling and the morphological operations including erosion and dilation [17, 18] are successively applied, to remove the spurious pixels and to fill small holes in the mapBM1. Then the map δ¯est(x,y) can be divided into several sub-regions according to the value of BM1(x,y). In each sub-region of the map δ¯est(x,y)the extreme-value points are searched by scanning along the row or column direction (the selection of scanning direction is dependent on the orientation of the sub-regions), and then the corresponding flip line can be fixed by linear fitting with the coordinates of the extreme-value points. The estimated flip line related to the border sub-region would be rechecked to ensure its validity, by referring to the fitting residual error and comparing the slope of it with the counterparts in other sub-regions. Supposing totally nestimated flip lines are found, the whole field will be divided inton+1sub-regions, which will be labeled assubi(i=1,2,,n+1)in order by their spatial locations. Then a new binary map BM2 can be produced by letting
BM2(x,y)={1,p(x,y)sub2m,m=1,2,...,[(n+1)/2]0,p(x,y)sub2m,m=1,2,...,[(n+1)/2]
or
BM2(x,y)={1,p(x,y)sub2m,m=1,2,...,[(n+1)/2]0,p(x,y)sub2m,m=1,2,...,[(n+1)/2]
where p(x,y)denotes the pixel point (x,y). Then, the mapδ¯est(x,y) can be unwrapped inside the approximate range (0,2π) as

δ˜est(x,y)={δ¯est(x,y),BM2(x,y)=12πδ¯est(x,y),BM2(x,y)=0

2.4 Estimation of the tilt phase shift parameters

From the expression δ(x,y)=kxx+kyy+drelated to the phase-shifting plane, we can deduce that the data in any row (column) of the complex signal exp[jδ(x,y)]is single-tone, with frequency kxrad/pixel(kyrad/pixel). Thus, we can estimate the phase-shifting plane parameterskx,kyfrom exp[jδ˜est(x,y)] using an efficient 2-D frequency estimation method given in the paper [16]. If the estimated results are denoted as kx,estandky,est, then the estimated parameter d can be further obtained as

dest=angle{x,yexp[jδ˜est(x,y)]exp[j(kx,estx+ky,esty)]}
where angle{}denotes the operation of extracting the phase angle from a complex number. Then the estimated phase-shifting plane between the interferograms I1,I2can be obtained as

δ˜2(x,y)=kx,estx+ky,esty+dest

2.5 Elimination of the sign ambiguity among the estimated phase-shifting planes

Supposing there are totallyNphase-shifted interferograms, then all the phase-shifting planesδi(x,y),2iN can be estimated in order with the procedures introduced in section 2.1-2.4. However, it can be noted that as to the specific phase-shifting plane, the selection of Eq. (12) or Eq. (13), will give rise to two different estimation results with the same phase-shifting distribution but with an opposite global sign. As a result, the estimated phase-shifting planes will have at most2N1possible arrangements {±δ˜2(x,y),±δ˜3(x,y),,±δ˜N(x,y)}, and the majority of them will give erroneous phase extraction results, i.e., only the two arrangement which are closest to the true phase-shifting arrangement {δ2(x,y),δ3(x,y),,δN(x,y)} or its opposite version {δ2(x,y),δ3(x,y),,δN(x,y)}, will achieve the accurate estimation of the modulation phase φ(x,y) or its opposite version φ(x,y). Here we put forward a simple way to find the interested arrangement. Considering the case of three phase-shifted interferograms Ii,Ij,Ik(1i,j,kN), where the estimated coarse phase distributions [see Eq. (6)], the intermediate binary maps [see Eqs. (12) or (13)], and the estimated phase-shifting planes from the two-step phase-shifted interferograms {Ii,Ij}, {Ii,Ik}are denoted as{φest,ij(x,y),φest,ik(x,y)}, {BM2ij,BM2ik}, and {δ˜ij(x,y),δ˜ik(x,y)}, respectively. Then we can define

φ˜est,ij(x,y)={φest,ij(x,y),BM2ij(x,y)=12πφest,ij(x,y),BM2ij(x,y)=0
The φ˜est,ik(x,y)is defined withφest,ik(x,y)and the mapBM2ikin the similar way. Subsequently, we further define
val1=|x,y{exp[jφ˜est,ij(x,y)]exp[jφ˜est,ik(x,y)]}|val2=|x,y{exp[jφ˜est,ij(x,y)]exp[jφ˜est,ik(x,y)]}|
if val1>val2, the arrangements {δ˜ij(x,y),δ˜ik(x,y)} or {δ˜ij(x,y),δ˜ik(x,y)}will be adopted, otherwise the arrangements {δ˜ij(x,y),δ˜ik(x,y)}or {δ˜ij(x,y),δ˜ik(x,y)}will be adopted. When there are more than three phase-shifted interferograms, the phase-shifting arrangement can be searched in the similar way.

2.6 Phase extraction with the least-squares method

Equation (1) can be rewritten as

In(x,y)=a(x,y)+b(x,y)cos[δn(x,y)]+c(x,y)sin[δn(x,y)]
Where A(x,y)=a(x,y),b(x,y)=B(x,y)cos[φ(x,y)],c(x,y)=B(x,y)sin[φ(x,y)], δn(x,y)=kxnx+kyny+dn. Supposing there are totallyNphase-shifted interferograms (N3), after all the phase-shifting planesδest,n(x,y)have been determined, the phase distributionφ(x,y)can be extracted with the least-squares method [1, 4, 15] as follows
[aest(x,y)best(x,y)cest(x,y)]=[Nn=1Ncos[δest,n(x,y)]n=1Nsin[δest,n(x,y)]n=1Ncos[δest,n(x,y)]n=1Ncos2[δest,n(x,y)]n=1Ncsn(x,y)n=1Nsin[δest,n(x,y)]n=1Ncsn(x,y)n=1Nsin2[δest,n(x,y)]]1×[n=1NIn(x,y)n=1NIn(x,y)cos[δest,n(x,y)]n=1NIn(x,y)sin[δest,n(x,y)]]
φ˜(x,y)=arctan[cest(x,y)/best(x,y)]
where in Eq. (20), csn(x,y)=cos[δn(x,y)]sin[δn(x,y)]. φ˜(x,y)is the phase extraction result; while the estimated contrast parametersB˜(x,y)=sqrt(cest2+best2)are dependent both onB(x,y)and the phase shifts related to the specific pixel(x,y), which can be taken as the reliability measure of the phase extraction result at that pixel.

3. Simulation results

As the accurate estimation of the phase-shifting plane parameters is crucial to the performance of the phase extraction procedure, a series of computation simulations have been carried out to verify the effectiveness of phase-shifting plane estimation by the proposed method. Some results will be given in the following paragraphs.

In the first simulation, the expressions for the two background-suppressed inteferograms are I˜1=B(x,y)cos[φ(x,y)]+n1(x,y)and I˜2=B(x,y)cos[φ(x,y)+δ2(x,y)]+n2(x,y), where φ(x,y)=2×peaks(256)+0.0586π×(x+y), δ2(x,y)=0.0179x+0.0164y+1.583, and B(x,y)=60exp(0.0052x20.0052y2); the units of φ(x,y)and δ2(x,y)are both in radians, and the image sizes of I˜1and I˜2are both 256×256, i.e., xandy[127,126,,128]; n1(x,y)andn2(x,y)represent zero-mean additive Gaussian white noise, the standard deviation of which are both equal to 6; while peaks()is a built-in function of MATLAB, which is obtained by translating and scaling Gaussian distributions. The simulated phase-shifting plane δ2(x,y)is shown in Fig. 2(c), the tilt-shift amount across the whole field is 1.40λ(8.78rad); the estimated phase-shifting plane δ˜2(x,y)is shown in Fig. 2(i), and the residual PV (peak to valley) error and RMS (root mean square) error are as small as 0.019λ(0.012rad) and 0.004λ(0.026rad), respectively [see Fig. 2(j)]. The total processing time of this simulation is1.32susing a2.5GHzlaptop and MATLAB, and the main allocations of the runtime are as follows: about0.50sis spent associated with the regularized optical flow method (that method is run two times to estimate the phase shift); about0.54sis spent associated with the spatial imaging processing of the intermediate binary maps BM1and BM2, including the connected-component labeling, the morphological operations, and division of sub-regions operations; and about 0.04sis spent on solving the parameters kx,ky by the frequency estimation method. We noted that the residual error distribution of the phase-shifting plane is linear proportional to the estimation errors of the parameterskx,kyand d, thus the residual error will accumulate to its highest values at the border region of the field. If we divide the field into 2×2blocks and solve the parameter d[see Eq. (15)] locally for each block, then the resultant residual PV error and RMS error related to this simulation will be decreased by 33% and 38%, respectively. However, such improvement ratios are case-dependent, so we will not have further discussions on it.

 figure: Fig. 2

Fig. 2 (a), (b): The simulated background-suppressed inteferograms I˜1andI˜2; (c) The simulated phase-shifting plane δ2(x,y); (d) The mapδ¯est(x,y); (e),(f): The binary mapBM1via Eq. (11) before and after being processed with the connected-component labeling as well as the morphological operations; (g) The binary mapBM2; (h) The mapδ˜est(x,y)computed via Eq. (14); (i) The estimated phase-shifting planeδ˜2(x,y); (j) The residual error of the estimated phase-shifting plane, i.e. the wrapped difference between (c) and (i). The data shown in (c), (d), (h)-(j) are all in radians. In (e)-(g), the black and gray color represent the values of zero and one, respectively.

Download Full Size | PDF

We also have evaluated the performance of our method by testing some other phase-shifting planes with different tilt-shift amounts and different orientations. The simulation parameters are all the same as in Fig. 2, except for the simulated phase-shifting plane. Specifically, four different phase-shifting planes are considered, whereδ3(x,y)=δ2(x,y)/2,δ4(x,y)=2δ2(x,y), δ5(x,y)=0.0011x+0.0243y+1.583, andδ6(x,y)=0.0164x+0.0179y+1.583, i.e., the corresponding tilt-shift amounts across the field are 0.70λ,2.80λ,1.40λand1.40λ, respectively. Actually, δ5(x,y) and δ6(x,y)are obtained by rotating δ2(x,y)anticlockwise in the x-y coordinate system with 0.25πand0.5πrad, respectively. Herein, the definition ofδ2(x,y)is the same as shown in Fig. 2(c). The residual errors of the estimation results in PV and RMS measures are shown in Table 1. As can be seen from Table 1, all the residual PV errors and RMS errors are within0.016λand0.004λ, respectively.

Tables Icon

Table 1. Residual Errors of the Estimated Phase-shifting Planes

The above simulations demonstrate that the proposed method can well estimated the tilt phase-shifting planes, and the estimation performance is robust to noise, the orientation and the amplitude of the phase-shifting plane. As explained before, the phase-shifting plane estimation result will face a global sign ambiguity problem. To reasonably evaluate the estimation performance of the proposed method, the simulation results given in Fig. 1 and Table 1 are all assigned with the correct signs.

4. Experimental results

For further verification of the feasibility of the proposed method, we have applied it to the experimental interferograms. Two set of experimental results will be provided below.

In the first experiment, the data of interferograms are obtained from the paper [10], which include totally nine real phase-shifted interferograms with unknown but constant phase shifts (the image size of them is481×641). In this case, we only need to estimate the parametersdi. In Fig. 3 we show the reference phase obtained by the AIA method [4] using all the nine interferograms (b), as well as the phase extraction results by the AIA method but with only eight interferograms (c) and the proposed method (d). As the extracted phase by the AIA method would include a trivial constant bias, the residual errors of the latter two results after removing the bias would be0.05λin PV measure, 0.006λin RMS measure [see Fig. 3(e)], and0.024λin PV measure, 0.007λin RMS measure [see Fig. 3(f)], respectively. Particularly, we think the residual errors shown in Fig. 3(e) can be taken as an indirect measure of the non-ideality of the testing environment. Then we can conclude that in this experiment the proposed method has similar accuracy to the AIA method, as the residual errors shown in Fig. 3(e) and Fig. 3(f) are comparable in PV and RMS measures. The typical “double-frequency” fringe error shown in Fig. 3(f) is due to the discrepancies in phase shift estimations between the AIA method and the proposed method, i.e., as to the AIA method, the phase shifts are estimated in an iterative way and the information of different interferograms will be coupled in the least-squares sense during the iteration process, while for the proposed method, the estimation of phase shift between two interferograms is independent of the information encoded in other interferograms.

 figure: Fig. 3

Fig. 3 (a) The first interferogram used in this experiment; (b) The wrapped reference phase map by the AIA method using all the interferograms; (c) The wrapped phase extraction result by the AIA method with eight interferograms (i.e., one of the interferograms is excluded); (d) The wrapped phase extraction result by the proposed method using all the interferograms; (e) The wrapped phase difference between (b) and (c), after removing the bias; (f) The wrapped phase difference between (b) and (d), after removing the bias. The data shown in (b)-(f) are in radians.

Download Full Size | PDF

In another experiment, a circular flat-based diffractive optical element (DOE) with aperture of 100mm is tested, by measuring its transmitted phase. This DOE has quasi-continuous surface fabricated by the ion-beam etching technology, and it was made on the substrate of K9 glass. Firstly, the measurements are performed with a ZYGO Fizeau interferometer with the active vibration isolation workstation turned off. The recorded interferograms include moderate amount of tilt phase shift induced by environmental vibration and considerably large tilt phase shift by purposely rotating the reflective mirror (located behind the DOE) in the test arm around the axis of the PZT. As the rotation operation is implemented manually, in our experiment the interferograms are captured with a time interval of about 3s, to eliminate the contrast reduction in them due to the sudden “rotation” actions. We randomly pick out six consecutively captured frames to extract the phase component, where two of the frames are shown in Figs. 4 (a) and 4(b), and the wrapped phase extraction result with the proposed method is given in Fig. 4(c). On the other hand, we have also tested the DOE by the Zygo interferometer under vibration-isolated conditions with a calibrated PZT. The phase of the test surface extracted by the Zygo’s MetroPro software is shown in Fig. 4(d). Since noticeable lighting scattering would take place at the staircase locations of the DOE surface, the phase data are missing somewhere with the Zygo’s software. The difference between the phase extraction results by the proposed method and the Zygo’s software is presented in Fig. 4(e), which amounts to0.028λin RMS measures with the piston and tilt components removed; it must be pointed out that from the sense of surface test the relative difference between them would decrease nearly by half, as all the phase extraction results given in the following figures relate to the “double-pass” transmitted wave-front. In addition, the estimated results of the phase-shifting planes with the proposed method are shown in Fig. 5. It can be seen that the tilt-shift amounts between consecutive interferograms reach4.15λ,3.80λ,2.33λ,4.07λ,6.18λacross the field, respectively; as to the accumulated phase shifts relative to the first interferogram, the highest tilt-shift amount across the field would be 14.34λ, corresponding to the phase shift distribution between the fifth and the first interferograms.

 figure: Fig. 4

Fig. 4 (a),(b) The first two interferogram used in this experiment with a image size of 240×240, corresponding to a part of the DOE area; (c) The wrapped phase extraction result by the proposed method; (d) The phase extraction result by the Zygo’s MetroPro software with the piston and tilt components removed, where the phase data related to the pixels in white color are missing; (e) The wrapped phase difference between (c) and (d), with the piston and tilt components removed; (f) The normalized estimated contrast parametersB˜(x,y)shown in thelog10scale (the definition of B˜(x,y)can be found in section 2.6). The data shown in (c)-(e) are in radians.

Download Full Size | PDF

 figure: Fig. 5

Fig. 5 The estimated phase-shifting planes: (a) between the first and the second interferograms; (b) between the second and the third interferograms; (c) between the third and the fourth interferograms; (d) between the fourth and the fifth interferograms; (e) between the fifth and the sixth interferograms. All the data shown in these figures are in radians.

Download Full Size | PDF

The residual error as shown in Fig. 4 is mainly due to the following factors: (1) the retrace errors; (2) the errors accompanied with the low estimated contrast parameters; (3) the fluctuations of the laser power; and (4) the dynamic variations of the testing environment, such as the air turbulence, which will be discussed in more detail. Firstly, as the wave-front of the DOE is far from the ideal plane wave, the introduction of tilt phase shifts will also give rise to non-negligible retrace errors. To demonstrate it, the DOE is measured at different status under vibration-isolated conditions. As shown in Fig. 6(c), the difference between the phase extraction results is mainly composed of tilt component, which equals to7.25λ in the PV measure, and the resultant relative retrace error between them is shown in Fig. 6(d), which equals to0.023λin the RMS measure. As the tilt shift amount increases, the retrace error will be larger, and vice versa. As to our experiment, different retrace errors have been introduced in the interferograms we used, and the phase extraction result by the Zygo’s software also include some retrace error itself. Therefore, the residual error as shown in Fig. 4(e) will inevitably contain the contribution of retrace errors. Secondly, as shown in Fig. 4(e) and (f) the areas with high residual error values are found to be related with extremely low estimated contrast parameters. It is easy to understand that in such areas the accuracy of phase extraction result is much more susceptible to the noise in the interferograms, as well as the residual estimation errors of the phase-shifting planes. This problem can be alleviated by using more interferograms for the phase extraction. Thirdly, as the interferograms used in this experiment were captured during a relative long time (about 16s), the parametersA(x,y),B(x,y)[see Eq. (1)] will be time dependent as a result of the fluctuations in laser power, which will introduce some error into the phase extraction result. Finally, the dynamic variations of the testing environment, such as the air turbulence will also make some differences between measurements.

 figure: Fig. 6

Fig. 6 Difference between phase extraction results by the Zygo’s MetroPro software at different testing status, to demonstrate the relative retrace error. (a) The typical interferogram related to the testing status 1; (b) The typical interferogram related to the testing status 2; (c) The difference in phase extraction results by the Zygo’s MetroPro software, between the testing status 1 and the testing status 2 (the phase data related to the pixels in white color are missing); (d) the same as (c), but with the piston and tilt components removed. The data shown in (c), (d) are in radians.

Download Full Size | PDF

As to the second experiment, since the tilt phase shift amounts related to the phase-shifting planes are considerably large, the AIA method will be failed; and the block-wise methods given in the papers [12, 14] also can hardly retrieve the correct phase, as in this case the interferograms should be divided into extremely many blocks with small sizes (probably more than one thousand blocks are required), so that the estimated phase shifts related to each small block are prone to be inaccurate. Besides, the methods given in the papers [1, 15] also can hardly be applied in this experiment, which rely on introduction of adequate spatial carriers into the interferograms. Additionally, we also found that if the tilt-shift amount across the field get too large (for instance, larger than10λ), the direct estimation of the phase-shifting plane by the proposed method is prone to be inaccurate, due to the probable failure in the correct sub-region division of the intermediate binary map BM2(as to the definition of the mapBM2, please referring to the section 2.3).

5. Conclusion

We have proposed a novel method for extracting phase distribution from interferograms with unknown tilt phase shifts. The proposed method can estimate the unknown tilt phase shift between two temporal phase-shifted interferograms, by extending the regularized optical flow method provided in the paper [10] with the spatial image processing and frequency estimation technology. With all the estimated tilt phase shifts, we can further obtain the phase extraction result with the least-squares method. Both simulation and experimental results have proved the feasibility of the proposed method. The proposed method is expected to be used in a testing environment with low frequency and high amplitude vibration, where costly and accurate phase-shifting devices are not longer required for steady-state measures.

Additionally, the proposed method has shown some more flexibility in comparison with the previous methods [1, 1115]. On one hand, the proposed method has relatively little restrictions on the amounts or orientations of the tilt phase shifts. The methods proposed in [11, 12, 14] are mainly adaptable to the phase-shifted interferograms with relatively low tilt-phase shift amounts (the typical amounts reported in those papers are much less than 1λ); and the method proposed in [13] is quite sensitive to the amounts and orientations of the unknown tilt phase shifts among interferograms, particularly, if the amount of the tilt phase shift is less than 1λ, that method will fail. While, the method proposed in this paper can handle with the unknown tilt phase shifts with amount up to several wavelength, no matter the orientations of them. On the other hand, unlike the methods proposed in the papers [1, 15], the method proposed in this paper requires no spatial carriers in the interferograms, i.e., it can work well with interferograms including open and closed fringes in any combination.

Acknowledgments

We are grateful to Professor Javier Vargas at Centro Nacional de Biotecnología-CSIC, Spain, for his kind help about the phase-shifting methods and providing us with partial experimental data of interferograms used in this paper. The research was partially supported by the National Basic Research Program of China under grant No. 2011CB706701. Thanks also go to the anonymous reviewers for their valuable comments and suggestions.

References and links

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

2. D. Malacara, Optical Shop Testing, 3rd ed. (Wiley, 2007), Chap. 14, pp. 547–666.

3. R. Juarez-Salazar, C. Robledo-Sánchez, C. Meneses-Fabian, F. Guerrero-Sánchez, and L. M. Arévalo Aguilar, “Generalized phase-shifting interferometry by parameter estimation with the least squares method,” Opt. Lasers Eng. 51(5), 626–632 (2013). [CrossRef]  

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

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

6. P. Gao, B. Yao, N. Lindlein, K. Mantel, I. Harder, and E. Geist, “Phase-shift extraction for generalized phase-shifting interferometry,” Opt. Lett. 34(22), 3553–3555 (2009). [CrossRef]   [PubMed]  

7. B. Li, L. Chen, W. Tuya, S. Ma, and R. Zhu, “Carrier squeezing interferometry: suppressing phase errors from the inaccurate phase shift,” Opt. Lett. 36(6), 996–998 (2011). [CrossRef]   [PubMed]  

8. J. Vargas and C. O. S. Sorzano, “Quadrature Component Analysis for interferometry,” Opt. Lasers Eng. 51(5), 637–641 (2013). [CrossRef]  

9. H. Guo, “Blind self-calibrating algorithm for phase-shifting interferometry by use of cross-bispectrum,” Opt. Express 19(8), 7807–7815 (2011). [CrossRef]   [PubMed]  

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

11. M. Chen, H. Guo, and C. Wei, “Algorithm immune to tilt phase-shifting error for phase-shifting interferometers,” Appl. Opt. 39(22), 3894–3898 (2000). [CrossRef]   [PubMed]  

12. A. Dobroiu, D. Apostol, V. Nascov, and V. Damian, “Tilt-compensating algorithm for phase-shift interferometry,” Appl. Opt. 41(13), 2435–2439 (2002). [CrossRef]   [PubMed]  

13. O. Soloviev and G. Vdovin, “Phase extraction from three and more interferograms registered with different unknown wavefront tilts,” Opt. Express 13(10), 3743–3753 (2005). [CrossRef]   [PubMed]  

14. J. Xu, Q. Xu, and L. Chai, “Iterative algorithm for phase extraction from interferograms with random and spatially nonuniform phase shifts,” Appl. Opt. 47(3), 480–485 (2008). [CrossRef]   [PubMed]  

15. J. Xu, Q. Xu, and L. Chai, “Tilt-shift determination and compensation in phase-shifting interferometry,” J. Opt. A, Pure Appl. Opt. 10(7), 075011 (2008). [CrossRef]  

16. S. Ye and E. Aboutanios, “Two dimensional frequency estimation by interpolation on Fourier coefficients,” in Proc. of IEEE Int. Conf. on Acoustics, Speech and Signal Processing, 3353–3356 (2012). [CrossRef]  

17. R. C. Gonzalez, R. E. Woods, and S. L. Eddins, Digital Image Processing Using MATLAB (Prentice Hall, 2004).

18. L. G. Shapiro and G. C. Stockman, Computer Vision (Prentice Hall, Upper Saddle River, New Jersey, USA, 2001).

Cited By

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

Alert me when this article is cited.


Figures (6)

Fig. 1
Fig. 1 The flow chart of the proposed method.
Fig. 2
Fig. 2 (a), (b): The simulated background-suppressed inteferograms I ˜ 1 and I ˜ 2 ; (c) The simulated phase-shifting plane δ 2 (x,y) ; (d) The map δ ¯ est (x,y) ; (e),(f): The binary map BM1 via Eq. (11) before and after being processed with the connected-component labeling as well as the morphological operations; (g) The binary map BM2 ; (h) The map δ ˜ est (x,y) computed via Eq. (14); (i) The estimated phase-shifting plane δ ˜ 2 (x,y) ; (j) The residual error of the estimated phase-shifting plane, i.e. the wrapped difference between (c) and (i). The data shown in (c), (d), (h)-(j) are all in radians. In (e)-(g), the black and gray color represent the values of zero and one, respectively.
Fig. 3
Fig. 3 (a) The first interferogram used in this experiment; (b) The wrapped reference phase map by the AIA method using all the interferograms; (c) The wrapped phase extraction result by the AIA method with eight interferograms (i.e., one of the interferograms is excluded); (d) The wrapped phase extraction result by the proposed method using all the interferograms; (e) The wrapped phase difference between (b) and (c), after removing the bias; (f) The wrapped phase difference between (b) and (d), after removing the bias. The data shown in (b)-(f) are in radians.
Fig. 4
Fig. 4 (a),(b) The first two interferogram used in this experiment with a image size of 240×240 , corresponding to a part of the DOE area; (c) The wrapped phase extraction result by the proposed method; (d) The phase extraction result by the Zygo’s MetroPro software with the piston and tilt components removed, where the phase data related to the pixels in white color are missing; (e) The wrapped phase difference between (c) and (d), with the piston and tilt components removed; (f) The normalized estimated contrast parameters B ˜ (x,y) shown in the log10 scale (the definition of B ˜ (x,y) can be found in section 2.6). The data shown in (c)-(e) are in radians.
Fig. 5
Fig. 5 The estimated phase-shifting planes: (a) between the first and the second interferograms; (b) between the second and the third interferograms; (c) between the third and the fourth interferograms; (d) between the fourth and the fifth interferograms; (e) between the fifth and the sixth interferograms. All the data shown in these figures are in radians.
Fig. 6
Fig. 6 Difference between phase extraction results by the Zygo’s MetroPro software at different testing status, to demonstrate the relative retrace error. (a) The typical interferogram related to the testing status 1; (b) The typical interferogram related to the testing status 2; (c) The difference in phase extraction results by the Zygo’s MetroPro software, between the testing status 1 and the testing status 2 (the phase data related to the pixels in white color are missing); (d) the same as (c), but with the piston and tilt components removed. The data shown in (c), (d) are in radians.

Tables (1)

Tables Icon

Table 1 Residual Errors of the Estimated Phase-shifting Planes

Equations (21)

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

I n ( x, y )=A( x, y )+B( x, y )cos[ φ( x, y )+ δ n ( x, y ) ] =A( x, y )+B( x, y )cos[ φ( x, y )+( k xn x+ k yn y+ d n ) ]
φ=arctan I 3 I 2 +( I 1 I 3 )cos δ 2 +( I 2 I 1 )cos δ 3 ( I 1 I 3 )sin δ 2 +( I 2 I 1 )sin δ 3
u k+1 = u ¯ k I x [ I x u ¯ k + I y v ¯ k + I t ] / ( ρ 2 + I x 2 + I y 2 ) v k+1 = v ¯ k I y [ I x u ¯ k + I y v ¯ k + I t ] / ( ρ 2 + I x 2 + I y 2 )
η=arctan( ν/u )
SPT{ }=F T 1 { ( ω x +i ω y ω x 2 + ω y 2 )FT{ } }
φ est =arctan( iexp( iη )SPT{ I ˜ 1 } I ˜ 1 )
φ est ( x,y ){ mod( φ( x,y ),2π ), 0<δ<π mod( φ( x,y ),2π ), π<δ<2π
I ˜ 1 '( x,y )= I ˜ 2 ( x,y )B(x,y)cos[ φ ( x,y ) ] I ˜ 2 '( x,y )= I ˜ 1 ( x,y )cos[ φ ( x,y )δ ]=cos[ φ ( x,y )+ δ ]
δ est ( x,y )=mod( φ est ( x,y )+ φ est ( x,y ),2π )2πδ
δ est ( x,y )=mod( φ est ( x,y )+ φ est ( x,y ),2π )δ
BM1( x,y )={ 1, | δ ¯ est ( x,y )π |π/2 0, | δ ¯ est ( x,y )π |>π/2
BM2( x,y )={ 1, p( x,y )su b 2m , m=1,2,...,[ ( n+1 ) /2 ] 0, p( x,y )su b 2m , m=1,2,...,[ ( n+1 ) /2 ]
BM2( x,y )={ 1, p( x,y )su b 2m , m=1,2,...,[ ( n+1 ) /2 ] 0, p( x,y )su b 2m , m=1,2,...,[ ( n+1 ) /2 ]
δ ˜ est ( x,y )={ δ ¯ est ( x,y ), BM2( x,y )=1 2π δ ¯ est ( x,y ), BM2( x,y )=0
d est =angle{ x,y exp[ j δ ˜ est ( x,y ) ]exp[ j( k x,est x+ k y,est y ) ] }
δ ˜ 2 (x, y)= k x,est x+ k y,est y+ d est
φ ˜ est,ij ( x,y )={ φ est,ij ( x,y ), BM 2 ij ( x,y )=1 2π φ est , ij ( x,y ), BM 2 ij ( x,y )=0
val1=| x,y { exp[ j φ ˜ est,ij ( x,y ) ]exp[ j φ ˜ est,ik ( x,y ) ] } | val2=| x,y { exp[ j φ ˜ est,ij ( x,y ) ]exp[ j φ ˜ est,ik ( x,y ) ] } |
I n ( x, y )=a( x, y )+b( x, y )cos[ δ n ( x, y ) ]+c( x, y )sin[ δ n ( x, y ) ]
[ a est ( x, y ) b est ( x, y ) c est ( x, y ) ]= [ N n=1 N cos[ δ est,n ( x, y ) ] n=1 N sin[ δ est,n ( x, y ) ] n=1 N cos[ δ est,n ( x, y ) ] n=1 N cos 2 [ δ est,n ( x, y ) ] n=1 N c s n ( x,y ) n=1 N sin[ δ est,n ( x, y ) ] n=1 N c s n ( x,y ) n=1 N sin 2 [ δ est,n ( x, y ) ] ] 1 ×[ n=1 N I n ( x, y ) n=1 N I n ( x, y )cos[ δ est,n ( x, y ) ] n=1 N I n ( x, y )sin[ δ est,n ( x, y ) ] ]
φ ˜ ( x, y )=arctan[ c est ( x,y ) / b est ( x,y ) ]
Select as filters


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