Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns

Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng

Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns Chenxing Wang a,b, Qian Kemao c,n, Feipeng Da a,d a

School of Automation, Southeast University, Nanjing 210096, China Multi-Platform Game Innovation Centre, Nanyang Technological University, 639798 Singapore, Singapore c School of Computer Engineering, Nanyang Technological University, 639798 Singapore, Singapore d Key Laboratory of Measurement and Control of CSE, Ministry of Education, Southeast University, Nanjing 210096, China b

art ic l e i nf o

a b s t r a c t

Article history: Received 6 January 2016 Received in revised form 12 April 2016 Accepted 19 April 2016

Fringe patterns are often produced from optical metrology. It is important yet challenging to reduce noise and remove a complicated background in a fringe pattern, for which empirical mode decomposition based methods have been proven useful. However, the mode-mixing problem and the difficulty in automatic mode classification limit the application of these methods. In this paper, a newly developed method named regenerated phase-shifted sinusoids assisted empirical mode decomposition is introduced to decompose a fringe pattern, and subsequently, a new noise-signal-background classification strategy is proposed. The former avoids the mode-mixing problem appearing during the decomposition, while the latter adaptively classifies the decomposition results to remove the noise and background. The proposed method is testified by both simulation and real experiments, which shows effective and robust for fringe pattern analysis under different noise, fringe modulation, and defects. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Image analysis Phase retrieval Fringe analysis

1. Introduction A fringe pattern from optical metrology is usually expressed as

I (x, y) = A (x, y) + B (x, y) cos [ωx + ϕ (x, y)] + n (x, y),

(1)

where A(x, y) and B(x, y) are the background and amplitude intensity, respectively; ω denotes the spatial frequency; ϕ(x, y) is the phase distribution; and n(x, y) is random noise. Normally the main purpose of fringe analysis is to retrieve ϕ(x, y). It is thus necessary to remove the irrelevant terms A(x, y) and n(x, y) in order to enhance the core part of the fringe pattern, B(x, y)cos[ωx þ ϕ(x, y)], which is called a phase-modulated (PM) signal. In the Fourier transform profilometry (FTP), when the phase is simple, the PM signal is band limited in the Fourier domain and can be easily separated from other parts if a high carrier frequency is provided [1,2]. However, when the phase is more complicated, a large error will appear because the spectrum of the PM signal is mixed with the spectra of other parts. In addition to this limitation, the Heisenberg's uncertainty principle also limits almost all the spectral analysis methods such as the windowed FTP [3,4], the ridge of wavelet transform (WT) [5] or the S-transform [6,7].

Empirical mode decomposition (EMD) based methods provide another approach for fringe pattern analysis. EMD is a data-driven technique that aims to decompose a non-stationary signal into a series of mono components (named as intrinsic mode functions, IMFs) [8]. It has been used to suppress noise [9] and eliminate the background to suppress the zero spectrum [10,11] in FTP. Furthermore, EMD combined with Hilbert transform (HT) is increasingly used for phase extraction in recent years, where EMD extracts the PM signal and HT constructs an analytic signal of the PM signal for phase retrieval [12–15]. A signal can be decomposed by EMD as follows,

I (x) =

∑

IMFk (x) + r (x) ,

(2)

1≤ k ≤ K

where k is the index of an IMF; K is the total number of IMFs; r(x) is the residue. As the IMFs range from high frequency to low frequency, they are classified into three groups as follows, I (x ) =

∑

IMFk (x ) +

k < k1

⎡ IMFk (x ) + ⎢ ∑ IMFk (x ) + r (x ) ⎢⎣ k1≤ k ≤ k 2 k 2
∑

⎤ ⎥ , ⎥⎦

(3)

where k1 and k2 are two critical indexes dividing all the IMFs into three groups; the first group ∑k < k IMFk (x ) represents the noise, the 1

n

Corresponding author. E-mail address: [email protected] (Q. Kemao).

IMFk (x ) is the desired PM signal and the last group means the background. The classification result directly second group ∑k

1≤ k ≤ k2

http://dx.doi.org/10.1016/j.optlaseng.2016.04.018 0143-8166/& 2016 Elsevier Ltd. All rights reserved.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

2

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

influences the accuracy of the phase retrieval. There are two technical challenges in getting Eq. (3). First, in EMD results, the socalled mode-mixing problem (MMP) often occurs, which gives an IMF uncertain physical meaning. Second, even without the MMP, it is difficult to correctly classify the three groups, i.e., k1 and k2 of Eq. (3) are not easy to be determined. For a fringe pattern, the MMP is usually caused by uneven distribution of noise, including either high-frequency noise or lowfrequency noise. Additionally, if the MMP exists in the high-frequency IMFs, it will spread to low-frequency IMFs. In [13], soft thresholding is used to eliminate the mixed PM signal existing in IMF1. Considering the possible MMP in any IMFs, in [14,15], local parts of each bidimensional IMF (BIMF) are picked out according to local amplitudes and then combined into the final result through a weighted mean. The methods show good results regardless of solving the MMP, but they still encounter challenges for various applications if the noise or background has the similar amplitude with the PM signals. Ensemble EMD (EEMD) [16] has been thought as the mainstream technique to solve the MMP. It processes an ensemble of the white-noise added signals and takes the average of all results as the end result. EEMD has been expanded to bidimensional EEMD [17] and multivariate EEMD [18,19], and further improved into Complete EEMD (CEEMD) to reduce the residue noise or spurious modes [20,21] in these years. All ensemble-based methods suffer from long computing time to process a big ensemble of data. In order to improve the efficiency, we have conducted some studies of adding a designed “noise” to achieve the results of EEMD [23,23]. However, these methods focus on the MMP in the high-frequency IMFs caused by noise, but ignore the MMP in other scales. As for the classification of noise and background, k1 ¼2 and k2 ¼ K are often set, postulating that only IMF1 is the noise and only the r(x) is the background [12,24]. To cope with more complicated α situations, in [25], the standard deviation of ∑k = 1 IMFk (x ) is computed as s(α), and k1 ¼ α  1 is set when s(α) becomes dramatically larger than s(α 1), where a threshold is needed. In [17], the power of autocorrelation parameters for each BIMF is computed and then the abrupt change of the computed values is found to set k1. As for background removal, k2 is often chosen manually [10,26]. In [27], the marginal entropies of BIMFs are computed to obtain mutual information, in order to evaluate the correlation between IMFs, and then the mutual information is used to determine k2 [28]. In [13], the mean of an IMF is tested whether the IMF belongs to the background through a threshold, providing that an IMF belonging to the background has a larger mean. Unfortunately, the MMP of the low-frequency IMFs always leads to unpredictable cases, thus, only the correlation between IMFs or the characteristic analysis of an IMF itself cannot work robustly to correctly determine k2. In [29], the method realizes k1 ¼ k2 ¼ K¼1 cleverly and presents a good result for background reconstruction based on the condition that the fringe pattern is clean. In another work [30], the time average subtraction method is used to eliminate the background, but this is unavailable for a single frame of fringe pattern analysis. We also developed some criteria to set k2 based on the overall changes of frequencies and amplitudes between IMFs [22,23], but more robust measures are needed in practice. In this paper, a novel regenerated phase-shifted sinusoids assisted EMD (RPSEMD) method is introduced to solve the MMP of EMD, and then a new strategy is proposed for mode classification to reduce noise and remove the background of a fringe pattern. The RPSEMD generates different sinusoids adaptively in different stages of decomposition to solve the MMP in all IMFs. The new noise-signal-background classification strategy determines k1 and k2 correctly for the resulted mono-component IMFs. Both the RPSEMD method and the classification strategy are automatic. The

proposed method is robust to cope with a variety of complex situations for fringe pattern analysis in optical metrology. Experiments also show high accuracy and efficiency of the proposed method.

2. The RPSEMD algorithm used for analyzing fringe patterns 2.1. The mode-mixing problem EMD is a data-driven method that aims to iteratively decompose a signal into a series of mono-component IMFs. We call each iteration, namely, the process producing an IMF, as a stage of decomposition. So IMFk(t) in Eq. (2) is the product after the kth stage of decomposition. Each stage of EMD starts from detecting the extremum points of the tested signal [8]. Based on the detected extrema, IMFk(t) is determined. We define an mono component of a signal as an intrinsic mode (IM) represented as IM(t)¼ a(t)cos[2πf(t)], where a(t) and f(t) mean the instantaneous amplitude and instantaneous frequency respectively [31]. Then the obtained result in the kth J stage can be written as IMFk (t ) = ∑j = 1 aj (t ) cos [2πf j (t ) t ], where J is the number of IMs. If the detected extrema belong to a single IM, namely, J¼ 1, then IMFk(t) is an ideal result. On the contrary, if the detected extrema belong to multiple IMs, namely, J 41, the decomposed IMF contains more than one IM, and the MMP occurs. 2.2. Solutions to the MMP and the RPSEMD algorithm EEMD [16] is a powerful method for solving the MMP, which can adjust the extrema of each scale in a signal by adding a large series of white noise because the white noise has scales populated throughout the time-frequency space. However, the ensemble size of white noise is required to be very large to make sure the added white noise is finally canceled out by average. Inspired by EEMD, we have proposed to add a designed sinusoid as “a noise” to the signal of a fringe pattern [22,23]. The methods successfully solve the MMP in high-frequency IMFs but they ignore the possible MMP in low-frequency IMFs. Moreover, in [23], we design the amplitude and frequency of the added “noise” only through a simple average method. Thus an advanced technique is needed to robustly cope with signals with high complexity. As a response to the above requirement, we developed the RPSEMD algorithm recently [32]. Within each stage of decomposition, a sinusoid is designed and added to make the extrema uniformly distributed in order to avoid the MMP. The procedure of RPSEMD is summarized as follows: 1. Apply EMD to a signal I(t) to get initial IMFs and design a sinusoid s(t| a, f, θp) where a, f, and θp are amplitude, frequency and an initial phase of the sinusoid, respectively. The amplitude and frequency are determined according to the initial IMFs, while the initial phase is 0; 2. Apply EMD to I(t) þs(t| a, f, θp) only to get a temporary IMF1(t); 3. Repeat Step 2 with θp increased by 2π/P for another (P  1) P times and get the final first IMF as c1 (t ) = ⎡⎣∑p = 1 IMF1 (t ) ⎤⎦ /P ; 4. Remove c1(t) from I(t), i.e., I(t)’I(t)-c1(t), then repeat Steps 1-3 on the residue I(t) to get all other IMFs. The final residue is denoted as r(t). The result of RPSEMD can also be represented by Eq. (2), but this time, the MMP is largely resolved. Details of each step can be

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

found in [32] where the MATLAB code is also provided. The RPSEMD is selected for fringe analysis in this paper because of its following advantages: First, it solves the MMP in each stage of decomposition, which is important for fringe analysis where the MMP often appears due to: (a) inevitable noise, including high-frequency noise and lowfrequency noise; (b) the complicity of a severely deformed fringe with uneven extrema; and (c) shadows and a complicated background. Furthermore, we introduce the clustering analysis to adaptively analyze the signal and design the added sinusoid according to a targeted IM. This way makes the designed signal effective to solve the MMP with a high probability, thus the RPSEMD algorithm is able to cope with a variety of complex fringe patterns in practice. Second, the phase shift of the sinusoid is useful for retaining the details of each separated IM. This is advantageous to improve the accuracy and expand the measurement range for optical metrology. Third, the RPSEMD provides an active means to solve the MMP, which improves the efficiency greatly compared with the EEMD. Typically, the RPSEMD can be tens or even hundreds of times faster than the EEMD because EEMD needs a large ensemble size. This is important when processing a 2D fringe pattern in a line-byline manner.

3. The strategies of classifying the results of RPSEMD for denoising and background removal As the next key step for fringe analysis, the mono-component IMFs successfully generated by the RPSEMD need to be further classified into a desired PM signal, and undesired noise and background. In other words, the critical indexes k1 and k2 in Eq. (3) should be determined. Correspondingly, two strategies will be proposed and elaborated in the next two subsections. 3.1. Determining the critical index k1 for denoising To determine k1 is to distinguish the signal IMFs from the noise IMFs. A simple way to characterize each IMF is to interrogate its autocorrelation function. The autocorrelation function R(τ) of a signal I(t) is defined as

R (τ ) =

∞

∫−∞ I (t ) I (t − τ ) dt,

(4)

where τ is a time delay. This function has the following properties: (PA1) R(-τ)¼ R(τ) when I(t) is a real signal; (PA2) |R(τ)| rR(0); (PA3) R(τ) has the same period with I(t) if I(t) is periodic; (PA4) if an IMF is full of noise, its R (τ) is random. In [16], the powers of the autocorrelation have been used, but setting qualitative strategies based on these power values is not easy. We propose a new strategy inspired by the properties PA3 and PA4. For each IMF, its R(τ) is calculated first. Since R(τ) is even (property PA1), it is sufficient to just consider non-negative τ, i.e., the right half of R(τ). Next, the first three extrema are detected, and denoted as e(i) (i¼1,2,3). Note that the first extremum is always a maximum located at τ[e(1)] ¼ 0 (property PA2). The second and third extrema are a minimum and maximum, respectively, which are numerically searched. These three extrema form a cycle, which will look like a sinusoid for a signal IMF (property PA3), while less so for a noise IMF (property PA4). As the period of fringe has to be larger than 3 pixels according to the sample limit [33], in our strategy, an IMF is determined as noise if

3

Dl = τ [e (2)] − τ [e (1)] = τ [e (2)] < 3 or Dr = τ [e (3)] − τ [e (2)] < 3,

(5)

where Dl and Dr mean the spans of the left and right halves of the first cycle, respectively. An example is given below. The decomposition of a signal is shown in Fig. 1(a). It is quite obvious that c1(t) is a noise IMF, while the signal is concentrated in c3(t). The IMF c2(t) is a little tricky. As it includes a part of weak signal around x ¼450, it should be classified as a signal IMF to protect the integrity of the signal. The autocorrelation of c1(t) and c2(t), denoted as R1(τ), and R2(τ), are then calculated and shown in Fig. 1(b) and (c), respectively. As expected, R1(τ) is random, while R2(τ) presents more regular oscillations. Next, for clarity, R1(τ) and R2(τ) (0 r τ r50) are shown in Fig. 1(d) and (e), respectively. Furthermore, their first three extrema of each correlation function are detected and highlighted by red stars. It is easy to calculate that for R1(τ), Dl ¼ Dr ¼1, which satisfies the inequalities in (5), while for R2(τ), Dl ¼9 and Dr ¼8, which do not satisfy the strategy. This leads to a desired decision of k1 ¼2. 3.2. Determining the critical index k2 for background removal After denoising, the remaining IMFs are mainly the PM signal and the background, called as PM group and background group respectively for simplicity. The next step is to determine the critical index k2 to separate these two groups. Some reliable and useful properties for the transition between the k2th and (k2 þ1)th IMFs, i.e, the last IMF in the PM group and the first IMF in the background group respectively, are listed below: (PB1) Global frequency: the main frequency of an IMF becomes smaller abruptly, supposing that background has slow spatial variation compared with the PM signal. The main frequency of IMFk(t) is estimated from the autocorrelation function, i.e., fmain (k)¼1/τ[e(3)], referring to (5). (PB2) Global energy: the energy of IMFk2 in the signal group usually tends to be the smallest, while that of IMFk2+ 1 becomes larger contrarily, forming an abrupt increase at the transition. The energy E(k) is also estimated from the autocorrelation function as Rk(0). (PB3) Local frequency and amplitude: the background often has locally flat brightness, which means most of the points in IMFk2+ 1 have non-ignorable amplitudes but have very small frequencies. To respond to such local features, a new parameter is defined as

ξ (k ) = mean ⎡⎣ ak ins (t )/fk ins (t ) ⎤⎦, ins

(6)

ins

where ak (t) and fk (t) are the instantaneous amplitude and instantaneous frequency, respectively. Consequently ξ(k2 þ1) is either the largest around the transition or increases the fastest from ξ(k2). (PB4) Orthogonality: the orthogonality between the two transition IMFs is weaker while the IMFs within each group are roughly orthogonal to each other. The orthogonality of two successive IMFs is estimated as their cross correlation, T 1 O (k ) = T ∑t = 1 ck (t ) ck + 1 (t ), where T is the length of the IMFs [8]. While all of the above properties are useful in identifying the transition IMFs, a strategy by organically combining all these properties is proposed in order to solve the challenging problem of signal-background separation in various practical fringe patterns. As a necessary condition for the transition, PB1 is thus considered to be used first. We scan the IMFs from IMFk1 (k1 has been determined as in Section 3.1) and compares the frequency difference between two successive IMFs. The scanning stops when fmain(kf)/fmain(kf þ1) 4 2 is met, where the twice difference of frequencies is a conservative limit for differentiating distinct

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

0.2 0 -0.2 0.5 0 -0.5 0.5 0 -0.5 0.2 0 -0.2 0.05 0 -0.05 0.05 0 -0.05 0.55 0.5 0.45

100

200 300 x /pixel

4

4

2

2

R2(τ)

R1(τ)

r

c6

c5

c4

c3

c2

c1

4

0 -2 -511

-200

0

200

511

400

500

0 -2 -511

-200

τ /pixel

0 τ /pixel

200

511

3

4

2

R1(τ)

R2(τ)

2 0

1 0 -1

-2

20

40 60 τ /pixel

80

100

-2

20

40 60 τ /pixel

80

100

Fig. 1. (a) The decomposition of a signal; (b) (c): the autocorrelation of c1(t) and c2(t) respectively; (d) (e): the local parts of autocorrelation in (b), (c) correspondingly.

components [34], such as the PM signal and background. We then conservatively reject the IMFs with k rkf for further processings as they are thought to belong to the PM group. This is, however, a rough location of the transition. Next, local features are utilized to more precisely locate the transition. We compute ξ(k) for (k Zkf) according to Eq. (6), where IMFkf is still involved to highlight the changes around IMFkf + 1 although we have rejected IMFkf as the PM signal. To find the most likely transition, three situations are considered successively: if maxima of ξ(k) exist, the first maximum at k ¼km is searched and kξ ¼km  1 is the output of this step; if no maximum exists but minima exist, similarly, the first minimum is searched and kξ ¼km is the output; if neither a maximum nor minimum exists, we compute the numerical difference of ξ(k) and kξ ¼km is the output

if the fastest increase appears from ξ(km) to ξ(km þ1). In our various tests, these two steps work very effectively in identifying k2 by just using kξ as k2. Nevertheless, properties PB2 and PB4 can be added as an optional step to decide k2, because there may appear low-frequency noise IMFs being adjacent to the transition. This optional step is executed within the IMFs with k Z kξ  1 where IMFkξ − 1 is involved also to highlight the changes of E (k) and O(k) around IMFkξ . We detect the minima of E(k) and O(k) (k Zkξ  1) respectively and record the corresponding indexes as two sets. If the interaction of the two sets is non-empty, k2 is determined as the smallest value in the intersection; otherwise, k2 keeps the result as kξ determined earlier. Note that in our strategy the residue r(t) is not involved.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4. Simulations and experiments

Table 1 The error comparison of three methods.

4.1. Effectiveness of the proposed method 4.1.1. Effectiveness of signal decomposition by RPSEMD A row of a fringe pattern, shown in Fig. 2(a), is processed to show the effectiveness of the proposed RPSEMD. Fig. 2(b) and (c) is the results from EEMD [15] and improved CEEMD (ICEEMD) [21], two mainstream methods for solving the MMP of EMD, with their codes downloaded from [35,36], respectively. The results are the best ones from those obtained by changing different parameters. We can see that the orthogonality of c2 and c3 are both weak for these two methods, which may cause low-frequency noise in the following IMFs. The result of our RPSEMD in Fig. 2(d) shows much better orthogonality with the least number of IMFs. 4.1.2. Effectiveness of the proposed classification strategies For Fig. 2(b) and (c), we manually set k1 ¼2 and k2 ¼5, which are set based on several attempts to minimize the errors of the extracted PM signal, by comparing it with the ideal signal. The errors are listed in Table 1. As for Fig. 2(d), we first apply the strategy in (5) to set k1. For each IMF, the autocorrelation is computed, from which e(1), e

5

serr γerr

EEMD

ICEEMD

RPSEMD

0.024 0.083

0.024 0.073

0.019 0.059

a

serr and γerr are the standard deviation of the errors and the largest error, respectively.

(2) and e(3) are detected to obtain both Dl and Dr in (5). As a result, we obtain Dl ¼1, Dr ¼4 for R1(τ) and Dl ¼7, Dr ¼7 for R2(τ), enabling us to set k1 ¼2 according to (5). We turn to determine k2. The main frequencies from c2(t) to c5 (t) are computed. The scanning of the frequency ratio, fmain(k)/fmain (k þ1), starts from k¼2. As fmain(2)/fmain(3) is evaluated to be 2.5, which is already larger than 2, we roughly delimit the background as the group from c3(t) to c5(t) based on PB1, and the scanning stops. Next, we compute ξ(k) (2 r kr 5) defined in Eq. (6), which is shown in Fig. 3(a). As a minimum at k¼ 4 is found, k2 ¼4 is determined. This decision can be further backup by the optional step using energy and orthogonality. We plot the energy E(k) and orthogonality O(k) in Fig. 3(b) and (c), respectively. Both of them

Fig. 2. (a) A signal of fringe pattern; the results of: (b) EEMD; (c) ICEEMD; (d) the proposed RPSEMD.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

3

10

2

ξ(k)

E(k)

15

5 0

(a)

2

3

4 k

5

1

O(k)

6

1 0

(b)

2

3

4 k

5

x 10

-3

0.5

0

(c)

2

3

4

5

k

Fig. 3. The change of: (a) the parameter ξ(k); (b) the power E(k); (c) the orthogonality O(k).

Fig. 4. (a) The error of PM signal reconstructed by: (a) EEMD, (b) ICEEMD, (c) RPSEMD; (d) the PM signal got by RPSEMD.

Fig. 5. Performance comparison of different methods for PM signal extraction. (a) and (b): Fringe patterns with standard deviations of noise at 0.08 and 0.2, respectively; (c) the ideal background; (d) the ideal phase. For the remaining sub-figures (e1)–(j3): row 2 to row 4 show the results of RPSEMD, ICEEMD and MOBEMD, respectively, while columns 1, 2, 3 show the separated noises, separated backgrounds, phase error from (a), respectively, and columns 4, 5, 6 show the separated noises, separated backgrounds, phase error from (b), respectively.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎ Table 2 Error comparison of different methods. Method

RPSEMD ICEEMD MOBEMD

Fig. 4(a) (Q¼ 0.81)

Table 3 Error comparison of different methods. Method

QPM

serr

γerr

0.927 0.918 0.876

0.183 0.187 0.227

0.847 0.915 1.310

7

RPSEMD ICEEMD MOBEMD

Fig. 4(b) (Q¼ 0.41)

Fig. 6(a1)

QPM

serr

γerr

0.817 0.781 0.482

0.246 0.265 3.94

1.63 5.68 44.6

a

serr and γerr are the standard deviation of the errors and the largest error, respectively, and QPM is the Q value.

have a minimum at k¼4, and supports the earlier decision of k2 ¼4. The results of our strategy, k1 ¼2 and k2 ¼4, are consistent with our manual analysis. Then we compute the error of the extracted PM signal by our method and compare it with the errors by EEMD and ICEEMD, as shown in Fig. 4(a)–(c) and Table 1. The extracted PM signal by RPSEMD is also shown in Fig. 4(d), where the results of EEMD and ICEEMD are not shown as they look similar. The RPSEMD shows the best accuracy but costs the least time using MATLAB with the same PC, i.e., it costs 0.66 s for RPSEMD but 31 s for EEMD and 20 s for ICEEMD, where the ensemble size is 200 for both the latter. This improvement of efficiency will be more apparent for processing a fringe pattern using the line-by-line operation. 4.2. The performance of the proposed method for analysis of different fringe patterns 4.2.1. For fringe patterns with different levels of noise We first test the performance of the proposed method for fringe patterns with different levels of noise. According to Eq. (1), two fringe patterns (512  512) are simulated as follows: the phase is ϕ(x, y)¼3  peaks where “peaks” is an MATLAB built-in function; the varying background is A(x, y) ¼2.5  ∂[peaks(y, x)]/ ∂xþ0.5; the amplitude is B(x, y) ¼0.25; the carrier frequency is ω ¼2π/16; and finally the random noises are additive, with standard deviations of 0.08 in Fig. 5(a) and 0.2 in Fig. 5(b). The index Q value expressed in [37] are computed as 0.81 and 0.41 for Fig. 5 (a) and (b), respectively, which means the quality of Fig. 5(b) is very low due to the heavy noise. Fig. 5(c) displays the ideal background and Fig. 5(d) shows the phase ranging from 19.65 rad and 24.31 rad.

RPSEMD ICEEMD MOBEMD

Fig. 6(a2)

serr

γerr

serr

γerr

0.113 0.120 0.125

1.19 1.67 5.60

0.098 0.141 0.178

0.768 2.12 6.50

a

serr and γerr are the standard deviation of errors and the largest error, respectively.

The proposed method processes Fig. 5(a) and (b) row by row. The removed noise and background are shown in Fig. 5(e1) and (f1) for Fig. 5(a), and in Fig. 5(h1) and (i1) for Fig. 5(b), respectively, which prove that our strategies work effectively under different noise environments. After denoising, the Q values are computed to quantify the denoised fringe patterns, as listed in Table 2, where the quality is improved twice when the noise is heavy. The background is smoothed using median filtering, which then will be very close to the ideal one. After denoising and background, the PM signals are obtained, which can be further smoothed by median filtering with a 5  5 window. The smoothed PM signals as not shown to save the space, but their phases by using HT and unwrapping are evaluated. Fig. 5(g1) and (j1) are phase errors for Fig. 5(a) and (b), respectively, from which we also compute the standard deviation and the largest value to quantify these errors, as listed in Table 2. Our method is seen to work satisfactorily even for very noisy fringe patterns. For comparison, ICEEMD is tested again, with standard deviations of added noises as 0.3 and 0.5 for Fig. 5(a) and (b), respectively, and the ensemble size is 200. The decomposed results are classified using our proposed strategies. Furthermore, another method called morphological operation-based bi-dimensional EMD (MOBEMD) [29] is also tested. Its code is downloaded from [38], where 2D digital wavelet transform (DWT) and a db9 wavelet are selected for pre-denoising as in [29]. Our method outperforms the two methods both visually from Fig. 5 and quantitatively from Table 2. 4.2.2. For fringe patterns with different modulations The influence of fringe modulation is examined. As shown in Fig. 6(a1) and (a2), the fringe patterns are simulated with the similar parameters to the ones in Fig. 6, but the amplitude is

Fig. 6. Influence of fringe modulation (a1): a fringe patterns with a varying modulations; (b1): the obtained noise; (c1): the obtained background; (d1) the phase error; (a2): a fringe pattern with a different modulation; (b2-d2): the obtained noise, background and phase error, respectively.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

8

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 7. Influence of a defect (a) a fringe pattern with a defect; (b) the obtained noise; (c) the obtained background; (d) the background after median filtering on (c); (e) the unwrapped phase; (f) a 3D view of (e).

Fig. 8. (a) The fringe pattern; (b) the removed noise; (c) the result by RPSEMD for a signal marked by a red line in Fig. 8(a); (d) the effect of the proposed method for a line; (e) a locally enlarged area of Fig. 8(d); (f) the whole reconstructed background; (g) the background after median filtering; (h) the wrapped phase; (i) the 3D map phase.

changed as B(x, y) ¼0.25  (x þ100)/256 for Fig. 6(a1) and B(x, y)¼ 2  {∂peaks(y, x)/∂xþ 0.15} for Fig. 6(a2) separately. The standard deviation of the noise in fringe patterns is 0.04. Fig. 6(b1)–(b2) show the removed noise and Fig. 6(c1)–(c2) are obtained background of the corresponding fringe patterns. The noise and background are successfully separated. The phase errors are shown in Fig. 6(d1)–(d2). These two fringe patterns are also processed by ICEEMD and MOBEMD, and the phase errors are shown in Table 3. Again, our method performs most satisfactory. 4.2.3. For a fringe pattern with defects We artificially set a black box in Fig. 6(a2) to simulate a defect. The fringe pattern is shown in Fig. 7(a). The removed noise by using our method is shown in Fig. 7(b). The obtained background and its version after median filtering are shown in Fig. 7(c) and (d), respectively. The unwrapped phase and the 3D view of the phase are shown in Fig. 7(e) and (f), respectively. It is interesting to see that the defect does not affect the success of demodulation and does not spread its error to the phase area. We also quantify the errors of the retrieved phase for the whole fringe pattern excluding the defect area in which the phase is meaningless. Then the standard deviation is 0.1005 and the largest error is obtained as 0.8423. We process the same fringe pattern that does not has this defect, namely Fig. 6(a2), and pick the same area excluding the defect area in Fig. 7(a). As a result, the standard deviation of errors is 0.0983 and the largest error is 0.7682. The quantitative results

further show that the influence of the defect is very small to the result of our method. 4.3. A real experiment A deformed fringe pattern obtained by projecting a fringe pattern onto a foam board is shown in Fig. 8(a). This fringe pattern suffers from the shadow, reflect light and uneven illumination. The removed noise is shown in Fig. 8(b). To show the details, we pick a row of the fringe pattern, which is marked by a red line in Fig. 8(a), and process it using the proposed method. Fig. 8(c) is the decomposition result of the proposed RPSEMD, which shows reasonable and clear physical meaning. After judged by the proposed strategies, c1(t) is recognized as the noise, while c4(t), c5(t) as well as the residue r(t) are assigned to the background. In Fig. 8(d), the red line is the original signal, the blue line is the denoised signal and the green line is the reconstructed background. A portion of Fig. 8(d) between the 600th pixel and the 800th pixel is enlarged in Fig. 8(e). The red line and blue line are not very differentiable as the noise is light. However, the green line representing the background is obviously satisfactory. The whole background is shown in Fig. 8(f) and the result after median filtering is displayed in Fig. 8(g), which describes the light change well even for the shadow area. Fig. 8(h) and (i) are the wrapped phase and the unwrapped one in 3D, both of which are seen successful.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i

C. Wang et al. / Optics and Lasers in Engineering ∎ (∎∎∎∎) ∎∎∎–∎∎∎

5. Conclusion The mode-mixing problem (MMP) and the automatic classification of the decomposed results are two challenges in EMD based fringe analysis methods. A novel RPSEMD method is selected to solve the MMP and new strategies are developed to properly classify and remove the noise and background from a fringe pattern. The proposed RPSEMD can result in reasonable and pure IMFs for the subsequent classification. The proposed classification strategies, by using properties of different types of components, are shown to be effective and robust. In this paper, the RPSEMD algorithm is implemented line by line, which avoids the extrema detection problem for saddle points in 2D space. However, extending the proposed method to 2D is of interest and being studied due to the need for analyzing fringe patters with varying fringe directions.

Acknowledgments This work is supported by the National Natural Science Foundation of China (61405034, 51475092), China Postdoctoral Science Foundation (2014M551486), the Post-doctoral Research Funding Program of Jiangsu Province (1302001 A), the international Postdoctoral Exchange Fellowship Program 2014, and the Multi-plAtform Game Innovation Centre (MAGIC) funded by the National Research Foundation, Prime Minister's Office, Singapore under its IDM Futures Funding Initiative.

References [1] Takeda M, Ina H, Kobayashi S. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1981;72:156–60. [2] Su X, Zhang Q. Dynamic 3-D shape measurement method: a review. Opt Laser Eng 2010;48:191–204. [3] Zhong J, Weng J. Generalized Fourier analysis for phase retrieval of fringe pattern. Opt Express 2010;18:26806–20. [4] Gao W, Qian K. Statistical analysis for windowed Fourier ridge algorithm in fringe pattern analysis. Appl Opt 2012;51:328–37. [5] Li S, Su X, Chen W. Wavelet ridge technique in optical fringe pattern analysis. J Opt Soc Am A 2010;27:1245–54. [6] Zhong M, Chen W, Jiang M. Application of S-transform profilometry in eliminating nonlinearity in fringe pattern. Appl Opt 2012;51:577–87. [7] Da F, Dong F. Windowed Fourier transform profilometry based on improved S-transform. Opt Lett 2012;37:3561–3. [8] Huang NE, Shen Z, Long S, Wu M, Shin H, Zheng Q, Yen N, Tung C, Liu H. The empirical mode decomposition and Hilbert spectrum for non-linear and nonstationary time series analysis. Proc R Soc A 1998;454:903–95. [9] Bernini MB, Galizzi GE, Federico A, Kaufman GH. Evaluation of the 1D empirical mode decomposition method to smooth digital speckle pattern interferometry fringes. Opt Laser Eng 2007;45:723–9. [10] Li S, Su X, Chen W, Xiang L. Eliminating the zero spectrum in Fourier transform profilometry using empirical mode decomposition. J Opt Soc Am A 2009;26:1195–201. [11] Wang C, Da F. Phase demodulation using adaptive windowed Fourier

9

transform based on Hilbert–Huang transform. Opt Express 2012;20:18459–77. [12] Rodriguez FA, Federico A, Kaufmann GH. Hilbert transform analysis of a time series of speckle interferograms with a temporal carrier. Appl Opt 2008;47:1310–6. [13] Zhang C, Ren W, Mu T, Fu L, Jia C. Empirical mode decomposition based background removal and de-noising in polarization interference imaging spectrometer. Opt Express 2013;21:2592–605. [14] Trusiak M, Patorski K, Wielgus M. Adaptive enhancement of optical fringe patterns by selective reconstruction using FABEMD algorithm and Hilbert spiral transform. Opt Express 2012;20:23463–79. [15] Trusiak M, Wielgus M, Patorski K. Advanced processing of optical fringe patterns by automated selective reconstruction and enhanced fast empirical mode decomposition. Opt Laser Eng 2014;52:230–40. [16] Wu Z, Huang NE. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Adv Adapt Data Anal 2009;1:1–41. [17] Zhou Y, Li H. Adaptive noise reduction method for DSPI fringes based on bidimensional ensemble empirical mode decomposition. Opt Express 2011;19:18207–15. [18] Wu ZH, Huang NE, Chen X. The multi-dimensional ensemble empirical mode decomposition method. Adv Adapt Data Anal 2009;1:339–72. [19] Rehman NU, Park C, Huang NE, Mandic DP. EMD via MEMD: Multivariate noise-aided computation of standard EMD. Adv Adapt Data Anal 2013;5 (1350007):1–25. [20] Torres ME, Colominas MA, Schlotthauer G, Flandrin P. A complete ensemble empirical mode decomposition with adaptive noise. In: Proceedings of IEEE ICASSP; 2011. p. 4144–7. [21] Colominas MA, Schlotthauer G, Torres ME. Improved complete ensemble EMD: a suitable tool for biomedical signal processing. Biomed Signal Process Control 2014;14:19–29. [22] Wang C, Da F. Phase retrieval for noisy fringe pattern by using empirical mode decomposition and Hilbert Huang transform. Opt Eng 2012;51(061306):1–12. [23] Wang C, Da F. Differential signal-assisted method for adaptive analysis of fringe pattern. Appl Opt 2014;53:6222–9. [24] Equis S, Jacquot P. The empirical mode decomposition: a must-have tool in speckle interferometry. Opt Express 2009;17:611–23. [25] Su W, Lee CK, Lee CW. Noise-reduction for fringe analysis using the empirical mode decomposition with the generalized analysis model. Opt Laser Eng 2010;48:212–7. [26] Bernini MB, Federico A, Kaufmann GH. Normalization of fringe patterns using the bidimensional empirical mode decomposition and the Hilbert transform. Appl Opt 2009;48:6862–9. [27] Osman S, Wang W. An enhanced Hilbert–Huang transform technique for bearing condition monitoring. Meas Sci Technol 2013;24(085004):1–13. [28] Trusiak M, Patorski K, Pakorski K. Hilbert–Huang processing for single-exposure two- dimensional grating interferometry. Opt Express 2013;21:28359– 79. [29] Zhou X, Podoleanu AG, Yang Z, Yang T, Zhao H. Morphological operation-based bi-dimensional empirical mode decomposition for automatic background removal of fringe pattern. Opt Express 2012;20:24247–62. [30] Zhou Y, Li H. A denoising scheme for DSPI fringes based on fast bi-dimensional ensemble empirical mode decomposition and BIMF energy estimation. Mech Syst Signal Process 2013;35:369–82. [31] Sharpley RC, Vatchev V. Analysis of intrinsic mode functions. Constr Approx 2006;21:17–47. [32] Wang C, Qian K, Da F. Regenerated phase-shifted sinusoid-assisted empirical mode decomposition. IEEE Signal Process Lett 2016;23:556–60. [33] Stevenson N, Mesbah M, Boashash B. A sampling limit for empirical mode decomposition. In: Proceedingsof the Eighth International Symposium on Signal Process and its Applications, vols 1–2; 2005. p. 647–50. [34] Rilling G, Goncalves Flandrin P. On empirical mode decomposition and its algorithm. IEEE-EURASIP Workshop on Nonlinear Signal Image Process, Grado, Italy; June 2003. p. 8–11. [35] 〈http://rcada.ncu.edu.tw/research1.htm〉. [36] 〈http://bioingenieria.edu.ar/grupos/ldnlys〉. [37] Wang Z, Bovik AC. A universal image quality index. IEEE Signal Process Lett 2002;9:81–4. [38] 〈http://gr.xjtu.edu.cn/web/zhouxiang〉.

Please cite this article as: Wang C, et al. Regenerated phase-shifted sinusoids assisted EMD for adaptive analysis of fringe patterns. Opt Laser Eng (2016), http://dx.doi.org/10.1016/j.optlaseng.2016.04.018i