1. Introduction

In digital holographic interferometry (DHI), the derivative of the phase of reconstructed interference field contains information about the displacement derivative or strain distribution, a parameter of great interest in the areas of experimental mechanics and non-destructive testing. Direct estimation of the phase derivative is important because the process of phase demodulation and subsequent differentiation is error-prone in the presence of noise [1]. Various techniques have been developed for the direct measurement of phase derivatives. A Fourier transform based method [2] was proposed for direct phase derivative estimation. Another approach which is quite popular for phase derivative measurement is by shearing the complex amplitude of the reconstructed interference wave and superimposing on the original wave [3–5]; though the sensitivity of the technique depends on the amount of shearing introduced. Besides, it provides wrapped and noisy estimates requiring unwrapping and filtering operations [6]. Some other approaches for direct phase derivative estimation include phase shifting windowed Fourier ridges method [7], grid method [8] and space-frequency distributions based methods [9,10]. Recently, methods [11,12] relying on piecewise polynomial approximation of phase have been proposed.

In this paper, we propose an elegant method to directly estimate the phase derivatives in DHI without requiring any unwrapping operation. The method relies on joint space-frequency representation of the complex reconstructed interference field in DHI by evaluating the corresponding two-dimensional pseudo Wigner-Ville distribution (2D-PSWVD). The phase derivative is subsequently estimated by tracing the peak of the distribution. The theory of the proposed method is discussed in the next section. Simulation and experimental results are presented in section 3 followed by conclusions and acknowledgments.

2. Theory

For deformation analysis using DHI, two holograms are recorded for object states before and after deformation and the corresponding complex amplitudes are obtained using numerical reconstruction [13]. Multiplying the post-deformation complex amplitude with the conjugate of pre-deformation complex amplitude, the reconstructed interference field is obtained which is given as [5, 14]:

where A(x,y) is the amplitude term; ϕ(x,y) is the interference phase and η(x,y) represents the noise assumed to be zero mean additive white Gaussian noise (AWGN). Here x and y refer to the pixel values along the N × N size image. The real part of the reconstructed interference field constitutes a fringe pattern. The 2D Wigner-Ville distribution corresponding to I(x,y) is given as [15]

where ‘*’ denotes the complex conjugate. To improve the robustness of the distribution against interference terms or artifacts, a smoothing window is introduced [16] and the resulting so called two-dimensional pseudo Wigner-Ville distribution can be equivalently written as

where w(τ ₁, τ ₂) is a real symmetric window usually having a low-pass behavior. From, Eq. (3), it is clear that 2D-PSWVD maps the reconstructed interference field I(x,y) in spatial domain to Γ(x,y,ω ₁,ω ₂), thereby providing a joint space-frequency representation. For the simplicity of analysis, the amplitude term is assumed constant within the window region and the noise terms are ignored. Hence from Eq. (1) and Eq. (3), we have

Assuming that the phase ϕ(x,y) is slowly varying in the neighborhood of (x,y) and using Taylor series expansion up to second order, we have the following approximations

where

${\varphi }_x(x,y)=\frac{\partial \varphi (x,y)}{\partial x},{\varphi }_{xx}(x,y)=\frac{{\partial }^2\varphi (x,y)}{\partial x^2},{\varphi }_{xy}(x,y)=\frac{{\partial }^2\varphi (x,y)}{\partial x\partial y}$

and

${\varphi }_y(x,y)=\frac{\partial \varphi (x,y)}{\partial y},{\varphi }_{yy}(x,y)=\frac{{\partial }^2\varphi (x,y)}{\partial y^2}$

Using these approximations in Eq. (4), we have

which is essentially the 2D Fourier transform of w(τ ₁,τ ₂) modulated by the term exp[2j(ϕ_x(x,y)τ ₁ + ϕ_y(x,y)τ ₂)]. Denoting Ŵ(ω ₁,ω ₂) as the 2D Fourier transform of w(τ ₁,τ ₂) and using the frequency shifting property of the Fourier transform, we have

As window w has a low-pass behavior, Γ(x,y,ω ₁,ω ₂) is maximum for Ŵ(0,0) i.e. when [ω ₁,ω ₂] = [ϕ_x(x,y),ϕ_y(x,y)]. In other words,

 

Fig. 1. (a) Simulated fringe pattern. (b) Estimated phase derivative along x in radians/pixel. (c) Wrapped phase derivative along x. (d) Error between original and estimated phase derivatives along x in radians/pixel. (e) Estimated phase derivative along y in radians/pixel. (f) Wrapped phase derivative along y. (g) Error between original and estimated phase derivatives along y in radians/pixel.

Download Full Size | PDF

Hence, the phase derivatives at pixel location (x,y) can be estimated by finding the values of (ω ₁,ω ₂) at which Γ(x,y,ω ₁,ω ₂) attains its maximum. Note that the phase derivatives along both directions i.e. x and y can be simultaneously determined using Eq. (10). By repeating the above analysis for all pixels, the phase derivatives for the entire fringe pattern can be estimated. The computational burden for evaluating the 2D Fourier transform can be greatly relieved by using 2D fast Fourier transform (FFT) algorithm.

3. Simulation and Experimental Results

A complex reconstructed interference field signal was simulated at signal to noise ratio (SNR) of 10 dB using ‘AWGN’ function of MATLAB 7.8.0(R2009a) and the corresponding fringe pattern (size 128 × 128) is shown in Fig. 1(a). To apply the 2D-PSWVD for phase derivative estimation, a 2D Gaussian window of the following form was used:

Table 1. Comparison of RMSEs (in radians/pixel) for phase derivative estimation along x at various SNRs (in dB)

View Table

with σ_x = σ_y = 16. Applying the proposed method, the obtained phase derivative estimate along x in radians/pixel and its wrapped form are shown in Fig. 1(b) and Fig. 1(c). The phase derivative estimate along y in radians/pixel and its wrapped form are shown in Fig. 1(e) and Fig. 1(f). Note that the proposed method directly provides unwrapped estimates and the wrapped forms are shown here for the sake of illustration only. The errors between original and estimated phase derivatives along x and y are shown in Fig. 1(d) and Fig. 1(g) respectively. The root mean square errors (RMSE) for phase derivative estimation along x and y are 0.0048 and 0.0036 radians/pixel respectively. The pixels near the borders of the fringe pattern were neglected in the analysis to ignore the errors at the boundaries.

In DHI measurements, noise is an important consideration and its various sources are discussed in [17]. The major advantage of the proposed method is the improved robustness against noise as compared to the phase derivative estimation method used in [10]. This is because a 2D window captures data samples along both x and y for the evaluation of the PSWVD in the proposed method and consequently provides smoothing in both directions effectively increasing the method’s immunity against noise. The increased noise robustness comes at the cost of increased computational burden due to the need of evaluating 2D Fourier transforms at each pixel location. On the contrary, the method in [10] uses a window along only one direction and though it is quite computationally efficient due to the use of 1D Fourier transform, it is more error-prone under high noise. To compare the performance of both methods with respect to noise-immunity, a reconstructed interference field signal similar to the one used for Fig. 1 was simulated for different SNRs. The RMSEs (in radians/pixel) for phase derivative estimation obtained by the proposed method and the method in [10] at various SNRs (in dB) are shown in Table 1. From the table, it is clear that the proposed method exhibits superior performance and provides phase derivative estimates with reasonable accuracy even for very low SNR such as 2 dB. The computational times required for the proposed method and the other method are about 8 minutes and 3 seconds on a 2.66 GHz Intel Core 2 Quad Processor machine with 3.23 GB random access memory. Though the computational cost of the proposed method is high, it can be significantly improved using techniques like fixed-point computations [18] and graphics processing units [19].

To test the practical applicability of the proposed method, a DHI experiment was conducted by subjecting a circularly clamped object to central load and recording two holograms before and after deformation using a Coherent Verdi laser (532 nm). The complex amplitudes of the object wave before and after deformation were obtained by numerical reconstruction performed using discrete Fresnel transform [13]. The reconstructed interference field was obtained as discussed earlier in the ‘Theory’ section. The experimental fringe pattern is shown in Fig. 2(a). The phase derivative estimate along x in radians/pixel and the corresponding wrapped form obtained using the proposed method are shown in Fig. 2(b) and Fig. 2(c). Similarly, the phase derivative estimate along y in radians/pixel and the corresponding wrapped form are shown in Fig. 2(d) and Fig. 2(e). For comparison, the digital shearing approach [5] was also applied for phase derivative estimation along y and the obtained wrapped estimate is shown in Fig. 2(f). It is clear from Fig. 2(e) and Fig. 2(f) that the proposed method shows superior performance. Moreover, since the obtained estimates are wrapped, unwrapping operation is required in digital shearing approach.

 

Fig. 2. (a) Experimental fringe pattern. (b) Estimated phase derivative along x in radians/pixel. (c) Wrapped form of phase derivative along x. (d) Estimated phase derivative along y in radians/pixel. (e) Wrapped form of phase derivative along y. (f) Wrapped estimate along y using digital shearing approach.

Download Full Size | PDF

4. Conclusions

A 2D-PSWVD based method is proposed in the paper for the direct estimation of the phase derivative in DHI. The method provides phase derivative estimates without requiring any unwrapping or filtering operations. The method also has high robustness against noise and enables simultaneous retrieval of phase derivatives along both spatial directions. The authors believe that the method could prove to be a potential choice for DHI applications involving measurement of displacement derivatives.