1. Introduction

In optical measurement, retrieving the phase map from a single fringe pattern enables measuring dynamic phenomena in real-time [1–3], or making measurement results immune from the influences of circumstantial disturbances and vibrations. In practice, however, analysis of a single fringe pattern usually suffers from a difficulty in isolating fringes from background intensities. Thanks to that, Takeda et al. [4] proposed in the early 1980s introducing a linear carrier into a fringe pattern so that its fringes may have much higher frequencies, gradients, or variations than the background intensities, thus providing a possibility of separating the fringes from background intensities mathematically. Since then, techniques of spatial carrier fringe pattern analysis have been extensively studied.

The mainstay technique for analyzing a spatial carrier fringe pattern is the Fourier transform (FT) method [4–10]. In the Fourier spectrum of the pattern, the fringe component corresponds to a pair of conjugately symmetric lobes, which locate at carrier frequency positions, and are disjointed from that of the background. Isolating one of the two lobes using a filter, taking its inverse FT (IFT), and then implementing arctangent operation yield the phase map of the pattern. Because the FT method is a global method, it has a limitation that the large deformation in the phase map increases the bandwidth of fringe component and thus may make the filtering process unreliable in separating the component of interest [7]. To solve this problem, the space-frequency analysis methods (i.e. the time-frequency analysis methods for temporal signals), which can provide spatial information and spectral information simultaneously, are usually employed. Among them, the windowed Fourier transform (WFT) method [11] determines spatial localization by sliding a window with constant size over the fringe pattern, and therefore it flaws in having a fixed resolution. By using scaled wavelets, the wavelet transform (WT) method [12–15] possesses the ability of multi-resolution analysis. It gives high spatial resolution for high-frequency components, and high frequency resolution for low-frequency components. As an extension of WFT or the continuous WT, the S-transform (ST) method [16] uses a varying Gaussian window with its width and height depending on frequency, and hence enables yielding a better signal clarity. Although some efficient algorithms such as the fast Fourier transform (FFT) and the fast wavelet transform (FWT) have been invented, the transform-based methods generally have relatively high computational complexities.

Compared with processing in transform domains, a more straightforward solution is to process the fringe pattern in the spatial domain directly. According to the theory about the linear time-invariant (LTI) system, the filtering process in FT method can be equivalently implemented in the spatial domain by convoluting the fringe pattern with a complex kernel. Based on this principle, typically the spatial synchronous detection methods [17–19] and some classical spatial carrier phase-shifting (SCPS) methods [20–28] have been developed. With them, convolutions of the fringe pattern with the real and imaginary parts of the kernels are implemented individually, so that their formulas are generally similar to the standard temporal phase shifting formulas with known and fixed phase shifts. These spatial LTI methods usually use convoluting kernels with very small sizes in order to improve computational efficiencies or to take advantage of hardware implementations, at the expense of low frequency resolutions. With them, detuning errors induced by frequency mismatching may occur, especially when the phase map has a large deformation. Applying non-LTI systems to the fringe pattern is helpful for solving this problem. For example, assuming the local phase variation to be linear within a small region allows deducing SCPS formulas insensitive to linear phase errors [29, 30]; Employing a least-squares algorithm enables us estimating the local spatial frequencies of a pixel from its neighborhood and further calculating its phase by use of a two-dimensional (2D) SCPS algorithm [31]; Using an iterative strategy can improve the measurement accuracy of phase map by updating the phase shifts [32]; and in [33], the principal component analysis (PCA) as an effective tool for revealing the internal structure of data is applied for calculating the phase of a point from its neighborhood. In addition, some algorithms combine the spatial operation and space-frequency analysis. In [34], the local frequencies are estimated using a space-frequency analysis method, and then the phases are calculated using a SCPS formula. From these analyses, we know that the challenging issue for analyzing a spatial carrier fringe pattern with large phase deformation is the estimation of local phase gradients (i.e. local spatial frequencies).

In this paper we present, to the best of our knowledge, a novel technique for analyzing a spatial carrier fringe pattern. First, this technique estimates the local phase gradients at a pixel from its neighborhood, simply by use of statistics of the intensity gradients. Second, it recovers the phase map of the fringe pattern from the phase gradients just estimated, by using numerical integration or by using a 2D adaptive SCPS algorithm. Both numerical simulation and experiment are carried out to verify the validity of this proposed technique.

2. Estimating phase gradients from intensity gradients

2.1. Phase gradients

Consider a spatial carrier fringe pattern of the form

As introduced in Section 1, calculating phase gradients is central to developing an adaptive technique for analyzing a spatial carrier fringe pattern with large phase variations. If the phase gradients are estimated, we can recover the phase map by implementing a numerical integration or by employing a SCPS algorithm with adaptively varying phase shifts. In the following subsections, we will deduce a method for estimating the phase gradients from the intensity gradients.

2.2. Principle of phase gradient estimation

Intensity gradients describe how steeply a fringe pattern varies. They depend mainly on the phase gradients, as well as on the background intensities and modulations. For a pixel (x, y), we can calculate its intensity gradient along the horizontal direction using the forward difference formula

We gain insight into the principle of estimating phase gradients by recalling the temporal phase shift estimation algorithm presented in [35]. With it, the intensity differences between three temporal phase-shifting interferograms are computed first, and then the standard deviations (SDs) of these differences are calculated. The phase shifts are estimated, by using the law of cosines, from a triangle whose lengths of sides are the SDs of the intensity differences. Emulating this temporal method in the spatial domain, the intensity differences between consecutive three pixels are ∇_h^{( + )}I(x, y), ∇_h^(-)I(x, y), and 2∇_HI(x, y). Because the background intensities, modulations, and phase gradients vary slowly across the fringe pattern, we can estimate the SDs of ∇_h^{( + )}I(x, y), ∇_h^(-)I(x, y), and 2∇_HI(x, y) within a sufficiently small neighborhood centered at the pixel (x, y). Assuming the neighborhood size to be (2K + 1) × (2K + 1), their SDs are calculated with

Based on the same principle, the phase gradient along the vertical direction can also be calculated. Using the forward and central difference formulas, we have

The phase gradients calculated using Eqs. (11) and (16) have a range from 0 to π radians, so the sign ambiguity arises. We have to calculate the phase gradients along the diagonal direction for determining the signs of phase gradients. The intensity gradients along the diagonal direction are

2.3. Two-dimensional implementation

The phase gradients calculated in the previous subsection, i.e. ω_h(x, y), ω_v(x, y), and ω_d(x, y) are not practically useful. First, the practical phase gradients along x and y directions are not always positive, but may have the same or opposite signs depending on the fringe direction. Second, it is possible that the fringes are parallel to the horizontal or vertical direction, in which case the phase gradients, calculated in the previous subsection, along this direction may have very large errors. For solving these problems, we calculate the averaging values of σ_h(x, y), σ_v(x, y), and σ_d(x, y) over the whole fringe pattern, i.e. $\overline{{\sigma }_{\text{h}}}$, $\overline{{\sigma }_{\text{v}}}$, and $\overline{{\sigma }_{\text{d}}}$, and then compare them for roughly determining the fringe direction. Among these three averaging values, the smallest one indicates a direction proximate to that of fringes. We should avoid directly calculating the phase gradients in this direction, and a switching mechanism as follows is used.

If $\overline{{\sigma }_{\text{h}}}\ge \overline{{\sigma }_{\text{d}}}$ and $\overline{{\sigma }_{\text{v}}}\ge \overline{{\sigma }_{\text{d}}}$,

If $\overline{{\sigma }_{\text{h}}}\ge \overline{{\sigma }_{\text{v}}}$ and $\overline{{\sigma }_{\text{d}}}\ge \overline{{\sigma }_{\text{v}}}$,

If $\overline{{\sigma }_{\text{d}}}\ge \overline{{\sigma }_{\text{h}}}$ and $\overline{{\sigma }_{\text{v}}}\ge \overline{{\sigma }_{\text{h}}}$,

3. Phase recovering algorithms

3.1 Phase recovering using numerical integration

The phase map can be recovered by numerically integrating the phase gradients just calculated. This task can be implemented by use of least squares techniques [36, 37]. In this subsection, we will present a method based on the 2D discrete cosine transform (DCT) for doing it. Assuming the calculated phase gradients to be M × N matrices, we denote the recovered phase at pixel (x, y) as φ(x, y), the phase gradient at a position between two consecutive pixels can be estimated using the central difference formula u(x + 1/2, y)≈φ(x + 1, y)-φ(x, y). This gradient can also be estimated by averaging u(x, y) and u(x + 1, y), so we have

By padding zeros outside the ranges of the calculated phase gradients and then using Eq. (29), we calculate a matrix of ρ(x, y) with its size becoming (M + 2) × (N + 2), in which case we can regard Eq. (28) as the discretization of a Poisson equation

In the above steps, we can subtract from the phase gradients their mean values, and then calculate the matrix of ρ(x, y) using the zero-mean normalized phase gradients. Doing so will result in a phase map with its carrier and tilt component having been removed.

In fact, the procedure of this phase recovering method is very much similar to that of phase unwrapping technique proposed by Ghiglia and Romero [38], with a difference that the matrix of ρ(x, y) on the right-hand side of Eq. (28) is calculated from the estimated phase gradients. This phase recovering method, which is based on 2D DCT, is suitable for dealing with a rectangular aperture. If the aperture is circular or has other different shapes, we can refer to [36, 37] or some literature regarding weighted phase unwrapping techniques, in order to find a solution for the system of equations in Eq. (28) with (x, y) having non-rectangular ranges.

3.2 Phase recovering using 2D least squares SCPS algorithm

As proposed in [31], estimating the phase gradients allows us recovering the phase map of a spatial carrier fringe pattern using a 2D SCPS algorithm. In the neighborhood of a pixel, the estimated phase gradients are taken as the relative phase steps, so that its phase can be calculated in the least squares sense.

We specify the size of neighborhood to be (2L + 1) × (2L + 1) pixels. By noting that the background intensities and modulations are nearly constant within this neighborhood, and defining c₀(x, y) = a(x, y), c₁(x, y) = b(x, y)cos[u_C x + v_C y + ϕ(x,y)], and c₂(x, y) = -b(x, y)sin[u_C x + v_C y + ϕ(x,y)], we have a system of linear equations

4. Numerical simulations

4.1. Phase gradient estimation

Here we perform numerical simulations for investigating the performances of our proposed approach. Figure 1(a) simulates a spatial carrier fringe pattern of 512 × 512 pixels, with its phase map shown in Fig. 1(b) and its carrier frequency being 0.5π radians/pixels along the direction of 40° from horizontal axis. We assume the background and modulation to be distributed as Gaussian functions with a 50% and 80% decrease in magnitude at the corners, respectively. Zero-mean Gaussian noise with the standard deviation (SD) of 0.02 is added into the intensities.

 

Fig. 1 (a) is a simulated spatial carrier fringe pattern and (b) is its phase map.

Download Full Size | PDF

We use the proposed method in Section 2 for estimating the phase gradients from the spatial carrier fringe pattern in Fig. 1(a). In this procedure, a neighborhood of 11 × 11 pixels is specified for calculating the SDs of intensity gradients. In Fig. 2, the two panels of the left column show the theoretical phase gradients, calculated from the simulated phase map, along x and y directions, and the middle column gives the estimated ones using the proposed method. Comparing them, we find that the estimated phase gradients have almost the same distributions with the theoretical values. The right column illustrates their differences, in which the maximum errors are 0.0320 radians in u and 0.0192 radians in v.

 

Fig. 2 (a) and (d) show the theoretical phase gradients in x and y directions, respectively. (b) and (e) are the estimated phase gradients along the two directions, respectively. (c) and (f) show the differences between the estimated and the theoretical phase gradients in x and y directions, respectively.

Download Full Size | PDF

Further simulations under different noise conditions are performed for investigating the accuracy of this method. As a result, the root-mean square (RMS) errors and the maximum errors in the estimated phase gradients are listed in Table 1, from which it is evident that the errors in the estimated phase gradients increase as the noise SDs increase. Even so, this algorithm still appears stable and robust. When relatively high noise with SD of 0.04 is added, the maximum error keeps being of the level of 0.01 radians, demonstrating that this method is less sensitive to random noise.

Table 1. Simulation Results of Estimating Phase gradients in the Presence of Gaussian Noise

View Table | View all tables in this article

We also see from Table 1 that, in the absence of noise, the estimated phase gradients still contains some errors. They are caused by inaccuracies in estimating the variances of intensity gradients of pixels from their neighborhoods. Using a large neighborhood is helpful for improving the resolution of phase gradients, at the expense of computational time and spatial resolution. In these simulations, we use an 11 × 11 neighborhood. It covers at least two fringe widths, enabling giving satisfying results.

4.2. Phase map recovery

Continuing the simulations, we recover the phase map of the fringe pattern in Fig. 1(a), using the methods presented in Section 3, and compare their results with those of traditional methods. In Fig. 3, the panels of the left column, from top to bottom, show the phase maps recovered using the numerical integration algorithm presented in Section 3.1, the classical 1D SCPS algorithm with the five-step formula, the 2D SCPS algorithm presented in Section 3.2, and the FT method, with their carrier phase components having been removed. By subtracting from these results the predefined phase map in Fig. 1(b), the errors of these techniques are calculated and shown in the right column of Fig. 3. Table 2 summarizes the RMS phase errors of these algorithms.

 

Fig. 3 Phase recovering results of the simulated fringe pattern. The left column, from top to bottom, shows the phase maps extracted using the numerical integration method, 1D SCPS method, 2D SCPS method, and FT method, respectively. The right column illustrates the phase errors corresponding to the above methods in turn.

Download Full Size | PDF

Table 2. Simulation Results of Phase Recovery in the Presence of Gaussian Noise

View Table | View all tables in this article

From Fig. 3 and Table 2, we know that the used four algorithms all recover the phase map successfully, but they achieve different accuracies. With the numerical integration method, additional phase unwrapping step is not required, and the procedure of removing carrier and tilt phase component can be performed in advance by making the phase gradients zero mean normalized. However, because numerical integration involves an inverse high-pass filtering process, by which the signal component within low-band will be magnified significantly, its result may contain large errors of low frequencies, like those shown in Fig. 3(e). The classical 1D SCPS algorithm is fast in computation, but it uses a constant relative phase shift over the whole fringe pattern, thus it has detuning errors in its results, as shown with zigzag-type errors like those shown in Fig. 3(f). The 2D SCPS algorithm adaptively uses the estimated phase gradients as the local relative phase shift, and calculates the phase of each pixel in the least squares sense from a 11 × 11 neighborhood. It achieves very high accuracies, but has a somewhat higher computational complexity in contrast with the 1D algorithm. The FT method, as well known, has nonuniform accuracies over the fringe pattern, and usually the larger errors occur at aperture boundaries where the intensity signal is not continuous. By cutting off the data of aperture boundaries, the FT method has a small error as listed in Table 2.

5. Experiment and discussions

An experiment is carried out for verifying the feasibility of the proposed method in practical applications. Figure 4 shows an interferogram of a plane mirror, with its size being 391 × 391 pixels. In it, each fringe nearly contains five pixels in both the horizontal and vertical directions.

 

Fig. 4 A practical interferogram with spatial carrier.

Download Full Size | PDF

By specifying a neighborhood of 11 × 11 pixels for calculating the variances of intensity gradients, the phase gradients along the horizontal and vertical directions are estimated using the method proposed in Section 2. Their distributions are shown in Fig. 5, from which we observe that the phase gradients vary across the fringe pattern, and have the largest absolute value near the corner with its coordinates being (391,1), corresponding to the top right corner of Fig. 4.

 

Fig. 5 The phase gradients in (a) the horizontal and (b) vertical directions estimated from the interferogram in Fig. 4.

Download Full Size | PDF

Furthermore, we recover the phase map from Fig. 4 using the same methods as those in simulations. Their results are illustrated in Fig. 6, with (a) through (d) in turn being of the numerical integration, the 1D SCPS algorithm, the 2D SCPS algorithm, and the FT method. The phase maps in Fig. 6 are very much alike in their profiles, implying that all these algorithms can recover the phase map successfully. The numerical integration method presented in Section 3.1 avoids phase unwrapping implementation. Its limitation lies in containing errors of low frequency, which is observable when comparing the phase map in Fig. 6(a) with others. The 2D SCPS algorithm presented in Section 3.2 gives a more accurate phase map shown in Fig. 6(c), in contrast with the one shown in Fig. 6(b) which is obtained using the 1D algorithm and contains noticeable zigzag-type errors. The result of 2D SCPS algorithm is also comparable with that of FT method displayed in Fig. 6(d).

 

Fig. 6 The phase maps of the interferogram in Fig. 4, recovered using (a) the numerical integration method, (b) 1D SCPS method, (c) 2D SCPS method, and (d) FT method.

Download Full Size | PDF

The FT method, as mentioned in Section 1, is the mainstay of techniques for spatial carrier fringe pattern analysis, and therefore it can be used as a benchmark for evaluating performances of other methods. In contrast with the FT method which involves a great number of complex multiplications despite the invention of FFT, spatial methods involve much fewer real computations thus being more efficient than the FT method, but generally they achieve lower accuracies. This paper proposes a method for estimating phase gradients from intensity gradients. Its significance for the spatial analysis is that, it can restrain detuning errors of the SCPS technique by taking advantage of its local adaptiveness. The 2D SCPS technique using local phase gradients, as demonstrated by both the simulation and experiment results, can achieve the same and even higher accuracy than the FT method. If using the numerical integration method, the same high accuracy cannot be achieved, but additional procedures for phase unwrapping and for carrier and tilt removing are cancelled.

Besides providing alternatives for the methods in the spatial domain, the proposed algorithm of phase gradient estimation is also potentially useful for enhancing the methods in the frequency domain. It can be combined with the FT method in order to improve its accuracy. An issue with the FT method is the design of filter. By using the filter, one of the lobes in the Fourier spectrum corresponding to fringes is isolated, and the position, size, and shape of the filter affect the measurement accuracies. Using the proposed algorithm allows us determine the average and range of the spatial frequencies of a fringe pattern, so it is helpful for exactly determining the carrier frequency and optimizing the filter parameters. Another issue is regarding boundary errors. The FT method has nonuniform accuracies over the fringe pattern, and usually the large errors appear at aperture boundaries where the intensity signal is not continuous, especially when the aperture is circular or of other shapes. In the situation that the aperture is not rectangular, we usually pad zeros outside the aperture to get a rectangular matrix for implementing 2D FT. Doing so may induce considerably large errors. For compressing the errors, multiplying the aperture with a 2-D Hanning window can partially soften the discontinuity at the aperture boundary [4], but a more effective solution is to perform extrapolation [8, 10, 39] instead of to pad zeros outside the aperture. Using the estimated phase gradients at boundaries enables predicting the intensities outside the aperture, thus extrapolating the aperture. In [31], aperture extrapolation using phase gradients has been demonstrated to be a very effective method for enhancing the accuracy of the FT method.

6. Conclusion

In this paper, we have proposed a novel spatial carrier fringe pattern analysis method. This method allows estimating the local phase gradients at a pixel from its neighborhood, by use of variances of the intensity gradients. Compared with the existing techniques for estimating phase gradients, e.g. those using least-square fitting, iterative strategy, or space-frequency analysis, the proposed algorithm has a simpler principle and is easier to implement. It is robust under noise conditions. In order to recover the phase map using the phase gradients just calculated, we presented a numerical integration method which recovers the phase map by solving a Poisson equation using 2D DCT. This phase recovering method does not require an additional phase unwrapping step, and can remove carrier and tilt phase component in advance by simply making the estimated phase gradients zero mean normalized. Its limitation lies in having low frequency errors in its result. We also introduced a 2D SCPS algorithm calculating the phase of each pixel in the least squares sense from its neighborhood. Because this 2D SCPS algorithm adaptively uses the estimated phase gradients as its local relative phase shift, high measurement accuracies are achieved. In addition, the proposed algorithm for estimating phase gradients also has potential uses in enhancing the methods in the frequency domain.