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:

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

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

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.

 

Fig. 1 Simulated vertical (a) and horizontal (b) fringe patterns and their additive superimposition (c).

Download Full Size | PDF

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)).

 

Fig. 2 Left - plot representing the relationship between coefficient a and Q index values for extracted patterns (blue line - horizontal, red line - vertical term); right - plot representing the relationship between normalized carrier frequency of both fringe families and Q index for extracted patterns (blue line - horizontal, red line - vertical term).

Download Full Size | PDF

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).

 

Fig. 3 On the left - plot representing relations between coefficient a and Q index values for extracted patterns in the case of a noisy interferogram; on the right - plot representing relations between added white Gaussian noise variance and obtained Q index values for extracted patterns, see text for explanations.

Download Full Size | PDF

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).

 

Fig. 4 Proposed algorithm diagram. Additive superimposition of noisy fringe sets is considered. The OEMD method is used for resolving two orthogonal fringe families, the ASR-EFEMD technique is applied for extracted single fringe pattern filtering and normalization, and the Hilbert spiral transform is employed for phase demodulation.

Download Full Size | PDF

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.

 

Fig. 5 Phase error maps obtained in the case of vertical (a-b) and horizontal (c-d) fringe set extraction from noisy crossed pattern, Fig. 4(a), using regular OEMD method (a,c) and additional directional averaging (b,d). Additional directional filtering decreases the error level in the final phase distribution.

Download Full Size | PDF

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”.

 

Fig. 6 (a) Simulated simple fringe pattern with uneven background illumination, (b)-(j) first eight BIMFs and the residue obtained using the EFEMD method, (k) background intensity term calculated using the MID approach, (l) fringe pattern normalized using the ASR-EFEMD method aided by MID.

Download Full Size | PDF

 

Fig. 7 Two plots illustrating the functions (a) f(x) = MI(BIMF_x, BIMF_{x + 1}) and (b) g(x) = f(x + 1)/f(x) calculated using MID for BIMF set depicted in Fig. 6. We established the flag number equal to 4.

Download Full Size | PDF

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).

 

Fig. 8 (a) Simulated multiplicative superimposition type crossed pattern, (b) vertical fringes extracted using OEMD with additional directional filtering, note strong uneven background term, (c) pattern (b) after normalization, (d) background intensity term extracted from (b) using MID, (e) vertical fringes enhanced and normalized with ASR-EFEMD aided by MID (with uneven background term removed), (f) horizontal fringe set extracted from (a).

Download Full Size | PDF

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.

 

Fig. 9 (a) Multiplicative type crossed pattern spoiled with white Gaussian noise of variance 0.02, (b) vertical and (c) horizontal fringe families extracted using OEMD with additional directional filtering (a = 7).

Download Full Size | PDF

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.

 

Fig. 10 (a) Simulated additive type binary crossed pattern, (b) vertical binary fringes extracted using the OEMD method, (c) cosinusoidal vertical fringes obtained filtering out high frequency BIMFs, (d) horizontal binary fringes extracted using OEMD, (e) cosinusoidal horizontal fringes obtained filtering out high frequency BIMFs, (f) unwrapped phase map of (e) - the final outcome of the HHGI algorithm.

Download Full Size | PDF

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.

 

Fig. 11 (a) Simulated multiplicative type binary crossed pattern, (b) vertical binary fringes extracted using the OEMD method, (c) cosinusoidal vertical fringes obtained filtering out high frequency BIMFs and the background term with MID, (d) horizontal binary fringes extracted using OEMD, (e) cosinusoidal horizontal fringes obtained filtering out high frequency BIMFs, (f) unwrapped phase map of (e) - the final outcome of the HHGI algorithm.

Download Full Size | PDF

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).

 

Fig. 12 (a) Experimental crossed interferogram obtained in the grating (moiré) interferometry setup (310x410 pixels image corresponding to the 4,5x6 mm² area of the sample under test; sample with transferred reflective diffraction grating 1200 lines/mm; carrier fringes introduced by tilting illumination beams [33]), (b) BIMF containing vertical fringes extracted from (a) using the OEMD algorithm (a = 10), (c) vertical fringes enhanced and normalized using modified ASR-EFEMD, (d) wrapped phase fringes obtained using Hilbert transform, (e-g) analogous processing results for the horizontal fringe set.

Download Full Size | PDF

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.

 

Fig. 13 Unwrapped phase maps corresponding to the u(x,y) displacement field obtained using the proposed HHGI approach (a) and the 5-frame temporal phase shift method (b); plot illustrating cross-sections along 50th row of the u(x,y) displacement field obtained using HHGI (red line) and TPS (blue line) methods (c). Unwrapped phase maps corresponding to the v(x,y) displacement field obtained using HHGI scheme (d) and TPS (e); (f) plot illustrating the cross-section along 200th row of the v(x,y) displacement field obtained using HHGI (red line) and TPS (blue line).

Download Full Size | PDF

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.

 

Fig. 14 Continuous Wavelet Transform processing results: extracted horizontal fringe set (a) and its unwrapped phase distribution with carrier term removed (b); extracted vertical fringe set (c) and its unwrapped phase distribution with carrier term removed (d).

Download Full Size | PDF

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].

 

Fig. 15 (a) Recorded self-image (750x750 pixels) of crossed binary grating (two 8 lines/mm binary gratings with orthogonal lines superimposed multiplicatively); extracted (b) vertical and (c) horizontal fringe families using OEMD - observe contrast variations corresponding to the uneven background illumination distribution of (a); (d)-(e) two fringe families normalized (note that the proposed method is robust to strong background illumination variations and applicable to binary structures); (f-g) two sinusoidal detecting gratings (reference patterns) with fringe spacing designed to match the spatial frequency of information carrying orthogonal terms (d-e). Numerically simulated reference gratings will be used for generating digital moiré fringes.

Download Full Size | PDF

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.

 

Fig. 16 Top: two moiré pattern pairs with pi/2 phase shift generated by multiplicative superimposition of the extracted term and simulated detecting (reference) gratings, Figs. 15(d)–15(g); moiré phase shifting was obtained by numerical error-free phase shifting of the reference grating. Bottom: moiré fringes filtered using the EFEMD decomposition and Hilbert transform normalization.

Download Full Size | PDF

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.

 

Fig. 17 Wrapped moiré phase maps obtained using well-established analytic signal tool: (a) calculated from Figs. 16(e) and 16(f) and (b) calculated from Figs. 16(g) and 16(h); unwrapped phase maps obtained from (c) vertical and (d) horizontal moiré fringes - note characteristic distributions of the orthogonal spatial derivatives of the spherical wave front under test.

Download Full Size | PDF

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.