Fringe Pattern Analysis in Wavelet Domain

Abstract

We present a full-field technique for single fringe pattern analysis based on wavelet transform. Wavelets technique is a powerful method that quantifies at different scales how spatial energy is distributed. In the wavelets domain, fringe pattern analysis requires spatial modulation by a high-frequency carrier. We realize the modulation process numerically by combining the fringe pattern and its quadrature generated analytically by spiral phase transform. The first application concerns the speckle denoising by thresholding the two-dimensional stationary wavelet transform (2D-swt) coefficients of the detail sub-bands. In the second application, the phase derivatives are estimated from the 1D-continuous wavelet transform (1D-cwt) and 2D-cwt analysis of the modulated fringe pattern by extracting the extremum scales from the localized spatial frequencies. In the third application, the phase derivatives distribution is evaluated from the modulated fringe pattern by the maximum ridge of the 2D-cwt coefficients. The final application concerns the evaluation of the optical phase map using two-dimensional discrete wavelet transform (2D-dwt) decomposition of the modulated fringe pattern. The optical phase is computed as the arctangent function of the ratio between the detail components (high-frequency sub-bands) and the approximation components (low-frequency sub-bands). The performance of these methods is tested on numerical simulations and experimental fringes.

1. Introduction

The single fringe pattern is an emerging technique for the analysis of full-field measurements in optical metrology. Fringe pattern analysis becomes a key technique in interferometric metrology which is more suitable for dynamic processes. The main purpose of this chapter is to exploit wavelet transform for single fringe pattern analysis to extract useful information. Wavelets are a powerful method allowing to know the spatial energy distribution at several scales.

Fringe pattern analysis using wavelets transform requires a spatial modulation by a high-frequency carrier as the case of Fourier domain [1]. We realize the modulation process digitally [2] by combining the fringe pattern and its quadrature generated analytically by spiral phase transform SPT [3, 4]. The advantage of this method lies in its implementation of numerical algorithms to reduce experimentation burden present in the phase shifting methods [5], namely, multiple fringe pattern generation and the experimental carrier introduction.

In the first application, we describe a residual speckle noise reduction technique by 2D-stationary wavelet decomposition of the speckle correlation fringe pattern [6]. The speckle noise is removed by thresholding the 2D-swt detail sub-bands coefficients [7].

The second application presents a method for phase derivative estimation based on 1D-continuous wavelet analysis of the modulated fringe pattern using Paul’s function as the mother wavelet [8]. The phase derivative is calculated by extracting the extremum scales from the localized spatial frequencies. The third application is phase derivatives evaluation by 2D-continuous wavelet analysis of the modulated fringe pattern using complex Morlet’s function as the mother wavelet [9]. The phase derivative is estimated by a maximum ridge of wavelets coefficients. Finally, we introduce an algorithm for optical phase extraction by 2D-Discrete wavelet decomposition of the modulated fringe pattern using Gabor’s function as the mother wavelet [10]. The optical phase is extracted as the ratio between the detail components (high-frequency sub-bands) and the approximation components (low-frequency sub-bands).

The remaining part of the chapter proceeds as follow: divided into three sections: the second section will examine a brief description of continuous, discrete and stationary wavelet transforms. In the third section, we briefly introduce the speckle correlation fringes obtained by digital speckle pattern interferometry (DSPI). Finally, the fourth section is devoted to applications of wavelets to fringe pattern analysis.

2. Wavelet transform analysis

The concept of the wavelets and its numerous fields of the application make them a useful tool in several studies concerning localized variations analysis for non-stationary or transient signal analysis. This concept is based on the multiresolution analysis that represents signal variations at different scales. A detailed review of wavelets theory has been published by Ingrid Daubechies in [11]. However, we present here a brief description of three wavelets families for completeness.

2.1. Continuous wavelet transform (cwt)

The analysis of a given signal by continuous wavelets transform concerns to decompose it into several basic functions named wavelets. They are oscillatory functions with a finite duration and having zero average value, also, they are characterized by irregularity and the good localization These properties of the continuous wavelets make them a superior basis for signals analysis with discontinuities. The wavelets are constructed by translating and dilating them other wavelet functions resulting in a self-similar wavelet’s families as follows.

$${\psi }_{s,\xi }\left(x\right)=\psi \left(s^{-1}\left(x-\xi \right)\right)$$ E1

where s and ξ are respectively the scale and translation parameters. The wavelet coefficients obtained by 1D-cwt decomposition of the signal f(x) are given by:

$$w\left(s,\xi \right)=\underset{-\infty }{\overset{+\infty }{\int }}f\left(x\right).{\psi }_{s,\xi }^{\ast }\left(x\right).\mathit{dx}=s^{-1/2}\underset{-\infty }{\overset{+\infty }{\int }}f\left(x\right).{\psi }^{\ast }\left(s^{-1}\left(x-\xi \right)\right).\mathit{dx}$$ E2

where ψ*_s,ξ(x) represents the conjugate of the wavelet ψ_s,ξ(x) defined for each shift ξ and each scale s. Wavelet coefficients are the output of the correlation product between a signal and the mother wavelet for different values of dilatation, and this is interpreted as a measure of the local similarity between them.

The 2D-continuous wavelet coefficients are obtained by the computation of the correlation product between the input image (signal 2D) and the mother wavelet. In the 2D case, a parameter of orientation angle is added, and the 2D-cwt wavelet coefficients of an input image f(x,y) are defined as:

$$w\left(\xi ,\eta ,s,\theta \right)=s^{-2}\iint f\left(x,y\right)\cdot {\psi }^{\ast }\left(s^{-1}{R}_{\theta }\left(x-\xi ,y-\eta \right)\right)\mathit{\text{dxdy}}$$ E3

where ξ and η are respectively the translation parameters along x- and y- directions, s and θ are respectively the scale vector and rotation angle; ψ denotes the 2D mother wavelet and R_θ represents the conventional 2×2 rotation matrix corresponding to θ given by

$$R_{\theta }=\left(x.\text{cos}\theta +y.\text{sin}\theta ,y.\text{cos}\theta -x.\text{sin}\theta \right)$$ E4

2.2. Discrete wavelets transform (dwt)

The discrete wavelet transform denoted by dwt is a fast computing algorithm; it was introduced by Mallat [12] in 1989 for signal or image decomposition into two important components called approximation and details.

The 1D-dwt decomposition is based on two functions called scaling and wavelets. These functions, i.e., the scaling and the mother wavelet functions are related to the impulse response of h[n] and g[n], respectively, as illustrated in Figure 1. The approximation (low-frequency) and details (high-frequency) obtained respectively by scaling and wavelets functions are the down-sampled outputs of the first filters h[n] and g[n].

Figure 1.

Block diagram of filter analysis.

The relation between the two functions and the two filters is expressed as:

$$\varphi \left(x\right)=\sqrt{2}\sum _{n}h\left[n\right]\phi \left(2x-n\right)$$ E5

$$\psi \left(x\right)=\sqrt{2}\sum _{n}g\left[n\right]\psi \left(2x-n\right)$$ E6

g is the conjugate mirror of h:

$$g\left[n\right]={\left(-1\right)}^n.h\left[1-n\right]$$ E7

In Eq. (7), conjugation and mirror effects are represented respectively by ${\left(-1\right)}^n$ and (−n). The internal orthogonality relation is satisfied by the low-pass filter h as follows:

$$\sum _n{h}_n{h}_{n+2j}=0$$ E8

and to have

$$\sum _n{h}_n^2=1$$ E9

The filter g satisfies h, the same internal orthogonality and both obey the mutual orthogonality relation expressed as following:

$$\sum _n{h}_n{g}_{n+2j}=0$$ E10

In 2D-dwt, the wavelet transform of an image involves recursive filtering and sub-sampling process. Figure 2 describes the procedure of dwt analysis, three detail images, and one approximation image are obtained at each level. Concerning detail images, they contain the high-frequency information we denote horizontal image sub-band by HD, vertical image sub-band by VD and the diagonal image sub-band by DD and the approximation image sub-band is denoted by AI which accommodates the low-frequency information.

Figure 2.

Diagrams of dwt image decomposition.

Figure 3.

Block diagram of filter analysis.

The 2D scaling function denoted by ϕ(x,y), and the 3D wavelet ψ^H(x,y), ψ^V(x,y), and ψ^D(x,y) are computed by the algebraic product between the one dimensional scaling and wavelet function as expressed in the following Eq. [13],

$$\begin{array}{l}\varphi \left(x,y\right)=\varphi \left(x\right).\varphi \left(y\right)\\ {\psi }^H\left(x,y\right)=\psi \left(x\right).\varphi \left(y\right)\\ {\psi }^V\left(x,y\right)=\varphi \left(y\right).\psi \left(x\right)\\ {\psi }^D\left(x,y\right)=\psi \left(x\right).\psi \left(y\right)\end{array}$$ E11

The three functions ψ^V, ψ^H, and ψ^D provide respectively the vertical, horizontal and diagonal variations, from there, the three details images are obtained.

2.3. Stationary wavelet transform (swt)

Stationary wavelets transform (swt) was presented in [14]. Swt has the same principle of decomposition as dwt transform, but the process of down-sampling is eliminated which means the swt is translation-invariant. The 2D-swt is founded on the idea of no down-sampling. Specifically, this transform is applied at each pixel of the image and saves the detail coefficients and exploits the low-frequency data at each level. A clear description of this technique is detailed in [15].

The principle of swt analysis is schematized in Figure 3. In brief, the swt technique consists to decompose an input sequence at each level and the given output is a new sequence that has the same length as the original sequences. In order to implement this transform, the process of original image decimation is removed, nonetheless, the size of the filter is modified at each level by zero paddings.

We introduce an operator denoted $\widehat{Z}$ that alternates an input sequence with zeros so that for all the integers j, we have ( ${\widehat{Z}}_x$ )_2j = x_j and ( ${\widehat{Z}}_x$ )_2j + 1 = 0. It is also assumed the filters H^r and G^r to have weights ${\widehat{Z}}^{r}h$ and ${\widehat{Z}}^{r}g$ , respectively. Thus, the filter H^r is characterized by a weight h^r₂^r_j = h_j and h^r_kj = 0 in the case where k is not a multiple of 2^r. The filter H^r is attained by introducing a zero between each adjacent pair of the filter H^(r−1) elements, and an identical trend should be followed for G^r. The following relations satisfy the above statements:

$$D_0^r{H}^{\left[r\right]}={\mathit{HD}}_0^r\text{and}{D}_0^r{G}^{\left[r\right]}={\mathit{GD}}_0^r$$ E12

We start the swt definition by setting u^J to be the original input sequence, the stationary wavelet transform can be described for j = J, J-1, …,1,

$$u^{j-1}=H^{\left[J-j\right]}{u}^j\text{and}{v}^{j-1}=G^{\left[J-j\right]}{u}^j$$ E13

The vectors u^j and v^j acquire the same length for the vector u^J with a length of 2^J.

3. Speckle correlation fringe analysis

3.1. Speckle effect

The speckle appears as a granular structure and it results from the self-interference of a large number of coherent waves randomly scattered from and/or transmitted through a rough object surface [16]. When we illuminate a porous surface with coherent light, the scattered light intensity has a random spatial variation called “speckle effect.” The “speckle pattern” appears chaotic and disordered; it is described by using statistics and probability. The speckle pattern structure is dependent on the coherence properties of the used illumination light and also on the characteristics of the object’s surface diffusion. Locally in space, the electric field of the speckle pattern is given by computing the contribution sum of all illuminated scattering elements of the rough object’s surface and is given by [16]

$$E\left(x,y,z\right)=\sum _{k=1}^n\left|a_k\right|\text{exp}\left(j{\varphi }_k\right)$$ E14

In Eq. (14), n represents the scattering number of elements; a_k and ϕ_k are respectively the amplitude and phase of the k^th scatterer contribution. The theorem of the central limit state that the random variable resulting from the sum of several independent random variables results in a Gaussian distribution when their number (number of an independent random variable) approaches infinity.

In the case where the number of scatterers is large, by applying the central Limit we deduce that:

1. The two components (real and imaginary) of the field are identically distributed Gaussian variables, independent and with zero means.

2. The intensity I(x,y,z) is characterized by a negative exponential probability distribution expressed as [16].

$$p\left(I\right)=\text{exp}\left(-I/<I>\right)/<I>$$ E15

The associated phase ϕ is uniformly distributed between −π and π [16]

$$p\left(\varphi \right)=1/2\pi$$ E16

3.2. Digital speckle pattern interferometry (DSPI)

Digital/electronic speckle pattern interferometry (DSPI/ESPI) is an optical interferometric technique demonstrated simultaneously by Macovski et al. [17], and Butters and Leendertz [18, 19] at the beginning of the 1970s. Later on, the technique was enriched by Beidermann and Ek [20] and Lokberg and Hogmoen [21] with several new investigations. Figure 4 shows the typical schematic of the DSPI system with a sensitivity vector responding to out-of-plane deformation/displacement measurements (Figure 4a), and the digital electronics and image processing unit (Figure 4b). The DSPI technique measures the phase changes introduced by the speckle intensity changes. In this technique, a spatially filtered reference beam is added to the speckle pattern which is scattered from or transmitted through the test object, to code its phase. Therefore, a speckled interferogram, resulting from the superposition of the speckle pattern and the reference beams, contains essential information about the random phase. In DSPI, as shown in Figure 4, two speckle interferograms, corresponding to two different states of the object are recorded and stored in the frame grabber card. For example, one speckle interferogram, A, is recorded when the object is in its initial state, say the reference state, while another speckle interferogram, B, is recorded in the deformed state of the object, say by a small distance “d,” as shown in Figure 4a. The DSPI fringe pattern is observed after subtraction of these two speckle interferograms A and B by using the appropriate software. The useful information about the object is encoded in the DSPI fringe pattern like the one shown in Figure 4b.

Figure 4.

Schematic of the DSPI system: (a) Speckle interferometer and (b) Digital electronics and image processing unit.

The DSPI system can be made sensitive to out-of-plane or in-plane deformations, or both, depending on the design of the optical setup. Numerous measurement applications of DSPI have been investigated with different proposed experimental configurations [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48] which have immense importance in the scientific, engineering, and industrial fields. The popularity of the technique, in optical metrology, is by virtue of the several advantages it offers such as: it is a non-contact/invasive type technique; it provides full-field of view information; it is faster in operation as almost real-time observations can be obtained; investigations on large-size objects can be observed with the use of proper optics, and the technique is less sensitive to environmental perturbations in comparison with its counterparts. The DSPI technique has proven its capability in many investigations, e.g., for the measurement of displacement/deformation with variable sensitivity to in-plane and out-of-plane direction [22, 23, 24, 25, 26, 27], measurement of three-dimensional shape [28, 29, 30], surface roughness measurement [31], vibration measurement/monitoring [32, 33, 34], measurement of material properties [35, 36, 37], flow visualization [38, 39], measurement of refractive index and temperature distribution [40, 41, 42], and the investigation of the magnetic fields on the temperature profile of gaseous flames [43, 44]. A Doctoral Thesis describing the comprehensive study of DSPI and some of its novel investigations is recommended for interested readers working in this field [48].

Mathematically, the fringe formation in DSPI is demonstrated here. The intensity distribution recorded before the object deformation is expressed as follows:

$$f_B\left(x,y\right)=b\left(x,y\right)+a\left(x,y\right).\text{cos}{\varphi }_s\left(x,y\right)$$ E17

In this equation, the term b(x,y) represents the background or bias, a(x,y) the visibility and term ϕ_s(x,y) is the original speckle phase appearing as a high frequency, and the intensity of pixel by pixel is randomly distributed. When the object is deformed, the intensity distribution of speckle-pattern becomes:

$$f_A\left(x,y\right)=b\left(x,y\right)+a\left(x,y\right).\text{cos}\left({\varphi }_s+\phi \right).$$ E18

where φ(x,y) represents the phase variation stemming from the object deformation. Assuming that the introduced deformations are sufficiently small, we can ignore the speckle decorrelation effects.

The speckle fringe correlation intensity distribution in the subtraction mode is expressed as:

$$f=f_{\mathrm{A}}-f_{\mathrm{B}}=2.a.\text{sin}\left(\mathrm{\phi }/2\right).\text{sin}\left({\varphi }_{\mathrm{s}}+\mathrm{\phi }/2\right)$$ E19

The above equation can be further expressed as:

$$f^2=4a^2.{\text{sin}}^2\left(\mathrm{\phi }/2\right).{\text{sin}}^2\left({\varphi }_{\mathrm{s}}+\mathrm{\phi }/2\right)$$ E20

The speckle phase ϕ_s(x,y) changes faster across the speckle pattern, the second sine squared term has an ensemble average expressed as:

$${\text{sin}}^2\left({\varphi }_{\mathrm{s}}+\mathrm{\phi }/2\right)\approx \left(1/2\mathrm{\pi n}\right)\underset{0}{\overset{2\pi }{\int }}{\text{sin}}^2\left({\varphi }_{\mathrm{s}}+\mathrm{\phi }/2\right)\mathrm{d}\left({\varphi }_{\mathrm{s}}+\mathrm{\phi }/2\right)=1/2$$ E21

Hence, the speckle fringe correlation intensity distribution becomes:

$$f_{\mathrm{C}}=f^2=2{\mathrm{a}}^2.{\text{sin}}^2\left(\mathrm{\phi }/2\right)={\mathrm{a}}^2.\left(1-\text{cos}\mathrm{\phi }\right)$$ E22

The obtained intensity distribution f_c is rendered as the classical form of a cosine fringe pattern.

3.3. Fringe pattern analysis

The analysis of the fringe patterns concerns to evaluate the coded phase distribution related to the physical magnitude. We classified the phase extraction techniques to methods that explore a shifted fringe pattern as the phase shifting techniques and methods that require the introduction of a spatial carrier such as the Fourier transform method and the wavelet method.

The phase shifting (see Ref. [5]) and the Fourier transform methods (see Ref. [1]) are the most common methods exploited for fringe pattern analysis. The goal of the next section is to present algorithms based on wavelets transform to analyze a single frame speckle fringe pattern.

4. Applications to speckle correlation fringe analysis

As we stated previously, the fringe pattern analysis in wavelet domain requires a modulation process by introducing a spatial carrier. In our case, the spatial carrier in generated numerically and the entire process is carried out directly using a computer. Before presenting the four applications, we present in the next subsection the modulation process required by wavelets transform.

4.1. Fringe pattern modulation

In order to modulate a given fringes pattern in a chosen direction, we combine the fringe pattern and its quadrature with a cos(m.x) and sin(m.x) matrix respectively where m means the modulation rate. We consider the fringe pattern intensity distribution expressed as follow:

$$f\left(x,y\right)=b\left(x,y\right)+a\left(x,y\right)\cdot \text{cos}\left(\phi \left(x,y\right)\right)$$ E23

A low-pass filter applied to the intensity distribution removes the background illumination and Eq. (23) becomes:

$$\tilde{f}\left(x,y\right)=a\left(x,y\right)\cdot \text{cos}\left(\phi \left(x,y\right)\right)$$ E24

Larkin proposed recently a transform called spiral phase quadrature denoted -SPT- for the 2D fringe pattern (see Refs. [3, 4]). The SPT of $\stackrel{\sim }{f}$ is computed as:

$$\mathit{SPT}\left(\tilde{f}\right)=j\cdot \text{exp}\left(j\cdot D\left(x,y\right)\right)\cdot a\left(x,y\right)\cdot \text{sin}\left(\phi \right)$$ E25

where, a(x,y).sinφ is the quadrature term, j² = −1, and D(x,y) represents the direction map. Then, we obtain the sine fringe pattern as:

$$a\left(x,y\right)\cdot \text{sin}\left(\phi \right)=-j\cdot \text{exp}\left(-j\cdot D\left(x,y\right)\right)\cdot \mathit{SPT}\left(\tilde{f}\right)$$ E26

The direction map is computed by the ratio between the horizontal and vertical gradient of the $\stackrel{\sim }{f}$ as and analytically expressed as:

$$\text{tan}\left(D\right)={\nabla }_y\tilde{f}/{\nabla }_x\tilde{f}$$ E27

So, we define the quadrature map as:

$$q\left(x,y\right)=b\left(x,y\right)\cdot \text{sin}\left(\phi \right)=-j\cdot \text{exp}\left(-j\cdot D\right)\cdot \mathit{SPT}\left(\tilde{f}\right)$$ E28

Then, the intensity distribution of the modulated fringe pattern is defined as:

$$f_m=\tilde{f}\cdot \text{cos}\left(m\cdot x\right)-q\cdot \text{sin}\left(m\cdot x\right)=b\left(x,y\right)\cdot \text{cos}\left(\phi +m\cdot x\right)$$ E29

4.2. Speckle noise reduction using various wavelet-based techniques: a simulation study

4.2.1. Speckle noise reduction using 2D-swt

As a reminder, the speckle fringe pattern intensity distribution in the subtraction mode is expressed in (Eq. (19)) as $f=2.a.\text{sin}\left(\mathrm{\phi }/2\right).\text{sin}\left({\varphi }_{\mathrm{s}}+\mathrm{\phi }/2\right)$ . The speckle fringes, presented in Eq. (19), are characterized by multiplicative residual speckle noise, therefore, a denoising step is necessary before the evaluation of the phase distribution and their derivatives. The high-frequency sin(ϕ_s + φ/2) noise should be removed with an appropriate filtering technique. After decomposition by swt of the speckle fringe pattern, the noise is reduced by thresholding the detail sub-bands coefficients using a soft threshold function. The choice of this type of thresholding stems from its characteristic to obtain near-optimal minimax rate. Generally, the threshold estimation requires the knowledge of the noise variance and the optimal threshold is estimated by taking into consideration the variance (σ²) of the coefficients in the highest wavelets decomposition. Figure 5 illustrates the capability of stationary wavelet transform thresholding technique to reduce the speckle noise. The left and the right halves present the intensity distribution of speckle fringe correlation before and after denoising. It is clearly shown in the line profile plot in Figure 5b along the line AB, the influence of swt to remove the frequency from the noised image.

Figure 5.

(a) Speckle fringe correlation before and after denoising; and (b) the 1D intensity profile along line AB.

4.2.2. Phase derivatives extraction using 1D-cwt

The idea is the phase extraction by spatial frequencies using the1D-continuous wavelet transform and Paul’s function as the mother wavelet. The phase derivatives are evaluated from the modulus of the cwt coefficients by extracting the extremum scales from the localized spatial frequencies and the phase is obtained by numerical integration.

Generally, the modulated fringe pattern intensity distribution, as represented by Eq. (18), can also be expressed as:

$$f_m\left(x,y\right)=b\left(x,y\right)+a\left(x,y\right)\text{cos}\left(m.x+\phi \left(x,y\right)\right)$$ E30

Using Eq. (2), the wavelet transform coefficients of the modulated fringe pattern, given in Eq. (30), can be represented as

$$w\left(x,s,\xi \right)=s^{-1/2}{\int }_{-\infty }^{+\infty }\left[b\left(x,y\right)+a(x,y)\text{cos}\left(m.x+\phi \left(x,y\right)\right)\right]{\left(\psi \left(\left(y-\xi \right)/s\right)\right)}^{\ast }\mathit{dy}$$ E31

The localization property of the wavelet is taken into account in order to write the Taylor series of the phase φ(x,y) near the central value ξ as:

$$\phi \left(x,y\right)=\phi \left(x,\xi \right)+\left(y-\xi \right)\partial \phi \left(x,y\right)/\partial y+\dots$$ E32

Assuming that the bias a and the visibility b have a slow variation, so the higher order of (y-ξ) is neglected with respect to the phase-modulated carrier because of the localization of the wavelet. With these considerations, the 1D-cwt coefficients are now computed as:

$$w\left(x,s,\xi \right)=s^{-1/2}b\left(x,\xi \right)a\left(x,\xi \right){\int }_{-\infty }^{+\infty }\text{cos}\left[\mathit{mx}+\phi \left(x,\xi \right)+\left(y-\xi \right)\partial \phi (x,\xi )/\partial y\right]{\left(\psi \left(\left(y-\xi \right)/s\right)\right)}^{\ast }\mathit{dy}$$ E33

and the Parseval’s identity gives

$$w\left(x,s,\xi \right)=\frac{\sqrt{s}}{2\pi }{\int }_{-\infty }^{+\infty }\widehat{f_a}\left(x,k\right){\left(\widehat{\psi }\left(\mathit{sk}\right)\right)}^{\ast }{e}^{\mathit{i\xi k}}\mathit{dk}$$ E34

where

$$\widehat{f_a}\left(x,k\right)=b\left(x,\xi \right)a\left(x,\xi \right)\mathit{\pi h}\left(x,k\right)$$ E35

with

$$h\left(x,k\right)=\delta \left(k-m_1\right)\text{exp}\left[i\left(\phi \left(x,\xi \right)-\xi \partial \phi (x,\xi )/\partial y\right)\right]+\delta \left(k+m_1\right)\text{exp}\left[-i\left(\phi \left(x,\xi \right)-\xi \partial \phi (x,\xi )/\partial y\right)\right]$$ E36

$$m_1=m+\partial \phi \left(x,\xi \right)/\partial x$$ E37

The wavelet transform reduces to.

$$w\left(x,s,\xi \right)=0.5.a\left(x,\xi \right)b\left(x,\xi \right)\sqrt{s}\left[{\left(\widehat{\psi }\left({\mathit{sm}}_1\right)\right)}^{\ast }\text{exp}\left(i\left(\mathit{\xi m}+\phi \left(x,\xi \right)\right)\right)+{\left(\widehat{\psi }\left(-{\mathit{sm}}_1\right)\right)}^{\ast }\text{exp}(-i(\mathit{\xi m}+\phi (x,\xi \left)\right))\right]$$ E38

Introducing Paul’s mother wavelet defined as:

$$\psi \left(x\right)=\left(2^{n}n!{\left(1-\mathit{ix}\right)}^{-\left(n+1\right)}\right)/\left(2\pi \sqrt{\left(2n\right)!/2}\right)$$ E39

we get

$$w\left(x,s,\xi \right)={\left(\left(2n\right)!\right)}^{-1}a\left(x,\xi \right)b\left(x,\xi \right)s^{n+1/2}{m}_1^n\text{exp}\left(-{\mathit{sm}}_1\right)\times \text{exp}\left(i\left(\mathit{\xi m}+\phi \left(x,\xi \right)\right)\right)$$ E40

whose modulus is computed as

$$\left|w\left(x,s,\xi \right)\right|={\left(\left(2n\right)!\right)}^{-1}a\left(x,\xi \right)b\left(x,\xi \right)m_1^n{s}^{n+1/2}\text{exp}\left(-{\mathit{sm}}_1\right)$$ E41

The extremum scale denoted by S is represented by

$$S\left(x,\xi \right)=\left(2n+1\right)/2m_1$$ E42

Equations (37) and (42) provide the phase derivative as

$$\partial \phi \left(x,\xi \right)/\partial y=\left(2n+1\right)/2S\left(x,\xi \right)-m$$ E43

where m is the modulation ratio.

The performance of this algorithm is illustrated in Figure 6, where the phase derivatives of the speckle fringe correlation along x- and y-directions are presented.

Figure 6.

The evaluated phase derivatives along (a) x-direction and (b) y-direction using 1D-cwt technique.

4.2.3. Phase derivatives extraction using 2D-cwt

The idea is to evaluate the phase derivative by ridge point using2D-continuous wavelet transform and Morlet’s function as the mother wavelet. The phase derivative is computed from maximums scales relating to the ridge point of the wavelet coefficient modules. From the modulated fringe pattern f_m(x, y), we compute the 2D-cwt wavelet coefficients as defined in Eq. (3) and represented again as

$$w\left(\mathrm{t},\mathrm{d},s,\theta \right)=s^{-2}\iint f_m\left(x,y\right)\cdot {{\psi }_M}^{\ast }\left(s^{-1}{R}_{\theta }\left(x-t,y-d\right)\right)\mathit{\text{dxdy}}$$ E44

The 2D complex Morlet is used as the mother wavelet: it is essentially a plane wave modulated by a Gaussian window, and is expressed as:

$${\psi }_M=\text{exp}\left(-\left(x^2+y^2\right)/2\right)\cdot \text{exp}({\mathit{jk}}_0\left(x\cdot \text{cos}\theta +y\cdot \text{sin}\theta \right)$$ E45

where k₀ is the specific spatial frequency that should be in the range 5–6 in order to satisfy the admissibility condition [49], and j² = −1.

A new matrix containing the maximum value of each column of the wavelet coefficient modulus array defines the wavelet ridge, and then, the corresponding scale value is determined from the ridge wavelet [50, 51]. The maximum scales s_max correspond to the maximum ridge of the obtained wavelet coefficients modulus:

$$\left(s_{\max },\theta \right)=\underset{s\in R^{+},\theta \in \left[0,2\pi \right]}{\text{arg}\text{max}}\left|w\left(t,d,s,\theta \right)\right|$$ E46

where s_max represents the scale value for maxima. The phase derivative is obtained by

$$\nabla \phi =\left(k_0+{\left(k_0^2+2\right)}^{1/2}/2s_{\max }\right)-m$$ E47

The horizontal and vertical phase derivatives evaluated from the speckle fringe correlation are presented in Figure 7.

Figure 7.

The evaluated phase derivatives along (a) x-direction and (b) y-direction using 2D-cwt technique.

4.2.4. Optical phase extraction using 2D-dwt

This application concerns to extract the optical phase distribution by dwt components using the 2D discrete wavelets transform and Gabor’s function as the mother wavelet. The 2D-dwt method decomposes the fringe pattern into two components: approximation (low-frequency) and details (high-frequency) using two principal quadrature mirror decomposition filters: low-pass and high-pass filters.

We consider the modulated fringes expressed in Eq. (30) as $f_m\left(x,y\right)=b\left(x,y\right)+a\left(x,y\right).\text{cos}\left(m.x+\phi \left(x,y\right)\right)$ . It depends on the desired phase φ(x,y). Figure 8 presents the four sub-bands of the modulated speckle fringe correlation after 2D-dwt decomposition.

Figure 8.

The four sub-bands images after modulated fringe 2D-dwt decomposition.

The optical phase distribution coded in the modulated fringe pattern is evaluated by computing the arctangent function of the ratio between the detail and approximation images as following:

$$\phi \left(x,y\right)=\text{atan}2\left(\left(f_{\mathrm{m}}\ast \mathrm{h}\right)/\left(f_{\mathrm{m}}\ast \mathrm{g}\right)\right)$$ E48

Inside the atan2 function, the numerator and denominator present the convolution product between f_m and the two filters h and g, and give the detail and approximation images respectively. An illustration of the two images is shown in Figure 9.

Figure 9.

Approximation and detail image of speckle fringe correlation.

The detailed image in the numerator is chosen according to the direction of modulation. The low-pass filter h and high-pass filters g are respectively the real and imaginary part of Gabor function expressed as:

$$h=\text{exp}\left[-\left(x^2/{\sigma }_x^2+y^2/{\sigma }_y^2\right)/2\right].\text{cos}\left(2\pi f_0\left(x\text{cos}\theta +y\text{sin}\theta \right)\right)$$ E49

$$g=\text{exp}\left[-\left(x^2/{\sigma }_x^2+y^2/{\sigma }_y^2\right)/2\right].\text{sin}\left(2\pi f_0\left(x\text{cos}\theta +y\text{sin}\theta \right)\right)$$ E50

The 2D Gabor’s function [52] is the extension of the 1D Gabor function introduced by Gabor [53]. It is defined by a Gaussian window that modulates a complex exponential centered at a fixed frequency It is formulated by:

$$\mathrm{G}=\text{exp}\left[-\left({\mathrm{x}}^2/{\mathrm{\sigma }}_{\mathrm{x}}^2+{\mathrm{y}}^2/{\mathrm{\sigma }}_{\mathrm{y}}^2\right)/2\right].\text{exp}\left(j.2\pi f_0\left(x\text{cos}\theta +y\text{sin}\theta \right)\right)$$ E51

Central frequency, orientation, spatial extent, and aspect are four essential parameters characterize Gabor function, and they can be adjusted by varying respectively f₀, θ, σ_x and σ_y. The use of the atan2 function in Eq. (48) provides the phase distribution modulo 2π as shown in Figure 10a. To remove this discontinuity, a phase unwrapping step is obligatory [54]; for this reason, we have implemented the PUMA algorithm [55]. The unwrapped phase illustration in three-dimensional representations is presented in Figure 10b.

Figure 10.

(a) Wrapped phase distribution retrieved by 2-dwt; and (b) continuous phase distribution volume visualization.

4.3. Discussion of wavelet techniques on the numerically simulated fringe correlations

In order to study the performance and validity of the algorithms, an image quality assessment (Q) is used [56]. This metric model any image distortion by combining three factors. The first factor is the loss of correlation between the ideal and the obtained images denoted O and E respectively, which measures the linear correlation degree between the two input images. This factor is defined as:

$$Q_1={\sigma }_{\mathit{OE}}/{\sigma }_O{\sigma }_E$$ E52

The second factor is luminance distortion, it measures the mean luminance between O and E, which its value is computed as:

$$Q_2=2\overline{O}.\overline{E}/\left({\left(\overline{O}\right)}^2+{\left(\overline{E}\right)}^2\right)$$ E53

The third factor is contrast distortion that measures the contrast similarity between the two images and it is defined as:

$$Q_3=2{\sigma }_O{\sigma }_E/\left({\sigma }_O^2+{\sigma }_E^2\right)$$ E54

And then, the Q coefficient is computed by the product of three factors:

$$Q=Q_1\times Q_2\times Q_3$$ E55

where $\overline{O}$ and $\overline{E}$ present average of O and E, respectively, ${\sigma }_O$ and ${\sigma }_E$ are respectively the standard deviations of the two images. The Q index takes values between −1 and 1 where 1 means that the retrieval characteristic is exact. The Table 1 below summarizes the computed Q index values.

Retrieved characteristic map	Q index value	Time in seconds
Horizontal phase gradient	0.95 (1-cwt) & 0.96 (2-cwt)	50.23 (1-cwt) & 20 (2-cwt)
Vertical phase gradient	0.95 (1-cwt) & 0.96 (2-cwt)	50.26 (1-cwt) & 20.11 (2-cwt)
Recovered phase distribution by 2-dwt	0.97	10.55

Table 1.

The measured Q index values.

As a note, the 1D-cwt analyses images line by line, whereas the 2D-cwt scans images with the help of the parameter of orientation. On the analysis time side, we used Matlab on a 2.93 GHz Intel Pentium processor machine with 4 GB RAM. According to the Q index values, the presented algorithms give the desired information with good accuracy (Q between 0.95 and 0.97).

4.4. Application on an experimental speckle fringe correlation

After the presentation of the four applications of the stationary, continuous and discrete wavelet transform for fringe pattern analysis, and the favorable effectiveness validated by computer simulation, we exploit this algorithmic arsenal to analyze an experimental speckle fringe pattern recorded by using the DSPI setup showed in Figure 11a.

Figure 11.

(a) Schematic of DSPI setup [57, 58], (b) Recorded DSPI fringe pattern, and (c) Denoised DSPI fringe pattern by using swt thresholding technique.

In DSPI, two specklegrams corresponding to the two different states (before and after deformation) of the object are recorded with an image sensor and then subtracted in order to get the speckle fringe correlation as shown in Figure 11b. Figure 11c shows the denoised fringe resulting from the swt thresholding technique. The wrapped phase distribution and corresponding three-dimensional continuous phase distribution are shown in Figures 12a and b, respectively. The horizontal and vertical phase derivatives are shown in Figure 13.

Figure 12.

(a) Wrapped phase distribution extracted using dwt algorithm. (b) Three-dimensional representation of continuous phase distribution.

Figure 13.

The estimated phase gradient distribution using 2-cwt technique along: (a) x-direction and (b) y-direction.

5. Conclusion

This chapter has introduced an important thematic in interferometric metrology namely, fringe pattern analysis. We have discussed the capability of wavelets transform families (stationary, continuous and discrete wavelet) to analyze the fringe patterns.

This analysis consists in extracting the optical phase distribution coded in the fringes and its horizontal and vertical derivatives. We have presented with the help of computer simulations four applications of the wavelet transform, in order to analyze the speckle fringes correlation obtained in digital speckle pattern interferometry (DSPI). We summarize these applications as:

• Speckle noise reduction using stationary wavelet transform thresholding technique.

• Phase derivative extraction using one and two-dimensional continuous wavelet transform.

• Optical phase distribution extraction using the two-dimensional discrete wavelet transform

Acknowledgments

Manoj Kumar and Osamu Matoba acknowledge the support by Japan Society for the Promotion of Science KAKENHI (17F17370, 18H03888) and JST CREST Grant Number JPMJCR1755. All the authors wish to thank the anonymous reviewers for their great insight, comments, and suggestions.