The 2D Continuous Wavelet Transform: Applications in Fringe Pattern Processing for Optical Measurement Techniques

Abstract

Optical metrology and interferometry are widely known disciplines that study and develop techniques to measure physical quantities such as dimensions, force, temperature, stress, etc. A key part of these disciplines is the processing of interferograms, also called fringe patterns. Owing that this kind of images contains the information of interest in a codified form, processing them is of main relevance and has been a widely studied topic for many years. Several mathematical tools have been used to analyze fringe patterns, from the classic Fourier analysis to regularization methods. Some methods based on wavelet theory have been proposed for this purpose in the last years and have evidenced virtues to consider them as a good alternative for fringe pattern analysis. In this chapter, we resume the theoretical basis of fringe pattern image formation and processing, and some of the most relevant applications of the 2D continuous wavelet transform (CWT) in fringe pattern analysis.

1. Introduction

Fringe pattern processing has been an interesting topic in optical metrology and interferometry; owing to its relevance nowadays, it is a widely studied discipline. Digital fringe pattern processing is used in optical measurement techniques such as optical testing [1, 2], electronic speckle pattern interferometry (ESPI), holographic interferometry, and moiré interferometry or profilometry [3, 4, 5]. They are quite popular for non-contact measurements in engineering and have been applied for measuring various physical quantities like displacement, strain, surface profile, refractive index, etc. In optical methods of measurement, the phase, which is related to the measured physical quantity, is encoded in an intensity distribution represented in an image which is, in general, the result of the interference phenomena. This phenomenon is used in classical interferometry, in holographic interferometry, and in electronic speckle pattern interferometry to convert the phase of a wave of interest into an intensity distribution. As the physical quantity to be measured is codified as the phase of a fringe pattern image, the main task of fringe pattern processing is to recover such phase.

The methods for phase recovery from fringe patterns can be classified mainly in three categories [2, 6]: (a) Phase-stepping or phase-shifting methods which require a series of fringe images to recover the phase information. (b) Spatial domain methods which can compute the phase from a single fringe pattern in the spatial domain. (c) Frequency domain methods which uses some kind of transformation to the frequency domain to compute the phase. In this category, the Fourier and Wavelet transforms are the most common mathematical tools to carry out the task.

Apart from the phase recovery, there are other important steps in fringe pattern processing. For example, many times the fringe patterns are corrupted by noise, such as the case of the electronic speckle pattern interferometry. Then, fringe image enhancement by means of low-pass filtering is usually required. Owing that most algorithms to retrieve the phase from a fringe pattern give the phase wrapped in the interval $\left[-\pi ,\pi \right)$, other important step is the well-known phase-unwrapping process [6, 7]. In the field of fringe image enhancement, such as fringe image denoising or phase denoising, there has been a wide research activity in the last years. Researchers have realized that improving the quality of fringe images and wrapped phase fields is of main relevance for a successful phase recovery or phase unwrapping. However, enhancing fringe images or wrapped phase fields has resulted to be a task that must be realized in a special manner, so that ordinary techniques for image enhancement are not always adequate. Owing that frequencies of fringes and noise usually overlap and normally cannot be properly separated, common filters for image processing have blurring effects on fringe features, especially for patterns with high density fringes. For these cases, the use of anisotropic filters is a better way for removing noise without the harmful blurring effects.

In the fields of fringe pattern denoising and wrapped phase map denoising, there have been many proposals to realize these tasks. Some of the first contributions in this field were mainly based on convolution filters using different kinds of anisotropic filtering masks [8, 9, 10, 11, 12]. Other set of the main contributions in the last years is based on the variational calculus approach by solving partial differential equations [13, 14, 15, 16, 17, 18], and by means of the regularization theory [19, 20]. The use of the Fourier transform for fringe or phase map denoising has also been proposed in [21, 22] (Localized Fourier transform filter and windowed Fourier transform, respectively). There have been other proposals that used different methodologies such as coherence enhancing diffusion [23], image decomposition [24], and multivariate empirical mode decomposition [25]. The great disadvantage of already reported methods for fringe and phase map denoising is that they require the previous estimation of the so-called fringe orientation which, as it uses the computation of the image gradient, could be an inaccurate procedure in the presence of noise and low modulation of fringes. This is not the case for the Fourier-based methods [21, 22]; however, as in the case of the Windowed Fourier transform technique, several parameters have to be adjusted depending on the particular image and it may require a long processing time.

In the field of phase recovery from fringe images, there have been a lot of researches along the last decades. For the case of phase-shifting algorithms, outstanding summaries of them can be found in [2, 26]. For the case of spatial and frequency domain methods from a single pattern image, two of the most popular techniques are the well-known Fourier Transform method reported by Takeda et al. [27] and the Synchronous detection method [28]. Other methods that use the regularization theory were also proposed [29, 30]. However, although these methods are efficient and easy to implement, they are limited to be used in fringe images with frequency carrier, which just in few experimental situations these kinds of images can be obtained. In most cases, experimental conditions in optical measurement techniques yields fringe images without a dominant frequency (i.e., closed fringes) which becomes the phase recovery problem difficult, therefore more complicated algorithms must be used. One of the first proposals for phase demodulation from single closed fringe images was reported by Kreis using a Fourier based approach [31]. In the last decade of the twentieth century, it was a boom in the research of closed fringe images, specially using the regularization theory. The Regularized phase-tracking technique was reported by Servín et al. [32]. Marroquín et al. reported the regularized adaptive quadrature filters [33] and the regularization method that uses the local orientation of fringes [34]. At the beginning of this century, Larkin et al. proposed the spiral-phase quadrature transform [35] and Servín et al. reported the General n-dimensional quadrature transform [36]. Also, we proposed the orientational vector-field-regularized estimator to demodulate closed fringe images [37].

As will be shown, closed fringe and wrapped-phase images have certain characteristics that make them to be treated in a special manner. First, it is common that this kind of images present structures with high anisotropy at the same time that many frequencies are dispersed over the entire image. For these reasons, in most situations, the use of linear-translation-spatial (LTI) filters, which are spatially invariant and independent of image content, do not give proper results. Furthermore, owing that the Fourier transform is a global operation, this technique is not always suitable for accurately model the local characteristics of closed fringe images.

It is widely known that the wavelet transform is a powerful tool that provides local, sparse, and decorrelated multiresolution analysis of signals. In the last years, 2D wavelets have been used for image analysis as a proper alternative to the weakness of LTI filters and linear transforms as the Fourier one. In particular, it has been shown that 1-D and 2D continuous wavelet transform (CWT) using Gabor atoms is a natural choice for proper analyses of fringe images. This kind of analysis has been used for fringe pattern denoising and fringe pattern demodulation showing several advantages, for example in laser plasma interferometry [38], in shadow moiré [39, 40, 41], in profilometry [42, 43, 44], in speckle interferometry [45], in digital holography [46], and other optical measurement techniques [47, 48, 49, 50, 51, 52, 53, 54, 55].

In this chapter, the theoretical basis of fringe pattern image formation and processing is described. Also, in general, the theory and advantages of the 2D continuous wavelet transform (CWT) for fringe pattern processing is described. We also explain some of the main applications in fringe pattern processing, such as phase recovery and wrapped phase map denoising, showing some examples of applications in different optical measurement techniques.

2. Digital fringe patterns

2.1. Elements of digital fringe image processing systems

Often, a digital fringe image processing system is represented by a sequence of devices, which typically starts with an imaging system that observes the target, a digitizer system which samples and quantizes the analog information acquired by the imaging system, a digital storage device, a digital computer that process the information, and finally, a displaying system to visualize the acquired and processed information (Figure 1).

Figure 1.

Typical sequence in a digital fringe pattern image processing system.

A typical imaging system is composed by an objective lens to form images in a photosensitive plane which is commonly a CCD (charge couple devices) array.

2.2. Fringe image formation

Fringe pattern images are present in several kinds of optical tests for the measurement of different physical quantities. Such tests are examples for the quality measurement of optical devices using optical interferometry, photoelasticity for stress analysis, or electronic speckle pattern interferometry (ESPI) for the measurement of mechanical properties of materials. The interference phenomena are usually used in many optical methods of measurement. We now describe a classical way to form a fringe pattern image using the two-wave interference.

Two-wave interference can be generated by means of several types of interferometers, and the interferograms or fringe patterns are produced by superimposing two wavefronts. An interferometer can accurately measure deformations of the wavefront of the order of the wavelength. Considering two mutually coherent monochromatic waves, as depicted in Figure 2, $W\left(x,y\right)$ represents the wavefront shape under study (i.e., the wave that contains the information of the physical quantity to be measured). The sum of their complex amplitudes can be represented as

$$E\left(x,y\right)=A_1\left(x,y\right)e^{\mathit{ikW}\left(x,y\right)}+A_2\left(x,y\right)e^{\mathit{ikx}\text{sin}\theta },$$E1

where $A_1$ and $A_2$ are the amplitudes of the wavefront under test and the reference wavefront (a flat wavefront), respectively, and $k=\frac{2\pi }{\lambda }$, being $\lambda$ the wavelength.

Figure 2.

Interference of two wavefronts. Solid line represents the wavefront under test and dashed line represents the reference wavefront.

The irradiance at a given plane perpendicular to $z$-axis is then represented as

$$\begin{array}{ll}I\left(x,y\right)& =E\left(x,y\right)E^{\ast }\left(x,y\right)\\ & =A_1^2\left(x,y\right)+A_2^2\left(x,y\right)+2A_1\left(x,y\right)A_2\left(x,y\right)\text{cos}\left[\mathit{kx}\text{sin}\theta +\mathit{kW}\left(x,y\right)\right].\end{array}$$E2

For simplicity, Eq. (2) is usually written in a general form as:

$$I\left(x,y\right)=a\left(x,y\right)+b\left(x,y\right)\text{cos}\left[u_{0}x+\varphi \left(x,y\right)\right],$$E3

where $a\left(x,y\right)$ and $b\left(x,y\right)$ are commonly called the background illumination and the amplitude modulation, respectively. The term $u_0=k\text{sin}\theta$ is the fringe carrier frequency and $\varphi \left(x,y\right)=\mathit{kW}\left(x,y\right)$ is the phase to be recovered from the fringe pattern image. It must be noted that if the reference wavefront is perpendicular to $z$-axis (i.e., $\theta =0$), the fringe carrier frequency is removed and Eq. (3) is simplified:

$$I\left(x,y\right)=a\left(x,y\right)+b\left(x,y\right)\text{cos}\left[\varphi \left(x,y\right)\right].$$E4

Equations (3) and (4) represent the mathematical expressions of fringe pattern images with and without fringe carrier frequency, respectively. Examples of these kinds of fringe images are shown in Figure 3.

Figure 3.

Examples of simulated fringe pattern images with (a) and without (b) fringe carrier frequency. The phase of modulation $\varphi \left(x,y\right)$ (c) is the same for both fringe images (phase shown wrapped and codified in gray levels).

3. Fringe pattern processing

3.1. Phase-shifting methods for phase recovery

One of the most popular methods for phase recovery is the well-known phase-shifting. This method requires a set of phase-shifted fringe patterns which are experimentally obtained in different ways depending on the optical measurement technique. For example, in interferometry the phase shifting is realized by moving some mirrors in the optical interferometer. The set of $N$ phase-shifted fringe patterns is defined as

$$I_n\left(x,y\right)=a\left(x,y\right)+b\left(x,y\right)\text{cos}\left[\varphi \left(x,y\right)+{\alpha }_n\right]n=1,2,\dots ,N.$$E5

The pointwise solution for $\varphi \left(x,y\right)$ from the non-linear system of equations is obtained by using the last-squares approach (see [2] for details):

$$W\left\{\varphi \left(x,y\right)\right\}={\text{tan}}^{-1}\left(-\frac{{\sum }_{n=1}^N{I}_n\text{sin}\left({\alpha }_n\right)}{{\sum }_{n=1}^N{I}_n\text{cos}\left({\alpha }_n\right)}\right)\in \left[-\pi ,\pi \right),$$E6

where $W$ is the wrapping operator such that $W\left\{\varphi \left(x,y\right)\right\}\in \left[-\pi ,\pi \right)$. Several algorithms can be used that require three, four, up to eight images.

3.2. Phase recovery from single fringe patterns with carrier

As previously mentioned, processing fringe patterns with fringe carrier frequency may be simple to carry out. The key point in the demodulation of fringe patterns with carrier is that the total phase function $u_{0}x+\varphi \left(x,y\right)$ represents the addition of an inclined phase plane $u_{0}x$ plus the target phase $\varphi \left(x,y\right)$. In this case, a monotonically increasing (or decreasing) phase function has to be recovered. If we analyze the Fourier spectrum of Eq. (3), for a proper separation between spectral lobes in the Fourier space, the following inequality must be complied:

$$\text{max}\left\{\Vert \nabla \varphi \Vert \right\}<\Vert u_0\Vert .$$E7

The analytic signal $g\left(x,y\right)$ to recover the phase $\varphi \left(x,y\right)$ can be computed with the Fourier transform method [27], which can expressed as

$$g\left(x,y\right)={\mathcal{F}}^{-1}\left\{H\left(u,v\right)\mathcal{F}\left\{I\left(x,y\right)\right\}\right\}=e^{i2\pi \left[u_{0}x+\varphi \left(x,y\right)\right]},$$E8

where $H\left(u,v\right)$ is a filter in the Fourier domain centered at the frequency $u_0$, $u$ the frequency variable along $x$ direction, and $v$ the frequency variable along $y$ direction. Finally, the wrapped phase is computed with

$$W\left\{\varphi \left(x,y\right)\right\}={\text{tan}}^{-1}\left(\frac{\text{Real}\left\{g\left(x,y\right)e^{-i2\pi u_0}\right\}}{\text{Imag}\left\{g\left(x,y\right)e^{-i2\pi u_0}\right\}}\right)\in \left[-\pi ,\pi \right).$$E9

Other technique to compute the phase from a carrier frequency fringe pattern is the synchronous detection technique [28], which is realized in the spatial domain. Using the complex notation, in this case, the analytic function $g\left(x,y\right)$ can be computed with

$$g\left(x,y\right)=h\left(x,y\right)\ast \left[I\left(x,y\right)e^{i2\pi u_0}\right]=e^{i2\mathit{\pi \varphi }\left(x,y\right)},$$E10

where ${}^{\ast }$ represents the convolution operator and $h\left(x,y\right)$ a low-pass convolution filter in the spatial domain. The wrapped phase can be computed with

$$W\left\{\varphi \left(x,y\right)\right\}={\text{tan}}^{-1}\left(\frac{\text{Real}\left\{g\left(x,y\right)\right\}}{\text{Imag}\left\{g\left(x,y\right)\right\}}\right)\in \left[-\pi ,\pi \right).$$E11

3.3. Phase recovery from single fringe patterns without carrier

As described in [34, 35, 36, 37], for the case in which $u_0=0$, the previous computation of the fringe direction is necessary to compute the analytic function $g\left(x,y\right)$, for example, using the quadrature transform [36]:

$$\text{Imag}\left\{g\left(x,y\right)\right\}=\text{sin}\left[\varphi \left(x,y\right)\right]={\mathbf{n}}_{\varphi }\left(x,y\right)\cdot \frac{\nabla I_n\left(x,y\right)}{\Vert \nabla \varphi \left(x,y\right)\Vert },$$E12

where $I_n\left(x,y\right)=\text{cos}\left[\varphi \left(x,y\right)\right]=\text{Real}\left\{g\left(x,y\right)\right\}$ is a normalized version of $I\left(x,y\right)$, and ${\mathbf{n}}_{\varphi }$ is the unit vector normal to the corresponding isophase contour, which points to the direction of $\nabla \varphi \left(x,y\right)$. It is well known that the computation of ${\mathbf{n}}_{\varphi }$ is by far the most difficult problem to compute the phase using this method.

Also, the modulo-$2\pi$ fringe orientation angle $\alpha \left(x,y\right)$ can be used to compute the quadrature fringe pattern by means of the spiral-phase signum function $S\left(u,v\right)$ in the Fourier domain [35]:

$$\text{Imag}\left\{g\left(x,y\right)\right\}=\text{sin}\left[\varphi \left(x,y\right)\right]=-{\mathit{ie}}^{-\mathit{i\alpha }\left(x,y\right)}{\mathcal{F}}^{-1}\left\{S\left(u,v\right)\mathcal{F}\left\{I_n\left(x,y\right)\right\}\right\},$$E13

where

$$S\left(u,v\right)=\frac{u+\mathit{iv}}{\sqrt{u^2+v^2}},$$E14

and $i=\sqrt{-1}$. However, the most difficult problem in this method is the computation of $\alpha \left(x,y\right)$. It can be deduced that Eqs. (12) and (13) are closely related because

$$\alpha \left(x,y\right)=\text{angle}\left\{{\mathbf{n}}_{\varphi }\left(x,y\right)\right\}\in \left(0,2\pi \right].$$E15

3.4. Wrapped phase maps denoising

The unwrapping process can be, in many cases, a difficult task due to phase inconsistencies or noise. In order to understand the phase unwrapping problem of noisy phase maps, we define the wrapped and the unwrapped phase as $\psi \left(x,y\right)$ and $\varphi \left(x,y\right)$ respectively. As it is known that $\psi \left(x,y\right)\in \left[-\pi ,\pi \right)$, the following relation can be established:

$$\psi \left(x,y\right)=\varphi \left(x,y\right)+2\mathit{\pi k}\left(x,y\right),$$E16

where $k\left(x,y\right)$ is a field of integers such that $\psi \left(x,y\right)\in \left[-\pi ,\pi \right)$. The wrapped phase-difference vector field $\mathrm{\Delta }\psi \left(x,y\right)$ which can be computed from the wrapped phase map, is defined as

$$\mathrm{\Delta }\psi \left(x,y\right)=\left[\psi \left(x,y\right)-\psi \left(x-1,y\right),\psi \left(x,y\right)-\psi \left(x,y-1\right)\right],$$E17

where $\left(x-1,y\right)$ and $\left(x,y-1\right)$ are contiguous horizontal and vertical sites, respectively. In a similar manner, we can also define the unwrapped phase-difference field:

$$\mathrm{\Delta }\varphi \left(x,y\right)=\left[\varphi \left(x,y\right)-\varphi \left(x-1,y\right),\varphi \left(x,y\right)-\varphi \left(x,y-1\right)\right].$$E18

It can be deduced that the problem of the recovery of $\varphi$ from $\psi$ can be properly solved if the sampling theorem is reached, that is, if the distance between two fringes is more than two pixels (the phase difference between two fringes is $2\pi$). In phase terms, the sampling theorem is reached if the phase difference between two pixels is less than $\pi$ or, in general

$$\Vert \mathrm{\Delta }\varphi \Vert <\pi ,\forall \left(x,y\right).$$E19

If this condition is satisfied, the following relation can be established:

$$\mathrm{\Delta }\varphi =W\left\{\mathrm{\Delta }\psi \right\}=\left[{\psi }_x,{\psi }_y\right],$$E20

where

$${\psi }_x=W\left\{\psi \left(x,y\right)-\psi \left(x-1,y\right)\right\}\text{and}{\psi }_x=W\left\{\psi \left(x,y\right)-\psi \left(x,y-1\right)\right\}.$$E21

Note that $W\left\{\mathrm{\Delta }\psi \right\}$ (the wrapped phase differences) can be obtained from the observed wrapped phase field $\psi$. Then, the unwrapped phase $\varphi$ can be achieved by two-dimensional integration of the vector field $W\left\{\mathrm{\Delta }\psi \right\}$.

A simple way to compute the unwrapped phase $\varphi$ from the wrapped one $\psi$ is by means of minimizing the cost function

$$U\left(\varphi \right)=\sum _{\left(x,y\right)\in L}\left\{{\left[{\psi }_x\left(x,y\right)-\left(\varphi \left(x,y\right)-\varphi \left(x-1,y\right)\right)\right]}^2+{\left[{\psi }_y\left(x,y\right)-\left(\varphi \left(x,y\right)-\varphi \left(x,y-1\right)\right)\right]}^2\right\},$$E22

where $L$ is the set of valid pixels in the image. Unfortunately, in most cases noise is present, therefore, inequality (19) is not always satisfied and the integration does not provide proper results. Therefore, denoising wrapped phase maps is a fundamental step before the phase unwrapping process.

4. The 2D continuous wavelet transform for processing fringe patterns

It is clear that the phase demodulation of fringe images with carrier may be easily realized. Owing that, in this case, the fringe image may represent a quasi-stationary signal along the direction of the frequency carrier, the use of classical linear operators such as the Fourier transform may be adequate. It works well mainly for few components in the frequency domain (i.e., for narrow spectrums); however, this is not the case for many signals in the real world. This dependence is a serious weakness mainly in two aspects: the degree of automation and the accuracy of the method specially when fringes produce spread spectrums due to localized variations or phase transients. Additionally, in the case of closed fringes there may be a wide range of frequencies in all directions. Then, evidently standard Fourier analysis is inadequate for treating with this kind of images because it represents signals with a linear superposition of sine waves with “infinite” extension. For this reason, an image with closed fringes should be represented with localized components characterizing the frequency, shifting, and orientation. A powerful mathematical tool for signal description that has been developed in the last decades is the wavelet analysis. Fortunately, for our purposes, a key characteristic of this type of analysis is the finely detailed description of frequency or phase of signals. In consequence, it can have a good performance especially with fringes that produce spread spectrums. Additionally, one of the main advantages using wavelets compared with standard techniques is its high capability to deal with noise. In particular, the 2D continuous wavelet transform have recently been proposed for the processing of interferometric images. Advantages of denoising and demodulation of interferograms using the 2D CWT has been discussed in [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55].

Considering an interferometric image (an interferogram or a wrapped-phase field) $G\left(r\right)$, where $r=\left(x,y\right)\in {\mathbb{R}}^2$, its 2D CWT decomposition can be defined as

$$G_{\mathcal{W}}\left(s,\theta ,\eta \right)=\mathcal{W}\left\{G\left(r\right)\right\}={\int }_{{\mathbb{R}}^2}G\left(r\right){\phi }_{s,\theta ,\eta }^{\ast }\left(r\right)\mathit{dr}.$$E23

In Eq. (23), $\phi$ represents the 2D mother wavelet and ${}^{\ast }$ indicates the complex conjugated. The variable $s\in {\mathbb{R}}^2$ represents the shift, $\theta \in \left[0,2\pi \right)$ the rotation angle, and $\eta$ the scaling factor. It has been shown that a proper mother wavelet for processing interferometric images is the 2D Gabor wavelet (see Figure 4). The mathematical representation of this kind of wavelet can be defined as

$${\phi }_{s,\theta ,\eta }\left(r\right)=\text{exp}\left[-\pi \frac{\parallel r-s{\parallel }^2}{\eta }\right]\times \text{exp}\left[i2\pi \frac{\nu }{\eta }\left(\left(r-s\right)\cdot \mathrm{\Theta }\right)\right],$$E24

where $\mathrm{\Theta }=\left(\text{cos}\theta ,\text{sin}\theta \right)$, $\left(\cdot \right)$ represents the dot product, and $\nu \in \mathbb{R}$ is the frequency variable.

Figure 4.

Example of a 2D Gabor wavelet. (a) Real part and (b) imaginary part.

Figure 5 shows that the 2D CWT is performed along different directions and frequencies.

Figure 5.

Frequency localization of the 2D wavelets in the Fourier domain ($f=\frac{\nu }{\eta }$).

4.1. Phase recovery with the 2D CWT

Owing that fringe pattern images with closed fringes generally contain elements with high anisotropy and sparse frequency components, the phase recovery is a complex procedure. Compounding the problem, the presence of noise makes the process even more complicated because noise and fringes are mixed in the Fourier domain.

Also, it has been shown that a single fringe pattern without carrier frequency, is not easy to deal with. Owing to ambiguities in the image formation process, a main drawback analyzing them is that several solutions of the phase function can satisfy the original observed image. Therefore, it is necessary to restrict the solution space of $\varphi$ in Eq. (4). Fortunately, as in most practical cases the phase to be recovered is continuous, the algorithm to process the fringe pattern usually seeks for a continuous phase function. However, the recovery of the continuous phase function is not a simple task to carry out as occur with fringe patterns with carrier frequency. It can be observed that the phase gradient represents the local frequencies of the fringe pattern in the $x$ and $y$ directions; however, the sign of $\nabla \varphi$ is ambiguous because negative and positive frequencies are mixed in the Fourier domain.

The following is a general description of the phase recovery method using the 2D CWT. First, it is necessary to consider a normalized version of the fringe pattern. The normalization procedure can be carried out using the method proposed in [56]. Consider we represent the normalized fringe pattern in complex form:

$$G\left(r\right)=\text{cos}\left[\varphi \left(r\right)\right]=\frac{\text{exp}\left[\mathit{i\varphi }\left(r\right)\right]}{2}+\frac{\text{exp}\left[-\mathit{i\varphi }\left(r\right)\right]}{2}.$$E25

In this particular case, the 2D CWT of $G\left(r\right)$ is

$$\begin{array}{l}\mathcal{W}\left\{G\left(r\right)\right\}={\int }_{{\mathbb{R}}^2}\frac{\text{exp}\left[\mathit{i\varphi }\left(r\right)\right]}{2}\text{exp}\left[-\pi \frac{\parallel r-s{\parallel }^2}{\eta }\right]\times \text{exp}\left[-i2\pi \frac{\nu }{\eta }\left(\left(r-s\right)\cdot \mathrm{\Theta }\right)\right]\mathit{dx}\\ +{\int }_{{\mathbb{R}}^2}\frac{\text{exp}\left[-\mathit{i\varphi }\left(r\right)\right]}{2}\text{exp}\left[-\pi \frac{\parallel r-s{\parallel }^2}{\eta }\right]\times \text{exp}\left[-i2\pi \frac{\nu }{\eta }\left(\left(r-s\right)\cdot \mathrm{\Theta }\right)\right]\mathit{dx}.\end{array}$$E26

Note that $\mathcal{W}\left\{G\left(r\right)\right\}$ represents a four-dimensional function depending on $x,y,\eta$, and $\theta$. The process to recover the phase $\varphi \left(r\right)$ using the 2D CWT consists on realizing the well-known ridge detection. To understand the phase recovery from the ridge detection, first it is necessary to know the meaning of Eq. (26). To do so, let $\tilde{r}=r-s$ and ${\nu }_{\theta }=\frac{\nu }{\eta }\left(\text{cos}\theta ,\text{sin}\theta \right)$, where ${\nu }_{\theta }\in {\mathbb{R}}^2$. Using Taylor’s expansion we know that

$$\varphi \left(\tilde{r}+s\right)\approx \varphi \left(s\right)+\nabla \varphi \left(s\right)\cdot \tilde{r}.$$E27

Then, we can now rewrite Eq. (26) as

$$\begin{array}{l}\mathcal{W}\left\{G\left(r\right)\right\}\approx \frac{\text{exp}\left[\mathit{i\varphi }\left(s\right)\right]}{2}{\int }_{{\mathbb{R}}^2}\text{exp}\left[i\left(\nabla \varphi \left(s\right)\cdot {\displaystyle \tilde{r}}\right)\right]\times \text{exp}\left[-\pi \frac{\parallel {\displaystyle \tilde{r}}{\parallel }^2}{\eta }\right]\text{exp}\left[-i2\pi \left(\tilde{r}\cdot {\nu }_{\theta }\right)\right]d\tilde{r}\\ +\frac{\text{exp}\left[-\mathit{i\varphi }\left(s\right)\right]}{2}{\int }_{{\mathbb{R}}^2}\text{exp}\left[-i\left(\nabla \varphi \left(s\right)\cdot {\displaystyle \tilde{r}}\right)\right]\times \text{exp}\left[-\pi \frac{\parallel {\displaystyle \tilde{r}}{\parallel }^2}{\eta }\right]\text{exp}\left[-i2\pi \left(\tilde{r}\cdot {\nu }_{\theta }\right)\right]d\tilde{r},\end{array}$$E28

or, which is the same

$$\begin{array}{l}\mathcal{W}\left\{G\left(r\right)\right\}\approx \frac{\text{exp}\left[\mathit{i\varphi }\left(s\right)\right]}{2}\mathcal{F}\left\{\text{exp}\left[i\left(\nabla \varphi \left(s\right)\cdot {\displaystyle \tilde{r}}\right)\right]\times \text{exp}\left[-\pi \frac{\parallel {\displaystyle \tilde{r}}{\parallel }^2}{\eta }\right]\right\}\\ +\frac{\text{exp}\left[-\mathit{i\varphi }\left(s\right)\right]}{2}\mathcal{F}\left\{\text{exp}\left[-i\left(\nabla \varphi \left(s\right)\cdot {\displaystyle \tilde{r}}\right)\right]\times \text{exp}\left[-\pi \frac{\parallel {\displaystyle \tilde{r}}{\parallel }^2}{\eta }\right]\right\}.\end{array}$$E29

The two terms in (29) contains Fourier transforms of complex periodic functions of frequencies $\nabla \varphi \left(s\right)/2\pi$ and $-\nabla \varphi \left(s\right)/2\pi$. Then, applying the Fourier’s similarity and modulation theorems this last equation can be finally written as

$$\begin{array}{l}\mathcal{W}\left\{G\left(r\right)\right\}\approx \eta \frac{\text{exp}\left[\mathit{i\varphi }\left(s\right)\right]}{2}\text{exp}\left[-\mathit{\eta \pi }{‖{\nu }_{\theta }-\frac{\nabla \varphi \left(s\right)}{2\pi }‖}^2\right]\\ +\eta \frac{\text{exp}\left[-\mathit{i\varphi }\left(s\right)\right]}{2}\text{exp}\left[-\mathit{\eta \pi }{‖{\nu }_{\theta }+\frac{\nabla \varphi \left(s\right)}{2\pi }‖}^2\right].\end{array}$$E30

In this case, ${\nu }_{\theta }$ is the two-dimensional frequency variable. Note that for a fixed $s$, $\mathcal{W}\left\{G\left(r\right)\right\}$ represents two Gaussian filters in the Fourier domain localized at polar coordinates $\left(\frac{\nu }{\eta },\theta \right)$. It can also be visualized as an orientation and frequency decomposition of the fringe pattern.

To detect the analytic function and consequently compute the phase $\varphi \left(s\right)$ at a given pixel $s$ (i.e., the ridge detection), we can choice one of two possibilities: at ${\nu }_{\theta }=\frac{\nabla \varphi \left(s\right)}{2\pi }$ or ${\nu }_{\theta }=-\frac{\nabla \varphi \left(s\right)}{2\pi }$. Owing that the sign of the phase gradient cannot be determined from the image intensity, there exists a sign ambiguity of the phase in the $\theta -\eta$ map. In Figure 6, it can be observed that in this situation, there are two maximum in each $\theta -\eta$ map. Also, it can be deduced that the magnitude of the coefficients map is periodic with respect to $\theta$ with period $\pi$. To solve the problem of sign ambiguity, Ma et al. [48] proposed a phase determination rule according to the phase distribution continuity. Also, Villa et al. [55] proposed a sliding 2D CWT method that assumes that the phase is continuous and smoothly varying, in this way, the ridge detection is realized assuming that the coefficient maps are similar in adjacent pixels, reducing the processing time too.

Figure 6.

(a) Example of noisy simulated fringe pattern. The square indicates a region around a pixel $s$ where the phase is estimated. (b) Magnitude of the $\theta -\eta$ map at the pixel $s$, codified in gray levels. Horizontal direction represents the rotation angle while the vertical direction represents the scale. The two white regions represent the two terms in Eq. (30).

Once detected the ridge $\mathcal{W}{\left\{G\left(r\right)\right\}}_{\mathit{\text{ridge}}}$ that represents a 2D function, the wrapped phase can be computed with

$$W\left\{\varphi \left(r\right)\right\}={\text{tan}}^{-1}\left(\frac{\text{Real}\left\{\mathcal{W}{\left\{G\left(r\right)\right\}}_{\mathit{\text{ridge}}}\right\}}{\text{Imag}\left\{\mathcal{W}{\left\{G\left(r\right)\right\}}_{\mathit{\text{ridge}}}\right\}}\right).$$E31

Figures 7 and 8 show examples of fringe pattern phase recovery using the 2D CWT method reported in [55]. It is important to remark that this method is highly robust against noise.

Figure 7.

Example of the 2D CWT method applied to phase recovery. (a) Synthetic noisy fringe pattern. (b) Recovered phase.

Figure 8.

Example of the 2D CWT method applied to phase recovery. (a) Experimentally obtained moiré fringe pattern. (b) Recovered phase.

A big advantage of using the 2D CWT method to compute the phase from fringe patterns without carrier is that the sign ambiguity of $\nabla \varphi$ can be easily solved, for example, with the method reported in [55]. The key idea of the method is the assumption that the phase $\varphi$ is smooth; in other words, the fringe frequency and fringe orientation are very similar in neighbor pixels, hence the ridge detection at each $\theta -\eta$ map is simplified registering the previous computation of neighbor pixels.

4.2. The 2D CWT for wrapped phase maps denoising

Other of the most relevant tasks in fringe pattern processing is the wrapped phase maps denoising. Owing that the phase unwrapping is a key step in fringe pattern processing for optical measurement techniques, the previous denoising of the wrapped phase is crucial for a proper measurement. Several optical measurement techniques, such as the electronic speckle pattern interferometry, use different phase recovery methods, inherently produces highly noisy wrapped phase maps. In these situations, the phase map denoising is a crucial pre-process for a successful phase unwrapping. Considering the problem of denoising wrapped phase maps, the drawback is that owing to $2\pi$ phase jumps of the wrapped phase $\psi$, direct application of any kind of filter is not always a proper procedure to solve it. For example, the application of a simple mean filter may smear out the phase jumps. In order to avoid this drawback, the wrapped phase filtering must be realized computing the following complex function:

$$G\left(r\right)=\text{exp}\left[\mathit{i\psi }\left(r\right)\right],$$E32

where $i=\sqrt{-1}$. As both imaginary and real parts are continuous functions, we can properly apply a filter over $G\left(r\right)$, and the argument of the filtered complex signal will contain the denoised phase map. Again, substituting (32) in (23), we now obtain

$$\mathcal{W}\left\{G\left(r\right)\right\}={\int }_{{\mathbb{R}}^2}\text{exp}\left[\mathit{i\psi }\left(r\right)\right]\text{exp}\left[-\pi \frac{\parallel r-s{\parallel }^2}{\eta }\right]\times \text{exp}\left[-i2\pi \frac{\nu }{\eta }\left(\left(r-s\right)\cdot \mathrm{\Theta }\right)\right]\mathit{dx}.$$E33

Following the same reasoning to obtain Eq. (30), for this case, we obtain:

$$\mathcal{W}\left\{G\left(r\right)\right\}\approx \eta \text{exp}\left[\mathit{i\psi }\left(s\right)\right]\text{exp}\left[-\mathit{\eta \pi }{‖{\nu }_{\theta }-\frac{\nabla \psi \left(s\right)}{2\pi }‖}^2\right].$$E34

The difference of this equation with the result shown in Eq. (30) is that at each $\theta -\eta$ map, there is only one maximum: at ${\nu }_{\theta }=\frac{\nabla \psi \left(s\right)}{2\pi }$ (see Figure 9). Thus, in this case, the ridge detection is simpler and the filtered wrapped phase map ${\psi }_f\left(r\right)$ can be computed with

$${\psi }_f\left(r\right)={\text{tan}}^{-1}\left(\frac{\text{Real}\left\{\mathcal{W}{\left\{G\left(r\right)\right\}}_{\mathit{\text{ridge}}}\right\}}{\text{Imag}\left\{\mathcal{W}{\left\{G\left(r\right)\right\}}_{\mathit{\text{ridge}}}\right\}}\right)\in \left[-\pi ,\pi \right).$$E35

Figure 9.

(a) Zoom of a small square region in a noisy wrapped phase map (around some pixel $s$). (b) Magnitude of the $\theta -\eta$ map at the pixel $s$, codified in gray levels. Horizontal direction represents the rotation angle while the vertical direction represents the scale.

Figures 10 and 11 are examples of the results applying the 2D CWT in wrapped phase map denoising. Note the outstanding performance removing the structures due to the gratings in the experimentally obtained wrapped phase map with moire deflectometry (Figure 11).

Figure 10.

(a) Simulated noisy wrapped phase map. (b) Filtered wrapped phase map.

Figure 11.

(a) Experimentally obtained moiré noisy wrapped phase map. (b) Filtered wrapped phase map.

The key step in the 2D CWT method for phase map denoising is the ridge detection. In this way, all the coefficients in the $\theta -\eta$ map contributed by the noise and spurious information are removed. A comparison of the performance of this method compared with the windowed Fourier transform method [22] and the localized Fourier transform method [21] is shown in Table 1. In this case, the normalized-mean-square-error (NMSE) was used as the metric applied over a synthetic noisy phase map $\psi$ (Figure 10). Although the performance against noise of the WFT is better that the 2D CWT method, this last is much simpler to implement, as discussed in [53].

$$\text{NMSE}=\frac{\parallel \psi -{\psi }_f{\parallel }^2}{\parallel \psi {\parallel }^2}.$$E36

2D-CWT	WFT	LFT
0.0692	0.0521	0.0747

Table 1.

Performance comparison of the 2D CWT, WFT, and LFT methods, using the NMSE.

5. Conclusions

It can be obviously deduced that often fringe patterns contain elements with high anisotropy, sparse frequency components, and noise, which makes the processing of this kind of images by means of classical LTI methods inadequate. Several authors have shown that the use of multiresolution analysis by means of the 2D CWT for processing fringe patterns has resulted a proper and interesting alternative for this task. The 2D CWT methods present some attractive advantages compared with other commonly used techniques. (1) The use of the Gabor mother wavelet for processing this kind of images is a natural choice to model them, as can be obviously deduced analyzing the physical theory of fringe image formation. (2) In most classical methods for processing fringe images, the previous estimation of the fringe direction or orientation is a must, especially for fringe patterns without a fringe carrier frequency. Owing that the multiresolution analysis using the 2D CWT methods models the image by means of the angle $\theta$, fringe direction or orientation is inherently computed through the ridge detection. (3) As the 2D CWT methods models the interferograms by means of scale and orientation, all spurious information and noise contributing in the $\theta -\eta$ map is efficiently removed through the ridge detection, resulting a powerful tool to remove the noise.