## Abstract

Single-shot crossed-type fringe pattern processing and analysis method called Hilbert-Huang grating interferometry (HHGI) is proposed. It consist of three main procedures: (1) crossed pattern is resolved into two fringe families using novel orthogonal empirical mode decomposition approach, (2) separated fringe sets are filtered using modified automatic selective reconstruction aided by enhanced fast empirical mode decomposition and mutual information detrending, and (3) Hilbert spiral transform is employed for fringe phase demodulation. Numerical and experimental studies corroborate the validity, versatility and robustness of the proposed HHGI technique. It can be successfully applied to multiplicative and additive type crossed patterns with sinusoidal and binary orthogonal component structures. Efficient adaptive filtering enables successful fast processing and analysis of complex and defected patterns.

© 2013 Optical Society of America

## 1. Introduction

Optical interferometry represents one of the most powerful tools in optical metrology. It is indispensable for quantitative studies of diverse scientific phenomena and engineering problems. As a full-field technique it enables simultaneous acquisition and parallel processing of experimental data for evaluating both macro and microscale object information in static and dynamic events. The measurand is usually encoded in the shape and position of fringes, i.e., the departure from the interferogram uniform intensity distribution (null-fringe detection mode) or the reference fringe departure from straightness (finite fringe detection mode). Quantitative analysis of fringe patterns is performed using computer-aided automatic fringe pattern analysis (AFPA) methods (see, for example [1,2],). Although AFPA has over 30 years of history, it is still under dynamic development dictated by new needs and challenges.

Two-beam interferometry provides the object phase information by employing three basic approaches:

- a) Interference of object and reference beams (reference beam interferometry, object phase maps are obtained [3]);
- b) Interference of object beam with its replica (shearing interferometry; first derivative of object phase is obtained under the assumptions of small displacement between the two beams and slow object phase variations [3]);

Interferometers use beam-splitting and recombining elements that can be made of bulk, diffraction and polarization optics. Diffraction gratings are especially attractive because they can be applied in a wide spectral range. Both linear and crossed gratings are employed.

This paper is devoted to automatic processing and analysis of interference patterns obtained in crossed (two-dimensional) grating interferometers with simultaneous, single-frame recording of information encoded in mutually orthogonal directions. In lateral shearing interferometry the representative example is Talbot interferometry based on the self-imaging phenomenon and developed independently by Yokozeki and Suzuki [6,7] and Lohmann and Silva [8–10]. Ten years later a similar technique was described by Kafri under the name of moiré deflectometry [11]. The diffraction approach was used by the first authors to explain the system performance, whereas Kafri’s description was based on ray optics. The first approach is more universal, taking into consideration various effects involved. The two treatments give equivalent results as signaled by Patorski [12] although Keren and Kafri [13] and their collegues [14] first mentioned the difference between the cases of placing the phase object in front of and behind the first grating. The origin of Talbot interferometry and moiré deflectometry can be traced to the paper of Oster and associates [15, 16]. In the visible spectrum range, where both gratings are mostly of amplitude type, an exemplary application of 2D Talbot interferometry (deflectometry), among many others [16], deals with the ophthalmic lens power mapping [17–21]. X-ray 2D Talbot interferometry represents the recent exciting and challenging topic of studies [22–29]. Here the first grating is a binary periodic phase structure and its best intensity contrast Fresnel image (located in-between the self-images) is detected by the amplitude type binary grating. Two implementations of the object phase derivative information acquisition using self-imaging can be distinguished: 1) recording the Fresnel diffraction field intensity pattern [19,21,22] and 2) recording moiré fringes formed just behind the second cross-type grating [17–20, 23–29]. The deciding factor is the detecting camera spatial resolution.

Moiré (grating) interferometry commonly used for in-plane displacement measurement and strain analysis serves as the most important example of conjugate beam 2D grating interferometry [4, 30–33].

Crossed fringe patterns obtained in the works cited above, as single acquisitions, are processed and analyzed using multiple (phase-shifting) [17,19,23,24,27,28,30,31] and single frame methods, e.g., Fourier transform (FT) [18,19,21,22,25,29–31], spatial carrier phase shifting (SCPS) [19,32], regularized phase-tracking [20], continuous wavelet transform (CWT) [33], and spatial mapping using a local cross-correlation between small regions of the specimen-free and specimen-present images [26] (applicable to the cases when the grating frequency is not well separable from the sample spectrum).

The accuracy and calculation speed of single fringe pattern analysis techniques depends mainly on the algorithmic solution applied. Note that a single fringe pattern usually contains noise, nonuniform background and intensity modulation variations. They greatly influence the phase retrieval accuracy, especially in the case of cross-type interference patterns. The problem is quite challenging since in most applications of two dimensional grating interferometry the orthogonal fringe families are mutually multiplied [17–20,22–29] rather than summed [19,21,30–33].

In this paper we present a novel solution to the grid-type fringe pattern analysis called Hilbert-Huang grating interferometry (HHGI) data processing. It combines empirical mode decomposition based single fringe family extraction (and further enhancement) technique with the Hilbert spiral transform (HS) phase demodulation. For the purpose of two orthogonal fringe family separation we introduce the Orthogonal Empirical Mode Decomposition (OEMD) method. It extracts single fringe family using one dimensional order-statistic and averaging filters for mean envelope estimation. The separation is conducted performing single-iteration sifting process to calculate single bidimensional intrinsic mode function (BIMF) by means very similar to the fast and adaptive bidimensional empirical mode decomposition (FABEMD [34],) approach but with the filter size modification. As the separated fringe set is the first BIMF it has zero-mean-value and can be successfully processed (phase demodulated) using the Hilbert transform. The extracted fringe pattern has perpendicular direction to the orientation of filter arrays employed in sifting process. In case of a low quality crossed pattern additional directional filtering is applied to the extracted set using Matlab 'motion' filter with the same orientation as the obtained pattern. It increases the outcome of the extracted pattern further processing (enhancement and normalization) using the ASR-EFEMD (Automatic Selective Reconstruction – Enhanced Fast Empirical Mode Decomposition) technique, recently reported in [34]. Here we propose a simple modification of the automatic selective reconstruction technique called the mutual information detrending (MID). MID automatically defines the number of BIMFs constituting the signal underlying trend and eliminates them from the reconstruction process (as they do not contain valuable fringe information). It is based on the mutual information calculation between adjacent BIMFs. The MID reduces computation time and boosts the ASR method efficiency.

Having our crossed pattern separated using the novel OEMD method and both fringe sets enhanced applying the modified ASR-EFEMD scheme we proceed with calculating their phase distributions taking advantage of the analytic signal approach using the Hilbert spiral transform for quadrature component computation. The HS can be employed efficiently due to the zero-mean-value of each extracted fringe set - it is an inherent feature of every EMD-based algorithm exhibited naturally by OEMD and ASR-EFEMD. The proposed HHGI technique is fast, local, adaptive, automatic and can be efficiently applied to very low quality crossed fringe patterns degraded with noise, low modulation and strong background illumination variations. Although the EMD method is naturally predestined for decomposition of additive superimpositions we will present numerical and experimental studies for both types of squared pattern realizations, i.e., sum and product of mutually orthogonal transmission linear gratings having both binary and cosinusoidal profiles.

The paper is organized as follows. In Section 2 we introduce and describe in detail the HHGI method with the emphasis put on OEMD and MID techniques as main novelties and practical features of the paper. Section 3 contains comprehensive numerical studies corroborating the validity of the proposed approach for crossed fringe pattern processing and analysis. We show the versatility of the HHGI technique analyzing crossed patterns formed by both additive and multiplicative superimpositions of cosinusoidal and binary one-dimensional patterns. Experimental verification of the proposed HHGI scheme is presented in Section 4. We efficiently demodulate real patterns obtained using grating (moiré) interferometry for in-plane displacement measurement (additive patterns with cosine fringe profile) and Talbot interferometry for the lens spherical aberration determination (multiplicative superimposition of binary gratings). Section 5 concludes the paper and points out directions of further studies.

## 2. Method and algorithm description

The input of our algorithm is a single crossed interferogram containing two additively superimposed fringe families with orthogonal phase-modulated carriers. This case is simpler than the product relationship between the component patterns but enables clear explanation of our processing approach without limiting, however, the method generality.

Our goal is a fast, adaptive, automatic and robust separation of the two fringe sets and their efficient phase demodulation. Undertaken task is performed employing a new decomposition approach called Orthogonal Empirical Mode Decomposition (OEMD) and recently reported automatic selective reconstruction using enhanced fast empirical mode decomposition (ASR-EFEMD) technique [34] for further filtering and normalization. Phase demodulation is conducted by means of the Hilbert spiral transform approach calculating the analytic signal. HT can be successfully applied thanks to previous EMD filtering employed for the bias term removal. In our approach we take full advantage of the Hilbert-Huang transformation but not in a classical way (analyzing Hilbert spectrum) but rather in a manner especially tailored for fringe pattern processing and analysis. In this spirit we have named the proposed technique the Hilbert-Huang grating interferometry algorithm and abbreviated it as HHGI.

Our novel OEMD algorithm for the fringe set separation is based on the empirical mode decomposition concept. The one-dimensional version of EMD was firstly developed by Huang and associates [35]. It is adaptive and data-driven approach. Unlike in the Fourier or wavelet methods, no predefined decomposition basis is used; it is rather derived from the signal itself. EMD extracts intrinsic mode functions (IMFs) from signal in the so-called sifting process. It is tailored to process efficiently nonlinear and nonstationary data. The bidimensional EMD (BEMD) is a fully two-dimensional method, which interpolates envelopes of the corresponding data extrema with functions such as bidimensional cubic splines [36,37] or radial based functions [38]. Some comparison of different interpolation techniques is given in [39]. Reported applications of BEMD in the fringe pattern processing and analysis deal with noise reduction in digital speckle interferometry [40], fringe pattern normalization [41], phase measurement in temporal speckle interferometry [42], evaluation of amplitude encoded fringe patterns [43] and, most recently, background and noise removal from fringe patterns [44–47]. The practical impact of BEMD is limited by the calculation time – spline interpolation on the irregular grid is the most expensive part of the algorithm. To overcome this limitation FABEMD (Fast Adaptive BEMD) approach was proposed [48,49] in which the envelope determination of conventional BEMD is modified by replacing the 2D surface interpolation by nonlinear order-statistics-based filtering followed by smoothing operation. Beside significantly shortening the computation time, more accurate estimation of the bidimensional intrinsic mode functions (BIMFs), representing image features at various spatial scales is obtained in many cases. Our group has successfully employed the FABEMD method for studying the moiré phenomenon [50,51], low quality fringe pattern adaptive enhancement [52] and multispectral image fusion [53]. The FABEMD method was also recently introduced to assist face recognition [54]. For further efficient computation time reduction EFEMD (enhanced fast EMD) scheme was reported in [34]. Combined with automatic selective reconstruction (ASR) it constitutes a very efficient fringe pattern adaptive filtering approach called the ASR-EFEMD [34].

In this work we modify the FABEMD algorithm to fulfill the fringe family separation condition. By now different fringe sets present in a single image could be recovered using the FABEMD approach under the assumption that their spatial frequencies differ considerably. High frequency part would be the extracted in one of the first BIMFs whereas the low frequency term would end up in one of the last BIMFs. The reason is that in the FABEMD/EFEMD method the squared/disk-shaped order-statistic and smoothing (averaging) filter masks are used for BIMFs extraction. If the filter window width is denoted as *w* the filter mask size is *[w, w]*. We propose a simple modification of the filter mask size determination to make our decomposition sensitive to single fringe set orientation. Here we apply one dimensional arrays, e.g. *[1,w]* and *[w,1]* in the FABEMD scheme for the orthogonal fringe family extraction. Horizontal filter enables vertical fringe set extraction whereas vertical filter is used for horizontal pattern separation. Our approach is not limited to 0 and 90 degree fringe orientation - it can be applied for different angles (i. e. + 45 and −45 degree) under assumption they are known a priori. This knowledge helps with appropriate “directional” filter mask design representing a novelty and practical feature of the proposed orthogonal empirical mode decomposition technique (OEMD).

In the OEMD algorithm the filter window width *w* is estimated using the approach introduced in EFEMD [34]. We assume that the image extrema are distributed uniformly (which is a natural assumption in cross grid-type interferograms). We detect them, count their number *n* and set

*a*is a parameter controlling the orthogonal decomposition. Its influence on the fringe family separation efficiency will be studied in detail in the next section. Proper selection of

*a*value is especially important in case of noisy patterns. Empirically derived appropriate range for

*a*values will be given. Moreover directional averaging using Matlab's so-called 'motion' filter will be introduced to boost separation efficiency when processing noisy interferograms or multiplicatively superimposed ones.

The proposed OEMD algorithm can be summarized in the following steps.

- 1. Obtaining the orientation angle for both fringe sets (a priori knowledge is needed; the case of grating (moiré) interferometry is simple – the orientation of the fringe families likely corresponds to the specimen grating lines and with the grid we often have 0 and 90 degree orientation angles).
- 2. Detecting the fringe pattern extrema, estimating the filter window width
*w*using Eq. (1) and setting proper*a*value. Designing the orthogonal 1D filter arrays. - 4. Calculating the mean envelope as an arithmetic mean of upper and lower envelopes.
- 5. Subtracting the mean envelope from the initial cross interferogram to obtain single fringe family extraction. Only one iteration is required for single fringe family separation (in a single BIMF). The extracted fringe set orientation is perpendicular to the one dimensional filter array direction. In case of noisy interferograms additional Matlab's 'motion' filter is employed to directionally smooth the separated fringe set. The 'motion' filter has the same direction as the extracted fringe set and the same window width
*w*as estimated before.

Proposed approach can be successfully applied to both additive and multiplicative orthogonal fringe sets superimpositions. The HHGI enables analysis of both binary and cosinusoidal carrier patterns. Although the extraction is straightforward for the additive superimposition, the multiplicative one requires more subtle approach with additional pattern background intensity term removal. For this purpose we develop a novel automatic detrending scheme based on calculating mutual information (MI) [55,56] between adjacent BIMFs to determine the so-called “flag”. This is the first BIMF of the residual term - once establishing the flag mode number we can neglect all the BIMFs with larger number therefore eliminating spurious background intensity term. Detailed computational scheme for both multiplicative and additive superimposition as well as the MI detrending (MID) technique will be described in the next section.

In case of real, noisy grid-like interferograms each fringe set, after extraction, is rather of low quality due to the aforementioned noise and low amplitude modulation (one fringe family is a single extracted mode therefore its contrast is low). To ensure the fringe pattern efficient phase demodulation using the HS we perform automatic adaptive interferogram enhancement and normalization employing the ASR-EFEMD technique [34]. Here we focus on a modification introduced to the automatic selective reconstruction scheme for computation time reduction and algorithm efficiency boost purpose. In [34] we automatically neglected spurious terms: first noisy mode and last decomposition component - the residue. In many cases low spatial frequency signal underlying trend (which we want to remove from our pattern) is “hiding” not only in the monotonic residue but in few last BIMFs [57–60]. This is because of the complexity of analyzed signal and mode mixing phenomenon [61]. Several techniques have been proposed to define which BIMFs (or IMFs in case of 1D signal processing methods) form the signal underlying trend [34,35,44,46,52,57–60]. In our work we use novel MID technique to decide how many last BIMFs should be neglected. Obtained BIMFs number reduction increases efficiency of the ASR-EFEMD and shortens its computation time.

To obtain phase distribution of the filtered fringe pattern (treated as a final outcome of the HHGI) we apply the well-established analytic signal approach. In our case the analytic signal is a complex signal with real part in the form of the filtered fringe pattern and the imaginary part equal to pi/2 phase-shifted fringe quadrature component obtained using the Hilbert transform. For quadrature component calculation we use Hilbert spiral phase transform (HS) introduced in [62, 63]. Its superiority over other Hilbert transform 2D extensions has been proven in [43]. Analytic signal was previously employed by our group for fringe analysis in [34,43,52]. In those papers we used the analytic signal for amplitude demodulation (calculating the modulus of the complex analytic signal). In this paper we calculate fringe phase distribution as an angle of the complex analytic signal (arc tangent Im/Re).

## 3. Numerical studies

This section containing numerical experiments is divided into four subsections. First we define optimal OEMD parameters for efficient separation of additive type crossed interferograms. Next we introduce the MID technique and test our algorithm in the case of multiplicative type grid patterns. Numerical part of the paper ends with investigating the applicability of the proposed algorithm to both cosinusoidal and binary carrier patterns.

#### 3.1 Additive type crossed interferogram processing

### 3.1.1 Noise free additive crossed patterns

We start with a simple additive superposition of noise and background free orthogonal patterns. Our goal is to define the appropriate range for selecting the value of parameter *a.* As mentioned in the Introduction this parameter is used to control the window width *w* (which depends on the detected extrema number, Eq. (1)) of orthogonal filters employed for extracting the individual fringe families. Simulated patterns are shown in Figs. 1(a) and 1(b) alongside with their additive superimposition, Fig. 1(c). One simulated fringe family has its phase modulated with the Matlab 'peaks' function whereas the second one represents a signal modulated with quadratic phase term. In this way both fringe shape and density variations are present in the analyzed superimposition.

For quantitative evaluation of the extraction efficiency we use the image quality index (Q) introduced in [64] and widely used ever since, e.g. in [34,40,41,47,52]. Due to the paper length limitation we will just shortly introduce Q as an index with maximum value of 1 corresponding to the perfect match between images (in our case we compare extracted family and the simulated 'ground truth'). For simplicity we assume that similarity between images drops together with decreasing Q value - corresponding to the algorithm efficiency drop.

In order to define the appropriate range for the parameter *a* selection we performed extraction of both fringe families using *a* values from 0.5 to 10 with the interval of 0.5. Extracted single BIMF (individual fringe family) was normalized using the HS and quantitatively compared with the simulated pattern using Q index. The plot showing extraction quality results is presented in Fig. 2(a). Blue line corresponds to the horizontal fringe family separation and red one relates to the vertical pattern. It is worthy to note that the larger the *a* coefficient value the stronger the filtering applied in OEMD. As the extraction is conducted upon subtracting the mean envelope from initial signal we observe quality saturation for horizontal fringes (blue line). It is a consequence of efficient mean envelope calculation (simple straight vertical pattern) not influenced by the filter window width. The opposite behavior occurs for the vertical fringe family extraction (red line) as we compute mean envelope in the form of horizontal fringes (twisted with 'peaks' function). In this case applying bigger filters corrupts extracted family because of inefficient mean envelope calculation, therefore we observe considerable Q drop - red line in Fig. 2(a). In general the best fit for the noise free cross interferogram separation is *a* = 1. It means that the algorithm sets the filter window width correctly basing only on the detected extrema number. Moreover the proposed OEMD approach is robust to the fringe pattern carrier frequency. As shown in Fig. 2(b) the extraction quality for both fringe families stays above very high 0.965 threshold for wide range of carrier spatial frequencies (normalized with respect to the carrier spatial frequency present in Fig. (1)). This statement is valid for every case studied throughout the paper. Due to very high Q values we do not enclose extracted images for simple noise-free case (they are very similar to Figs. 1(a) and 1(b)).

### 3.1.2 Noisy additive crossed patterns

Our next step involves studying noisy crossed patterns. Due to the presence of noise our OEMD algorithm will detect much more extrema, therefore setting a filter window width much smaller than in the noise free case (for a = 1). The OEMD coefficient *a* adjustment helps tuning the right filter window width in the case of strong noise level (increasing *w* applying a>1). Additionally we propose “directional averaging” performed on the extracted family using the Matlab's 'motion' filter oriented like the separated fringe set (window width equal to the one applied for BIMF calculation, Eq. (1)). The efficiency of the proposed scheme is proved applying *a* values in the range 1 to 20 with the interval of 0.5 for crossed interferogram spoiled with white Gaussian noise of variance 0.02 and extracting both fringe families with and without additional averaging. In case of noisy input pattern the separated families undergone enhancement and normalization using the automatic ASR-EFEMD technique [34]. Plot in Fig. 3(a) shows the results ('ver' denotes vertical fringes, 'hor' denotes horizontal fringes and +/− corresponds to the presence/absence of additional averaging).

The case without additional averaging is similar to the noise free pattern processing - we observe quality saturation for the horizontal pattern extraction and strong quality drop for vertical fringes separation (black and blue lines, respectively). It is interesting to note that additional filtering improves the OEMD efficiency and provides the best results for *a* value between 5 and 10. In this case we also observe the extraction quality saturation but now it occurs for vertical fringes (red line). Straight vertical pattern extraction quality simply increases with the size of the vertical averaging filter applied to the separated individual fringe set (more efficient denoising). Horizontal fringes behave in the opposite way - strong horizontal filtering corrupts the fringe shape (modulated by the 'peaks' function). In consequence there is a window of proper coefficient *a* value between 5 and 10. This is the most appropriate solution for the real noisy patterns processing. It is worth to notice that we obtained rather wide range for *a* parameter selection ensuring flexibility of the method. In general the *a* value needs to be high enough to clear out the noise and low enough to preserve the fringe shape.

The OEMD robustness to noise was tested by processing cross interferograms spoiled with white Gaussian noise with variances from 0.002 to 0.04. Obtained results are shown in Fig. 3(b). Please note that additional directional averaging improves the extraction quality and “slows down” the Q drop rate associated with the increasing noise variance. We set a = 7 and obtained satisfactory results for both fringe families separation.

Exemplary results of HHGI processing are presented in Fig. 4 - it can be seen as a flow chart of the proposed method. Let us consider the simulated interferogram, Fig. 1(c), spoiled white Gaussian noise (variance 0.02). We extracted both fringe families without and with additional directional averaging. Additional filtering reduces noise and enhances fringe shape. Due to low quality of extracted patterns they were enhanced and normalized (only results for horizontal pattern are presented for the sake of keeping paper length reasonable) with the ASR-EFEMD technique [34] (coefficients for ASR-EFEMD: a = 2, additional modulation averaging with *avg_param* = 0.3, see [34] for details). Index Q for the vertical pattern increased from 0.8055 to 0.9464 after additional directional filtering; for horizontal set we obtained the increase from 0.8887 to 0.9537. Important modification of the ASR-EFEMD scheme called mutual information detrending is employed and described in details in the next subsection - it constitutes one of the main novelties proposed in the paper (alongside with the OEMD technique and the HHGI method).

Once an individual fringe family is separated and enhanced we apply Hilbert transform for the quadrature component calculation based on the well-established analytic signal approach; modulated phase can be calculated as the arc tangent of HS(pattern)/pattern [35, 57, 58]. An exemplary wrapped phase map obtained applying HS to the fringe set, separated using OEMD and enhanced using ASR-EFEMD is presented in the middle row of the flow chart, Fig. 4. Unwrapped phase distribution with carrier frequency removed as well as phase error map are shown.

It should be noted that the novel OEMD technique is a first step of the proposed HHGI scheme for the crossed fringe patterns analysis. Its outcome has the form of two resolved orthogonal fringe sets. Second step consist of applying the ASR-EFEMD filtration to the extracted patterns - their quality needs to be enhanced. Although the ASR-EFEMD technique was proposed by us before, its modification (based on the newly developed Mutual Information Detrending technique) introduced in the next section is an original solution and important feature of this paper. Third and final step of the HHGI is the Hilbert spiral transform processing for phase demodulation.

HHGI quantitative phase demodulation evaluation is presented in Fig. 5. Using HS we obtained phase distribution from four extracted families - two vertical fringe patterns (before and after directional averaging), and two horizontal sets. To assess the algorithm performance we subtracted each unwrapped phase distribution from the ideal phase obtaining error maps for vertical fringes ver-, Fig. 5(a), and ver + , Fig. 5(b) as well as horizontal patterns hor-, Fig. 5(c), and hor + , Fig. 5(d). Root-mean-squared phase error calculated for quantitative comparison purposes is equal to 0.1856, 0.0733, 0.1285 and 0.0606 for Figs. 5(a)–5(d), respectively. RMS values prove the importance of the additional directional averaging and in general support the HHGI technique in comparison with other well-established methods employed for the crossed interferogram processing: Fourier transform, Continues Wavelet Transform and Spatial Carrier Phase Shifting.

For example, the Fourier transform processing yields the following results in terms of the image quality index: Q = 0.9154 and Q = 0.9472 for vertical and horizontal fringes, respectively, and phase distribution RMS are equal to 0.1211 and 0.0585, respectively. In comparison with our approach (using additional directional averaging) the numbers for horizontal fringes are very similar whereas the outcome for vertical pattern promotes our local approach over global Fourier processing [1–3,18,19,31,32]. Additional feature of FT is rather sophisticated spectrum filtering required for efficient demodulation (even for single-fringe-family image), especially in case of real crossed interferograms, with noise, background variation, modulation changes and wide range of fringe spatial frequency [18,19,31,32]. Moreover it is hard to automate - the spatial filtering strongly depends on a pattern processed. It complicates greatly when the multiplicative superimposition of two orthogonal families is encountered or binary gratings are studied [19,20].

The continuous wavelet transform [33] gives excellent results, i.e., Q values around 0.9550 and 0.9570 for vertical and horizontal pattern, and RMS values of 0.0913 and 0.0539 for phase distributions. CWT does not share the global character with Fourier transform, it is local and very accurate. Despite local processing the phase distribution of vertical fringes still exhibits RMS higher than that produced with our approach. The CWT method in its state-of-the-art [33] is very efficient and accurate method but the calculations are rather complicated (several parameters need to be defined). It is extremely time consuming to perform full separation of two fringe sets (with high accuracy).

The SCPS method is characterized by fast and rather simple computations but has stringent requirements on introducing appropriate carrier spatial frequency and processing only the additive superimposition between two orthogonal linear fringe families [31,32].

The proposed OEMD technique is very fast and straightforward - it requires three simple statistical and averaging filtering operations (one statistical and averaging filtering for single BIMF calculation and one directional averaging for BIMF filtering to boost algorithm efficiency). OEMD is local, adaptive and robust to noise, uneven background illumination, modulation variations and fringe spatial frequency changes. Additional unique feature of the proposed OEMD (and HHGI) scheme is its direct applicability to the case of multiplicative superimposition of fringe families including binary (in general noncosinusoidal) profile of superimposed components. These statements will be proven below in Section 3.3. It is important to note that our approach is very fast - full separation and enhancement for most complicated patterns (noisy, low amplitude modulation, uneven background intensity, nonlinear profile) lasted 6 seconds at most (512x512 patterns, medium class PC).

#### 3.2 Fringe pattern detrending using mutual information

Automatic selective reconstruction using the enhanced fast empirical mode decomposition [34] plays crucial role in the proposed Hilbert-Huang Grating Interferometry data processing. Once we have extracted individual fringe family from crossed interferogram its quality is rather poor, especially in case of severely defected fringe pattern. The separated fringe set needs to be enhanced and normalized – this task is realized employing ASR-EFEMD [34]. The main idea behind the ASR-EFEMD technique is taking into account only valuable information from each BIMF, i.e. sharply extracted fringe regions with high amplitude modulation values. Inherent feature of every EMD-based algorithm is a narrow spatial frequency range of sharply extracted fringes in a single BIMF - this spatial frequency decreases over the set of BIMFs. ASR analyzes each pixel of each BIMF separately to established the best fit for the reconstruction. It automatically neglects last decomposition component called the residue and the first BIMF in case of noisy pattern [34]. Due to the mode mixing phenomenon (or general complexity of a signal under study) low spatial frequency spurious component (background intensity term to be neglected) is often stored not only in the residue but in last several BIMFs. In this paper we present a novel approach called Mutual Information Detrending (MID) for automatic determination of so called “flag” - the last valuable information carrying BIMF number. Having defined the “flag” we can neglect every BIMF with higher number as they constitute the image underlying trend.

The principle of MID is a calculation of mutual information (MI) between adjacent BIMFs. MI can be an indicator of the correlation between two variables [55,56]. Detailed explanation of the MI processing is reported in [55,56] and references therein. What differs MI from general autocorrelation calculation is sensitivity to nonlinear correlations, robustness to different size of the partition elements and rather simple and straightforward implementation. MI can be considered superior to the general joint entropy in measuring similarity between variables [55,56] because it takes into account joint data information. In [56] authors use MI for efficient intrinsic mode functions selection and formulation for enhanced bearing condition monitoring.

In our case compared variables are two consecutive BIMF. In Fig. 6(a) simulated fringe pattern with uneven background intensity distribution is presented. Figures 6(b)–6(j) show first eight BIMFs and the residue obtained after decomposition using the EFEMD algorithm (in this study we use filter size selection like in FABEMD OSFW 2 algorithm [49,50] but with fast and efficient implementation of EFEMD [34]). We discard the residue and compute MI between two adjacent BIMFs starting with BIMF1 and BIMF2 calculating f(x) = MI[BIMF(x),BIMF(x + 1)] for x = 1,2...7. The plot illustrating our results is depicted in Fig. 7(a). We can see that MI between adjacent BIMFs monotonically drops and has its minimum in x = 4. It means that BIMF4 is the last mode similar in terms of MI to the previous valuable information (fringe) carrying BIMFs. After BIMF4 we compared BIMF5 and BIMF6 and have found that the “similarity” increases. This is a simple way to determine in which BIMF (called “flag”) the information carrying fringe component ends and background intensity component starts. To make it more visible and robust we computed the function g(x) in the form of g(x) = f(x + 1)/f(x) obtaining the plot in Fig. 7(b). Now strong peak is an indicator of the “flag”.

We summed up the last four BIMFs and the residue forming the fringe pattern background, Fig. 6(k). Although some fringes can be discerned this fact does not affect the fringe pattern successful reconstruction using first four BIMFs for automatic selective reconstruction, Fig. 6(l). The obtained background term is not perfect but it was never the goal of MID. Using mutual information we can automatically narrow down the “ASR suspects” to first four BIMFs. In this case the MID method reduced calculation time by a factor of 2. The purpose of MID is to neglect more BIMFs than just the residue. It is purely empirical approach and in case of fringe patterns, which can be divided into modulated cosine term and the background term, is an efficient way to determine which part of the set of BIMFs corresponds to the cosine term and which constitutes the background component.

Similar approach was previously employed in [44,46] where authors used general autocorrelation for determining how many first BIMFs form the high frequency noise component. Here we use the mutual information (differences from the autocorrelation approach were listed in the previous paragraph) and decide how many last BIMFs constitute the underlying trend. We detect the noise by means of noise-to-signal ratio (NSR) calculation for the first BIMF, as proposed in [34]. If NSR>1 we neglect first few BIMFs from the MID 'flag' determination scheme (up to 3-4) as they clearly store high frequency components (due to the detected noise presence). Please note that for the automatic selective reconstruction we still neglect only BIMF1 (when NSR>1).

#### 3.3 Multiplicative type crossed interferogram processing

Having proven the validity of the proposed approach for additive type crossed interferogram processing, Section 3.1, we proceed with studying multiplicatively superimposed structures. This kind of superimposition is simpler to realize, therefore more popular in experimental practice [17–20,22–29]. Efficient processing of this type of crossed patterns is more challenging, however. Often well-established tools (e.g. FT, SCPS) fail to demodulate multiplicatively formed crossed fringes, therefore narrowing their usability only to additive patterns with clearly separated spectrum [20, 31, 32].

The HHGI approach can be directly applied to the multiplicative type crossed pattern separation and further processing in the form described for additive type patterns. Exemplary noise-free multiplicative type fringe pattern is presented in Fig. 8(a). Vertical fringes separated with a = 1 are shown in Fig. 8(b). Multiplicative superimposition yields low spatial frequency background term corresponding to the peaks function modulating phase of horizontal fringes. This term is responsible for low quality (Q = 0.7643) normalization outcome, Fig. 8(c). Background intensity component is successfully determined using the MID technique, Fig. 8(d). Neglecting low spatial frequency underlying trend we obtained almost ideal vertical pattern extraction with Q = 0.9842, Fig. 8(e). Using MID for horizontal fringes was unnecessary - straightforward OEMD method application gave excellent outcome (Q = 0.9579).

To corroborate the algorithm robustness to noise we spoiled simulated multiplicative type crossed pattern with white Gaussian noise of variance 0.02, Fig. 9(a), and performed successful separation of both fringe families, Figs. 9(b) and 9(c), obtaining satisfactory Q index values, i.e., Q = 0.9548 and 0.9579, respectively. Again, the choice *a* = 7 was found to be very effective in the case of noisy pattern separation. Calculations for other variance values exhibit the same character as those for the additive type superimposition, thus they are not repeated here. In previous section we showed that such a high Q value corresponds to very low RMS of HS-demodulated phase distribution - we do not repeat phase calculations to keep the paper length reasonable.

#### 3.4 Binary versus cosinusoidal profile

In experimental practice we often encounter fringe pattern nonlinearities in the form of, e.g., nonlinear fringe projection and/or detection. Moreover we might simply use binary amplitude gratings instead of cosinusoidal ones. Noncosinusoidal profile of examined structures often complicates the calculations causing the need for additional filtering in Fourier/Hilbert transform based methods [65,66]. In order to corroborate the validity of the HHGI data processing for binary crossed patterns we simulated additive and multiplicative superimposition type binary patterns. Additive superimposition results are presented in Fig. 10. It is worthy to note that separated vertical and horizontal sets, Figs. 10(b) and 10(d) exhibit binary profile. For efficient Hilbert transform based phase demodulation it is crucial to obtain cosinusoidal fringe patterns.

This task is realized decomposing binary patterns and reconstructing them using the ASR-EFEMD technique neglecting first BIMFs containing high spatial frequency information (responsible for binary profile). Filtered cosinusoidal vertical and horizontal fringes, Figs. 10(c) and 10(e) achieved Q index values around 0.95. Phase demodulation yielded RMS values 0.1172 and 0.0536, respectively, corroborating efficiency of ASR-EFEMD binary pattern filtering. Very similar results were obtained for multiplicative type superimpositions of binary orthogonal structures, Fig. 11. Again the MID technique was used to filter out the spurious background term from vertical set. We report Q index values around 0.95 for both fringe families and phase RMS values 0.1212 and 0.0653 for vertical and horizontal terms, respectively. Presented quantitative analysis corroborate the HHGI method ability to efficiently process binary crossed patterns.

## 4. Experimental data processing

In this Section we present experimental verification of the proposed HHGI technique. First we process low quality crossed interferogram obtained in the grating (moiré) interferometry setup. Extracted vertical and horizontal fringe patterns encode information about in-plane displacements fields *u(x,y)* and *v(x,y)* in the x and y directions, respectively. This is the case where two orthogonal cosinusoidal patterns are superimposed additively [33]. Second real crossed pattern to be processed was obtained illuminating the binary amplitude crossed grating (mutliplicative superimposition) and recording its distorted self-image in the Fresnel diffraction region. This self-image encodes information about first derivative (in two orthogonal directions) of the wave front with spherical aberration introduced by a tested lens. We do not use a second real grating as an analyzer (placed at the self-image distance) - we record the self-image, extract two orthogonal fringe families and then numerically design the digital reference pattern to be superimposed (by multiplication) with each set to create moiré fringes. Our approach can be seen as digital Talbot interferometry.

Efficient analysis of two aforementioned experimental fringe patterns comprehensively corroborates the versatility of the proposed HHGI technique. It is well suited for processing of both additive and multiplicative superimpositions of two orthogonal, low quality structures (noisy, uneven modulation and considerable background illumination variation). Both binary and cosinusoidal component structures can be successfully analyzed.

#### 4.1 Grating (moiré) interferometry data processing

In Fig. 12(a) low quality crossed interferogram obtained in the grating (moiré) interferometry setup is presented. Two fringe patterns with encoded in-plane displacement information *u(x,y)* and *v(x,y)* were additively superimposed in the stretched fabric loading experiment [33] (for the description of moiré interferometry simultaneous pattern recording see [30–33, 67] and references therein).

First we extract the vertical fringe set, Fig. 12(b), using the OEMD method with additional filtering (a = 10). Employing the modified ASR-EFEMD scheme we obtain enhanced, filtered and normalized vertical pattern, Fig. 12(c). It is worthy to note that high noise level, very low modulation and fringe spacing variations present in the initial crossed pattern were successfully overcame by the proposed OEMD algorithm aided with the modified ASR-EFEMD scheme. Global approaches (e.g., FT) produce errors connected with spatial frequency variations (see Section 3). Wrapped phase distribution containing in-plane displacement field *u(x,y)* was obtained using the Hilbert spiral transform, Fig. 12(d). The horizontal fringe set extraction is shown in Figs. 12(e)–12(g).

After unwrapping and plane-fitting for the carrier phase removal we obtained continuous *u(x,y)* and *v(x,y)* in-plane displacement fields (from vertical and horizontal fringe sets, respectively). They are presented in Figs. 13(a) and 13(d), respectively. For comparison purposes we implemented 5-frame temporal phase shifting algorithm (TPS [1–3],) obtaining unwrapped phase maps corresponding to the in-plane displacement fields *u(x,y)* and *v(x,y).* They are outlined in Figs. 13(b) and 13(e), respectively. High noise level and very low modulation present in each recorded frame resulted in quite noisy phase maps calculated using the TPS method. Nevertheless results obtained using 5 frames correspond very well to those produced with the proposed single-frame HHGI approach. Additionally, robustness to noise, low modulation and fringe spatial frequency can be confirmed comparing cross-sections in Figs. 13(c) and 13(f). We can observe that HHGI provides efficient filtering combined with high accuracy.

The HHGI processing results can be compared with the ones obtained using the recently proposed [33] continuous wavelet transform approach, Fig. 14. It can be clearly noted that CWT outcome contains fringe discontinuities as a consequence of wrong selection of coefficients dictated by severe fringe defects and spatial frequency variation.

#### 4.2 Talbot interferometry data processing

In Fig. 15(a) we show the self-image of the binary crossed grating (multiplicative superimposition of linear gratings). We can observe that the self-image is distorted by spherical aberration present in the illuminating beam. Obtaining the phase distribution would aid collimating lens testing. The spherical wave front aberration was generated by turning the collimator objective (focal length of 200 mm) by 180 degrees about the vertical axis. In this way the objective was set opposite its aberration correction direction [68,69].

We start with the vertical and horizontal fringe set extraction, Figs. 15(b) and 15(c), respectively. Due to considerable background illumination variations we can observe low contrast regions corresponding to the dark regions of the initial image, Fig. 15(a). Our approach (taking advantage of the Hilbert transform normalization) is robust to the uneven background illumination, however. Normalized patterns are shown in Figs. 15(d) and 15(e).

In Talbot interferometry the second grating is placed in the self-image plane producing moiré fringes with low spatial frequency. Here we propose digital realization of this solution for phase extraction from vertical and horizontal fringes. We numerically simulate sinusoidal detecting gratings with fringe spacing close to the spacing of extracted terms and parallel lines, Figs. 15(f) and 15(g). We digitally superimpose the simulated detecting gratings with the extracted fringe set obtaining vertical and horizontal moiré fringes, Figs. 16(a)–16(d). For each fringe set we create two moiré patterns with simulated reference pattern phase shifted by π/2.

In this way, having filtered moiré fringes using the EFEMD decomposition and the Hilbert transform normalization (moiré fringe pattern processing with FABEMD algorithm was described in detail in [51]), Figs. 16(e)–16(h), we can take advantage of the analytic signal and compute wrapped moiré fringe phase maps, Figs. 17(a) and 17(b). We simply define the wrapped phase as an angle of a complex signal with its real and imaginary parts in the form of two π/2 phase shifted patterns (quadrature components). After unwrapping and the linear phase term removal we obtain characteristic phase distribution corresponding to orthogonal spatial derivatives of the spherical wave front under test, Figs. 17(c) and 17(d). This experimental work corroborates the proposed HHGI approach applicability to binary structures processing, its robustness to considerable background illumination variations and ability to support digital Talbot interferometry.

## 5. Conclusions

In the paper we proposed the Hilbert-Huang grating interferometry data processing technique for single-frame crossed fringe pattern efficient filtering and phase analysis. It consists of three main computational routines: (1) Orthogonal Empirical Mode Decomposition employed for the orthogonal fringe set separation, (2) ASR-EFEMD enriched by the Mutual Information Detrending procedure applied for further single-family pattern fast and adaptive filtering and (3) Hilbert spiral transform for enhanced fringe phase demodulation. The OEMD algorithm and the MID technique are considered main novelties and practical features. Presented comprehensive numerical studies corroborate the validity, robustness and versatility of the proposed HHGI scheme. It can be successfully applied to crossed patterns obtained using both multiplicative and additive type superimpositions of binary and cosinusoidal unidirectional component structures. Experimental verification by means of efficient processing and analysis of crossed patterns recorded using grating (moiré) interferometry and Talbot interferometry techniques has been provided. HHGI is fast, adaptive, automatic, local and enables successful analysis of complex low quality fringes due to very efficient filtering provided by the modified ASR-EFEMD technique.

We would like to point out that our cross-grating interferometry data processing method obviates the need to use additional coherent optical filtering system [66]. Developing novel algorithmic solutions for fringe pattern analysis makes measurement system hardware simplifications possible. We believe that the solution presented in the paper will find interest among researchers working in different spectrum regions, especially in visible and X-ray ones. Extended analysis should be performed for X-ray applications, however. Separated paper will be devoted to this topic.

## Acknowledgments

This work was supported, in part, by Ministry of Science and Higher Education budget funds for science in the years 2012-2015 as a research project in the “Diamond Grant” programme, partly by National Science Center, Poland, under the project DEC-2012/07/B/ST7/01392, in part by statutory funds and grant founded by the Dean of Faculty of Mechatronics, Warsaw University of Technology. Work of MT was partially supported by the European Union in the framework of European Social Fund through the Warsaw University of Technology Development Programme.

## References and links

**1. **D. W. Robinson and G. Reid, *Interferogram Analysis: Digital Fringe Pattern Measurement* (Institute of Physics Publishing, 1993).

**2. **D. Malacara, M. Servin, and Z. Malacara, *Interferogram Analysis for Optical Testing* (Marcel Dekker, 1998).

**3. **D. Malacara, ed., *Optical Shop Testing* (John Wiley, 2007).

**4. **D. Post, B. Han, and P. Ifju, *High Sensitivity Moirè* (Springer, 1994).

**5. **K. Patorski, *Handbook of the Moirè Fringe Technique* (Elsevier, 1993).

**6. **S. Yokozeki and T. Suzuki, “Shearing interferometer using the grating as the beam splitter,” Appl. Opt. **10**(7), 1575–1580 (1971). [CrossRef] [PubMed]

**7. **S. Yokozeki and T. Suzuki, “Shearing interferometer using the grating as the beam splitter. Part 2,” Appl. Opt. **10**(7), 1690–1693 (1971). [CrossRef] [PubMed]

**8. **A. W. Lohmann and D. E. Silva, “An interferometer based on the Talbot effect,” Opt. Commun. **2**(9), 413–415 (1971). [CrossRef]

**9. **A. W. Lohmann and D. E. Silva, “A Talbot interferometer with circular gratings,” Opt. Commun. **4**(5), 326–328 (1972). [CrossRef]

**10. **D. E. Silva, “Talbot interferometer for radial and lateral derivatives,” Appl. Opt. **11**(11), 2613–2624 (1972). [CrossRef] [PubMed]

**11. **O. Kafri, “Noncoherent method for mapping phase objects,” Opt. Lett. **5**(12), 555–557 (1980). [CrossRef] [PubMed]

**12. **K. Patorski, “Diffraction effects in moiré deflectometry: comment,” J. Opt. Soc. Am. A **3**(5), 667–668 (1986). [CrossRef]

**13. **E. Keren and O. Kafri, “Diffraction effects in moiré deflectometry: reply to comment,” J. Opt. Soc. Am. A **3**(5), 669–670 (1986). [CrossRef]

**14. **E. Bar-Ziv, S. Sgulim, O. Kafri, and E. Keren, “Temperature mapping in flames by moire deflectometry,” Appl. Opt. **22**(5), 698–705 (1983). [CrossRef] [PubMed]

**15. **G. Oster, M. Wasserman, and C. Zwerling, “Theoretical interpretation of moiré patterns,” J. Opt. Soc. Am. **54**(2), 169–175 (1964). [CrossRef]

**16. **K. Patorski, “The self-imaging phenomenon and its applications,” in *Progress in Optics*, E. Wolf ed., vol. 27, 1–108 (Elsevier, 1989).

**17. **H. Canabal, J. A. Quiroga, and E. Bernabeu, “Improved phase-shifting method for automatic processing of moiré deflectograms,” Appl. Opt. **37**(26), 6227–6233 (1998). [CrossRef] [PubMed]

**18. **J. A. Quiroga, D. Crespo, and E. Bernabeu, “Fourier transform method for automatic processing of moiré deflectograms,” Opt. Eng. **38**(6), 974–982 (1999). [CrossRef]

**19. **H. Canabal and E. Bernabeu, “Phase extraction methods for analysis of crossed fringe patterns,” Proc. SPIE **3744**, 231–240 (1999). [CrossRef]

**20. **J. Villa, J. A. Quiroga, and M. Servín, “Improved regularized phase-tracking technique for the processing of squared-grating deflectograms,” Appl. Opt. **39**(4), 502–508 (2000). [CrossRef] [PubMed]

**21. **J. L. Flores, B. Bravo-Medina, and J. A. Ferrari, “One-frame two-dimensional deflectometry for phase retrieval by addition of orthogonal fringe patterns,” Appl. Opt. **52**(26), 6537–6542 (2013). [CrossRef] [PubMed]

**22. **H. H. Wen, E. E. Bennett, R. Kopace, A. F. Stein, and V. Pai, “Single-shot x-ray differential phase-contrast and diffraction imaging using two-dimensional transmission gratings,” Opt. Lett. **35**(12), 1932–1934 (2010). [CrossRef] [PubMed]

**23. **I. Zanette, T. Weitkamp, T. Donath, S. Rutishauser, and C. David, “Two-dimensional x-ray grating interferometer,” Phys. Rev. Lett. **105**(24), 248102 (2010). [CrossRef] [PubMed]

**24. **I. Zanette, C. David, S. Rutishauser, T. Weitkamp, M. Denecke, and C. T. Walker, “2D grating simulation for X-ray phase-contrast and dark-field imaging with a Talbot interferometer,” AIP Conf. Proc. **1221**, 73–79 (2010). [CrossRef]

**25. **H. Itoh, K. Nagai, G. Sato, K. Yamaguchi, T. Nakamura, T. Kondoh, C. Ouchi, T. Teshima, Y. Setomoto, and T. Den, “Two-dimensional grating-based X-ray phase-contrast imaging using Fourier transform phase retrieval,” Opt. Express **19**(4), 3339–3346 (2011). [CrossRef] [PubMed]

**26. **K. S. Morgan, D. M. Paganin, and K. K. Siu, “Quantitative single-exposure x-ray phase contrast imaging using a single attenuation grid,” Opt. Express **19**(20), 19781–19789 (2011). [CrossRef] [PubMed]

**27. **S. Rutishauser, I. Zanette, T. Weitkamp, T. Donath, and C. David, “At-wavelength characterization of refractive x-ray lenses using a two-dimensional grating interferometer,” Appl. Phys. Lett. **99**(22), 221104 (2011). [CrossRef]

**28. **S. Berujon, H. Wang, I. Pape, K. Sawhney, S. Rutishauser, and C. David, “X-ray submicrometer phase contrast imaging with a Fresnel zone plate and a two dimensional grating interferometer,” Opt. Lett. **37**(10), 1622–1624 (2012). [CrossRef] [PubMed]

**29. **H. Wang, S. Berujon, I. Pape, S. Rutishauser, C. David, and K. Sawhney, “At-wavelength metrology using the moiré fringe analysis method based on a two-dimensional grating interferometer,” Nucl. Instrum. Methods Phys. Res. A **710**, 78–81 (2013). [CrossRef]

**30. **L. Salbut, “Multichannel system for automatic analysis of *u,v,w* displacements in grating interferometry,” in *Physical Research*, W. Juptner and W. Osten, eds. 19, 282–287 (Akademie, 1993).

**31. **M. Kujawinska and M. Pirga, “Fringe pattern analysis for the investigation of dynamic processes in experimental mechanics,” in *Physical Research*, W. Juptner and W. Osten, eds. 19, 391–394 (Akademie, 1993).

**32. **M. Pirga and M. Kujawinska, “Two directional spatial-carrier phase-shifting method for analysis of crossed and closed fringe patterns,” Opt. Eng. **34**(8), 2459–2466 (1995). [CrossRef]

**33. **K. Pokorski and K. Patorski, “Separation of complex fringe patterns using two-dimensional continuous wavelet transform,” Appl. Opt. **51**(35), 8433–8439 (2012). [CrossRef] [PubMed]

**34. **M. Trusiak, M. Wielgus, and K. Patorski, “Advanced processing of optical fringe patterns by automated selective reconstruction and enhanced fast empirical mode decomposition,” Opt. Lasers Eng. **52**, 230–240 (2014). [CrossRef]

**35. **N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for non-linear and non-stationary time series analysis,” Proc. R. Soc. Lond. A **454**(1971), 903–995 (1998). [CrossRef]

**36. **C. Damerval, S. Meignen, and V. Perrier, “A fast algorithm for bidimensional EMD,” IEEE Signal Process. Lett. **12**(10), 701–704 (2005). [CrossRef]

**37. **C. B. Barber, D. P. Dobkin, and H. Huhdanpaa, “The quickhull algorithm for convex hulls,” ACM Trans. Math. Softw. **22**(4), 469–483 (1996). [CrossRef]

**38. **J. C. Nunes, Y. Bouaoune, E. Delechelle, O. Niang, and P. Bunel, “Image analysis by bidimensional empirical mode decomposition,” Image Vis. Comput. **21**(12), 1019–1026 (2003). [CrossRef]

**39. **S. M. A. Bhuiyan, N. O. Attoh-Okine, K. E. Barner, A. Y. Ayenu-Prah, and R. R. Adhami, “Bidimensional empirical mode decomposition using various interpolation techniques,” Adv. Adapt. Data Anal. **1**(2), 309–338 (2009). [CrossRef]

**40. **M. B. Bernini, A. Federico, and G. H. Kaufmann, “Noise reduction in digital speckle pattern interferometry using bidimensional empirical mode decomposition,” Appl. Opt. **47**(14), 2592–2598 (2008). [CrossRef] [PubMed]

**41. **M. B. Bernini, A. Federico, and G. H. Kaufmann, “Normalization of fringe patterns using the bidimensional empirical mode decomposition and the Hilbert transform,” Appl. Opt. **48**(36), 6862–6869 (2009). [CrossRef] [PubMed]

**42. **M. B. Bernini, A. Federico, and G. H. Kaufmann, “Phase measurement in temporal speckle pattern interferometry signals presenting low-modulated regions by means of the bidimensional empirical mode decomposition,” Appl. Opt. **50**(5), 641–647 (2011). [CrossRef] [PubMed]

**43. **M. Wielgus and K. Patorski, “Evaluation of amplitude encoded fringe patterns using the bidimensional empirical mode decomposition and the 2D Hilbert transform generalizations,” Appl. Opt. **50**(28), 5513–5523 (2011). [CrossRef] [PubMed]

**44. **Y. Zhou and H. Li, “Adaptive noise reduction method for DSPI fringes based on bi-dimensional ensemble empirical mode decomposition,” Opt. Express **19**(19), 18207–18215 (2011). [CrossRef] [PubMed]

**45. **X. Zhou, T. Yang, H. Zou, and H. Zhao, “Multivariate empirical mode decomposition approach for adaptive denoising of fringe patterns,” Opt. Lett. **37**(11), 1904–1906 (2012). [CrossRef] [PubMed]

**46. **X. Zhou, H. Zhao, and T. Jiang, “Adaptive analysis of optical fringe patterns using ensemble empirical mode decomposition algorithm,” Opt. Lett. **34**(13), 2033–2035 (2009). [CrossRef] [PubMed]

**47. **X. Zhou, A. G. Podoleanu, Z. Yang, T. Yang, and H. Zhao, “Morphological operation-based bi-dimensional empirical mode decomposition for automatic background removal of fringe patterns,” Opt. Express **20**(22), 24247–24262 (2012). [CrossRef] [PubMed]

**48. **S. M. A. Bhuiyan, R. R. Adhami, and J. F. Khan, “A novel approach of fast and adaptive bidimensional empirical mode decomposition,” in *Proc. IEEE Int. Conf. Acoustic, Speech and Signal Process.* (Institute of Electrical and Electronics Engineers, 2008), 1313–1316. [CrossRef]

**49. **S. M. A. Bhuiyan, R. R. Adhami, and J. F. Khan, “Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation,” EURASIP J. Adv. Signal Process. **2008**(21), 728356 (2008). [CrossRef]

**50. **K. Patorski, K. Pokorski, and M. Trusiak, “Fourier domain interpretation of real and pseudo-moiré phenomena,” Opt. Express **19**(27), 26065–26078 (2011). [CrossRef] [PubMed]

**51. **M. Trusiak and K. Patorski, “Space domain interpetation of incoherent moiré superimpositions using FABEMD,” Proc. SPIE 8697, 18th Czech-Polish-Slovak Optical Conference on Wave and Quantum Aspects of Contemporary Optics, 869704 (December 18, 2012). [CrossRef]

**52. **M. Trusiak, K. Patorski, and M. Wielgus, “Adaptive enhancement of optical fringe patterns by selective reconstruction using FABEMD algorithm and Hilbert spiral transform,” Opt. Express **20**(21), 23463–23479 (2012). [CrossRef] [PubMed]

**53. **M. Wielgus, M. Bartys, A. Antoniewicz, and B. Putz, “Fast and adaptive bidimensional empirical mode decomposition for the real-time video fusion, ” in Proc. 15th Int. Conf. Inform. Fusion (FUSION*,*2012), 649–654.

**54. **Y. Zhou, S.-T. Zhou, Z.-Y. Zhong, and H.-G. Li, “A de-illumination scheme for face recognition based on fast decomposition and detail feature fusion,” Opt. Express **21**(9), 11294–11308 (2013). [CrossRef] [PubMed]

**55. **C. Studholme, D. L. G. Hill, and D. J. Hawkes, “An overlap invariant entropy measure of 3D medical image alignment,” Pattern Recognit. **32**(1), 71–86 (1999). [CrossRef]

**56. **S. Osman and W. Wang, “An enhanced Hilbert-Huang transform technique for bearing condition monitoring,” Meas. Sci. Technol. **24**(8), 085004 (2013). [CrossRef]

**57. **Z. Wu, N. E. Huang, S. R. Long, and C. K. Peng, “On the trend, detrending, and variability of nonlinear and nonstationary time series,” Proc. Natl. Acad. Sci. U.S.A. **104**(38), 14889–14894 (2007). [CrossRef] [PubMed]

**58. **Z. Yang, B. W.-K. Ling, and C. Bingham, “Trend extraction based on separations of consecutive empirical mode decomposition components in Hilbert marginal spectrum,” Measurement **46**(8), 2481–2491 (2013). [CrossRef]

**59. **P. Flandrin, P. Goncalves, and G. Rilling, “Detrending and denoising with empirical mode decomposition,” In EUSIPCO 2004. September 6–10, Vienna, Austria (2004)

**60. **G. Rilling, P. Flandrin, and P. Goncalves, “On empirical mode decomposition and its algorithms,” IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03, Grado (I), 2003.

**61. **G. Rilling and P. Flandrin, “One or two frequencies? The empirical mode decomposition answers,” IEEE Trans. Signal Process. **56**(1), 85–95 (2008). [CrossRef]

**62. **K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform,” J. Opt. Soc. Am. A **18**(8), 1862–1870 (2001). [CrossRef] [PubMed]

**63. **K. G. Larkin, D. J. Bone, and M. A. Oldﬁeld, “Natural demodulation of two-dimensional fringe patterns. II. Stationary phase analysis of the spiral phase quadrature transform,” J. Opt. Soc. Am. A **18**(8), 1871–1881 (2001). [CrossRef] [PubMed]

**64. **Z. Wang and A. C. Bovik, “A universal image quality index,” IEEE Signal Process. Lett. **9**(3), 81–84 (2002). [CrossRef]

**65. **L. Xiong and S. Jia, “Phase-error analysis and elimination for nonsinusoidal waveforms in Hilbert transform digital-fringe projection profilometry,” Opt. Lett. **34**(15), 2363–2365 (2009). [CrossRef] [PubMed]

**66. **N. Sun, Y. Song, J. Wang, Z.-H. Li, and A.-Z. He, “Volume moiré tomography based on double cross gratings for real three-dimensional flow field diagnosis,” Appl. Opt. **51**(34), 8081–8089 (2012). [CrossRef] [PubMed]

**67. **J. Schmit, K. Patorski, and K. Creath, “Simultaneous registration of in- and out-of-plane displacements in modified grating interferometry,” Opt. Eng. **36**(8), 2240–2248 (1997). [CrossRef]

**68. **K. Patorski, “Grating shearing interferometer with variable shear and fringe orientation,” Appl. Opt. **25**(22), 4192–4198 (1986). [CrossRef] [PubMed]

**69. **K. Patorski, “Conjugate lateral shear interferometry and its implementation,” J. Opt. Soc. Am. A **3**(11), 1862–1870 (1986). [CrossRef]