1. Introduction

The information about a physical quantity being measured interferometrically is contained in the phase term ϕ of the fringe pattern described by (1). The goal of analysis [1] consists in extracting this phase term. Many widely used phase-shifting methods are based on phase detection of a sinusoidal signal [2, 3, 4] and use at least three phase-shifted interferograms. The phase shifts, however, should be controlled with high accuracy. Any error in linearity of an actuator response can introduce an error into the calculated phase. A number of methods referred as self-calibrating techniques has been developed to deal with linear and nonlinear piston-shift errors and tilt-shifts errors. Lai and Yatagai [5] directly calculate phase shifts from additional linear fringes due to an additional tilted reference mirror moved together with the test object. In the methods suggested by Okada et al [6] and by Kong and Kim [7], the phase distribution and phase shifts are calculated alternatively by iteration of the least-square algorithm. Chen et al [8] have presented a modified method which uses the first-order Taylor series expansion of the phase-shifting errors to compensate also the phase error caused by tilt-shift errors. Dobroiu et al [9] cancel the effects of both phase shift errors and the presence of tilts by calculating average phase shifts for sufficiently small blocks of the interferogram space, so that the phase shifts inside every block can be assumed constant. Based on a contrast-level map, the algorithm iteratively adjusts the phase shifts and tilts.

Another popular approach is based on a spatial analysis of the interferogram, e.g. with the discrete Fourier transform [1, 10, 11, 12]. This method needs only one interferogram but works well only when the interferogram has a narrow bandwidth, a high-frequency carrier and a low noise level. The method cannot be used if the interferogram contains closed fringes, though different improvements have been proposed to overcome this obstacle [13]. Goldberg and Bokor [14] have proposed to use the Fourier-transform method together with the phase-shifting to determine unknown phase-step increments. The spatial analysis of the interferograms allows to find the phase steps deterministically, not iteratively, without any a priori information; while time-domain based phase shifting provides ease of implementation. Both phase-shifting and Fourier transform-based methods calculate the phase wrapped in the interval (-π,π]. Some methods — such as the local phase-tracking technique and genetic algorithms [15, 16, 17] — return unwrapped phase values, but these are computationally expensive.

This paper presents a new algorithm which allows one to reconstruct the phase from three interferograms registered with three different arbitrary unknown wavefront tilts. Unlike the methods described in [8, 9], the tilts are purposely introduced, and do not represent a phase-shifter error. On the contrary, the phase shifts are introduced in every interferogram point due to the tilts. The magnitude of the tilts is larger then in tilt-eliminating methods [8, 9] (several fringes per interferogram space), but not as large as in FFT-based methods [14]. Our method uses both the spatial analysis of the interferogram and the phase shifting. This makes it possible to alleviate usual for spatial- or time-domain-only methods restrictions related to the bandwidth, carrier frequency, noise level, background illumination and intensity variation, and, moreover, accuracy of the phase shifts. In fact, introduced phase shifts are calculated by the algorithm and thus need not be known a priori. The only requirement for these is that they change linearly over the detector region. Spatial analysis of two interferograms with such a phase difference will detect its parameters, and then a general phase-shifting algorithm can be used to retrieve the phase.

2. The algorithm

We will describe the interferogram as

where x=(x,y) denotes the position in the recorded image, and square brackets are used to emphasize its discrete nature. We assume the functions a and b to be dependent only on the detector position x. Sometimes we will omit the argument x for simplicity.

The method starts with recording three interferograms, I ₀, I ₁, I ₂, which differ from each other only by a small linear term in the phase

where τ_i[x]=t i·x are the tilt terms and σ_i the phase shifts. The interferograms can be obtained by, for instance, tilting the reference mirror in the interferometer by a small random angle. Again, we assume that a and b can be considered the same for all three interferograms.

For every pixel with the coordinates x, we have three equations (2), and we have five unknowns to find: a,b,ϕ0,δ₁=τ ₁+σ ₁,δ ₂=τ ₂+σ ₂. In phase shifting interferometry, δ ₁ and δ ₂ with τ ₁,τ ₂=0, are known, and equations (2) can be solved for every pixel to obtain the (wrapped) phase ϕ ₀. For instance, the following identity:

(I ₁-I ₂)cosϕ+(I ₂-I₀)cos(δ ₁+ϕ)+(I ₀-I ₁)cos(δ ₂+ϕ)=0

can be used to find the phase as

In our case, the phase shifts δ ₁ and δ ₂ are supposed to be unknown and different for every pixel. Moreover, for those pixels where one or both of the phase shifts are equal or close to integer numbers of 2π, the number of equations in system (2) is reduced from three to two, or one. However, as will be shown later, for the rest of the object area, the set of three interferograms contains enough information to determine the tilt and piston terms, and thus to make system (2) solvable.

The following method is proposed to find the unknown terms. Consider the differences of the interferogram I_d,1 and I_d,2:

I_d,1=I₁-I₀=b(cos(ϕ+τ₁+σ₁)-cos(ϕ))

and

I_d,2=I₂-I ₀=b(cos(ϕ+τ ₂+σ ₂)-cos(ϕ)),

which we can write as

Then, the set of zero-crossing points for each of the differences consists of two groups: the first one is determined by the term sin $\frac{{\tau }_i+{\sigma }_i}{2}$ when τ_i+σ_i=2kπ, and the second contains the points where $\varphi +\frac{{\tau }_i+{\sigma }_i}{2}=2k\pi$, where k is an integer. The first group is not dependent on a,b,ϕ and forms parallel lines. It can be fully described by three parameters:

where θ is the common normal, λ is the separation distance, and s is the distance from the origin. The algorithm proceeds by finding these parameters. Then t can be found as $\frac{2\pi }{\lambda }(\cos \theta ,\sin \theta )$ and $\sigma =\frac{2\pi s}{\lambda }$. This gives the phase shift values δ_i for every pixel, which can be substituted in equations (2) or directly in the phase formula (3).

 

Fig. 1. Three original interferograms I ₀, I ₁, I ₂ obtained with a Twyman-Green interferometer. The first interferogram corresponds to the original object; additional tilts were introduced in the second and third interferograms. Note that all three interferograms contain closed fringes.

Download Full Size | PDF

The next section will describe the details of the practical implementation of the algorithm and the method for extracting tilt and piston values from the set of zero-crossing points.

3. Practical implementation and examples

In this section we describe and illustrate one of possible implementations of the algorithm.

We begin with the three interferograms shown in Fig. 1. Subtracting the interferograms, we obtain the images shown in Fig. 2. Note the regions of middle gray, corresponding to zero differences.

To extract the zero-crossing lines, we used a threshold-based method in which the point is considered to be zero-crossing if the corresponding value is closer to zero than some chosen value. Fig. 3 shows the zero lines extracted in this step. Points extracted from each differential image form two sets: one with parallel lines, and one repeating the fringes of a “half-tilted” interferogram.

To determine the angle of parallel lines, their separation and shift relative to the origin, we used an algorithm based on the Hough transform [18]. The Hough transform maps any point (x,y) onto a sinusoidal function ρ=xcosθ+ysinθ. The key property of the Hough transform is that sinusoidals in Hough space associated with points lying on the same line have a common point of intersection, say (ρ ₀,θ ₀). This line has a normal (cosθ ₀, sinθ ₀) and is shifted by distance ρ ₀ from the origin.

In practice, the Hough space is limited to (ρ,θ)∊[-R,R]×[0,π) for some R (the maximum distance of any image point from the origin). Then it is quantized with finite steps ${\delta }_{\rho }=\frac{R}{m}$ and ${\delta }_{\theta }=\frac{\pi }{n}$ and represented by a matrix of size (2m+1)×n

A=(Ai,j), i=-m,…,m, j=0,n-1,

which is called an accumulator. The rows of A represent the angle θ ; the columns, the radius vector ρ. The sinusoidals for every feature pixel in Fig. 3 are “drawn” in the accumulator as follows. Initially, every cell of the accumulator is set to zero. Then, for every row number j, the column number i _(x,y)(j) is calculated by

$i_{(x,y)}\left(j\right)=\left[\frac{x\cos \left(j{\delta }_{\theta }\right)+y\sin \left(j{\delta }_{\theta }\right)}{{\delta }_{\rho }}-\frac{1}{2}\right],$,

where the square brackets denote the integer part. The element $A_{i_{(x,y)}\left(j\right),j}$ is then incremented by one. The result of the Hough transform is shown in Fig. 4. The intersections of a large number of sinusoidals are clearly visible.

 

Fig. 2. Differences of the interferograms I_d,1=I ₁-I ₀ and I _d,2=I ₂-I ₀. Black and white correspond to minimum (-255) and maximum (255) levels, respectively. The circle marks the mask. Masking the object edges helps to avoid spurious zero level points.

Download Full Size | PDF

 

Fig. 3. Points in interferogram differences which are close to zero are used as zero levels. The threshold was chosen to be 1, the minimal possible level that is greater then 0.

Download Full Size | PDF

After the Hough transform has been calculated, the parameters (5) are found as follows. Find the maximum element M in the A and use θ_M and ρ_M corresponding to its position as lines normal θ and shift from origin s. To find the separation distance λ, consider the row of the accumulator corresponding to θ_M (Fig. 5) and find the average distance between its local maxima, which have values greater than α_M for some threshold α (typical values are 0,3–0, 7, see Fig. 7). To reduce influence of the noise, the algorithm replaces the values of nearby positioned maxima with their average.

The maxima at the center can potentially reach higher values than maxima near the edges, due to the circular shape of the aperture. To compensate the effect, we can normalize the accumulator by dividing each cell by its maximum possible value (see Fig. 6). This makes maxima more uniform and facilitates locating them with the threshold method.

To test the values found for the added tilts one can divide the differential interferogram (4) by the first term $-2\sin \frac{{\tau }_i+{\sigma }_i}{2}$, to obtain

shown in Fig 8. With proper values for the tilt the picture should be as smooth as the original interferogram except in the vicinity of the parallel lines, where sin $\frac{{\tau }_i+{\sigma }_i}{2}$ is close to zero.

Now the phase ϕ can be found. Replacing the formula (3), one can use a simpler identity

sin(ϕ+α)=cos(ϕ+β) sin(α-β)+cos(α-β) sin(ϕ+β)

to write

$\cos \left(\varphi +\beta \right)=\frac{\sin \left(\varphi +\alpha \right)-\cos \left(\alpha -\beta \right)\sin \left(\varphi +\beta \right)}{\sin \left(\alpha -\beta \right)}.$.

Thus, with $\alpha ={\delta }_1⁄2=\frac{{\tau }_1+{\sigma }_1}{2}$ and $\beta ={\delta }_2⁄2=\frac{{\tau }_2+{\sigma }_2}{2}$,

This, together with (6b), allows us calculate the phase ϕ+δ₂/2 as

The result is presented in Fig. 9.

 

Fig. 4. The differential images I _d,1 and I _d,2 in the Hough space. R=354, δ_ρ=1, δ_θ=π/360. One can see that zero-level points of I _d,1 should contain about 8 lines at angle θ about 170π/360, and I _d,2 — about 14 lines at angle approximately equal to 70π/360.

Download Full Size | PDF

 

Fig. 5. Accumulators’ rows containing the maximum elements. The local maxima corresponds to the sinusoidal intersection points. The highest maxima corresponds to parallel lines in zero-crossing curves. The second raw of pictures shows the maxima values for normalized accumulator.

Download Full Size | PDF

 

Fig. 6. Normalized accumulators’ rows containing the maximum elements and a norm function we have used. The norm function is the same for every horizontal line in the Hough space.

Download Full Size | PDF

4. Discussion

4.1. Artifacts in the extracted phase

The extracted phase in Fig. 9 contains visible errors along the lines where either δ ₁, or δ ₂, or δ ₁–δ ₂ is equal to 2kπ, k∊ℤ, that is, where the system (2) is badly defined. Although these artifacts introduce a small rms error, they can seriously affect the unwrapped phase. They can be removed either by introducing a fourth interferogram with tilt δ ₃ in the algorithm and calculating the phase as the median of the results of three possible tilt combinations (see Fig. 10), or by obtaining a and b from the system (2) and then calculating the phase ϕ by substituting their smoothed by low-pass filter values in equation (1) or just by using robust and noise-immune phase unwrapping algorithms.

To illustrate the error caused by the artifacts, we used the phase extracted with the Fourier transform method from an interferogram with a large tilt (Fig. 11). Note the good correspondence between the two extracted phases.

 

Fig. 7. Lines detected by the Hough transform. The threshold α for maxima selection was set to 0.45 for the first set of zero-level points and to 0.25 for the second one.

Download Full Size | PDF

 

Fig. 8. Interferogram differences divided by $-2\sin \frac{{\tau }_i+{\sigma }_i}{2}$.

Download Full Size | PDF

4.2. Initial phase

If the phase originally has a large linear carrier, its fringes are close to the set of parallel lines themselves, and this can interfere with phase shift detection. In this case some initial estimates on introduced tilts are needed.

4.3. Optimal tilt range

As the tilt information is derived from spatial analysis based on periodicity, the linear term should change at least several wavelengths over the detector area. The more lines are contained in the zero-level points, the more accurate are the calculated parameters (5). For very small tilt values, the lines in the set of zero-level points will become too thick, even for minimal threshold values (Fig. 12), and the Hough transform can give a wrong value for the line slope. On the other hand, for large tilts, the zero-level points due to the second terms in (4) begin to approach parallel lines with the same slope, making tilt detection difficult. Thus the algorithm works best for those tilt values for which Fourier-transform based methods fail.

 

Fig. 9. Extracted phase $\varphi +\frac{{\tau }_2+{\sigma }_2}{2}$.

Download Full Size | PDF

4.4. Ambiguity in the detection of phase-shift signs

The Hough transform detects the normal θ only in the range from 0 to π. It is clear from the equations (4), that four different pairs of phase shifts {±δ ₁,±δ ₂} produce parallel lines with identical parameters (5). Thus the algorithm does not see the difference between interferograms with “positive” and “negative” tilts, and formula (8) provides, generally speaking, four different phase distributions. Two of them are the “real” phase accurate to the sign, and two the “faulty” phase, which can be recognized by its almost binary shape (see Fig. 13). The algorithm cannot determine the sign of the real phase, as the interferograms (2) are the same for ϕ ₀,δ ₁,δ ₂ and -ϕ ₀,-δ ₁,-δ ₂, and some additional information is necessary. In practice, one can limit all possible tilts only to those increasing in y-direction, which correspond to θ∊(0,π) and exclude the tilt signs ambiguity.

4.5. Computational effectiveness and accumulator resolution

The accuracy of the extracted tilts depends on the angle and radius resolution of the accumulator. To fill the accumulator with angle step δ_θ=π/n, one needs O(n×N) operations, where N is the number of points in the zero-level set. To speed up the calculation one can first estimate the tilt with low resolution, and then fill only the small region of the accumulator near the maximum with high resolution. Note that too high resolution in ρ makes the algorithm too sensitive to noise and to quantization of the interferograms. For the zero-level points from Fig. 12, for instance, low resolution in ρ and high resolution in θ is necessary.

 

Fig. 10. The fourth and fifth interferograms, the phase calculated for other tilt values, and the resulting phase obtained by median.

Download Full Size | PDF

 

Fig. 11. An interferogram with a linear carrier, the phase extracted with Fourier-transform method used as a reference, and wrapped difference between the reference and the median phase from Fig. 10.

Download Full Size | PDF

4.6. Robustness

The algorithm does not depend on a and b as long as they are the same for all the interferograms. In regions where b is close to zero or in cases with high noise level of the detector, spurious zero-level points appear and decrease the accuracy. The algorithm also fails on overexposed interferograms or on interferograms recorded with a non-linear camera.

5. Conclusions

We have proposed a new algorithm for phase retrieval from three interferograms that have different unknown tilt terms. The method is independent of the carrier frequency and the exact control of the phase shifts. The method can be used in applications which use inexpensive or non-linear actuators or a control system with such characteristics.

 

Fig. 12. A small tilt will give a wrong slope.

Download Full Size | PDF

 

Fig. 13. If one of the tilt signs is incorrect, formula (8) produces unrealistic phase map.

Download Full Size | PDF