Full-field vibration profilometry using time-averaged interference microscopy aided by variational analysis

Maria Cywińska; Maciej Trusiak; Adam Styk; Krzysztof Patorski

doi:10.1364/OE.28.000435

1. Introduction

Full-field vibration profilometry is a relevant issue in the case of microelectromechanical system (MEMS or MOEMS) manufacturing. In the case when micromechanical component (i.e., micromembrane or microcantilever) is excited to harmonic oscillatory motion, its vibration amplitude distribution can provide the information about material properties in microscale. Vibration analysis can also help with validation of the design and simulation or the manufacturing process optimization. The full-field optical measurement techniques allowing for non-invasive measuring the vibration amplitude distribution are stroboscopic and time-averaged interferometry (TAI) [1–8]. In this paper we focus on advancing the second approach applied to MEMS elements with specular reflection. We propose the algorithmic solution for time-averaged interference microscopy called Hilbert-Variational Vibration Amplitude Demodulation (HV2AD) enabling complete end-to-end numerical analysis and reconstruction of the spatial distribution of vibration amplitude.

As we are dealing with interferometry the measurement result is received in the form of the interferogram (fringe pattern). The difference between conventional interferometry and the time-averaged one is that in the case of the first one (majority of interferometric measurements) the measurand is encoded in the phase function (optical path difference), while in the time-averaged interferometry the modulation of the interferogram stores information about investigated object harmonic oscillations (vibration amplitude distribution). In the case of classical two-beam interferometry and measurement of static objects with specular reflection the intensity distribution of the interferogram can be written as [2]:

(1)$${I_{stat}}({x, y} )= K({x, y} )\{{1 + {C_{stat}}({x, y} ) \cos ({{\varphi_{stat}}({x, y} )} )} \},$$

where K(x, y) describes background intensity, C_stat(x, y) describes fringe contrast distribution and φ_stat(x, y) denotes phase distribution due to the initial object surface shape (its deviation from flatness). In the case of the object harmonic sinusoidal vibration and its TAI measurement, where the recording time is much longer than the induced vibration period, the interferogram intensity distribution can be described as [2]:

(2)$${I_{vibr}}({x, y} )= K({x, y} )\left\{ {1 + {C_{stat}}({x, y} ){J_0}\left( {\frac{{4\pi }}{\lambda }{a_0}({x, y} )} \right)\cos ({{\varphi_{vibr}}({x, y} )} )} \right\},$$

where J₀ denotes Bessel function of zero order (first kind) which is modulating the amplitude of the TAI interferogram; Bessel function argument a₀(x, y) denotes the spatial distribution of vibration amplitude – the information we are mainly interested in vibration testing. If the average position of the vibrating object is the same as its shape in the static state then φ_vibr(x, y)= φ_stat(x, y).

The complete process of information retrieval in the case of TAI measurement can be divided in two steps:

(1) evaluation of the Bessel function (obtaining the Besselogram, i.e., TAI amplitude modulation distribution), which comes down to amplitude demodulation of the recorded interferogram;
(2) estimation of the spatial distribution of vibration amplitude, which comes down to the Bessel function argument calculation or in other words the Besselogram phase retrieval.

The classical approach to the problem of assessing the spatial distribution of the vibration amplitude is to use the temporal phase shifting algorithm [9–11]. To employ this method one needs to introduce phase shift (Δφ₁, Δφ₂) both to the TAI-recorded interferogram (Eq. (3)) and Besselogram (Eq. (4)):

(3)$${{I^{\prime}}_{vibr}}({x, y} )= K({x, y} )\left\{ {1 + {C_{stat}}({x, y} ){J_0}\left( {\frac{{4\pi }}{\lambda }{a_0}({x, y} )} \right)\cos ({{\varphi_{vibr}}({x, y} )+ \Delta {\varphi_1}} )} \right\},$$

(4)$${{I^{\prime\prime}}_{vibr}}({x, y} )= K({x, y} )\left\{ {1 + {C_{stat}}({x, y} ){J_0}\left( {\frac{{4\pi }}{\lambda }{a_0}({x, y} )+ \Delta {\varphi_2}} \right)\cos ({{\varphi_{vibr}}({x, y} )} )} \right\}.$$

First phase shifting is used to calculate the interferogram amplitude distribution – the Besselogram. Then, second phase shifting formula is used to calculate the Besselogram phase values connected with the vibration amplitude. Although the phase shifting method is considered the most accurate one in the case of fringe pattern analysis, one needs to record at least 9 interferograms to employ 3-frame Besselogram phase shifting and 25 interferograms in case of 5-frame Besselogram phase shifting. It should be noted, however, the former case is realized quicker experimentally (less data acquired), the later approach provides considerably more accurate results in both aforementioned steps [9,12,13]. Another major requirement associated with the use of the phase-shifting in the case of TAI-measurement is a specialized optical setup hardware – phase shifting modules are needed. Additionally, due to the difference between the Bessel function and cosine function the phase correction routine is necessary [11]. On the other hand the challenging task of the Besselogram estimation and demodulation can be solved, if required conditions are fulfilled, with the use of the single-frame automatic fringe pattern analysis techniques, for example, spatial carrier phase shifting (SCPS) [14], Fourier transform (FT) [15], continuous wavelet transform (CWT) [16] and the Hilbert-Huang transform (HHT) [17].

Many alternatives to the temporal phase shifting in Bessel pattern analysis were presented. Most of them are adapted from classical fringe pattern processing techniques. For this purpose Regional Inversion of Bessel Function [18], Hilbert Transform [19], Directional Spatial Carrier Phase Shifting [20], Peak Direct Inverse and Optimization Method [21], Derivative-based regularized phase tracker model [22] or most recently Linear Regression [23] have been proposed. An interesting idea, similar up to some extent with [18], is tracing Bessel fringe centers representing equidistant contours of equal amplitude. The map of vibration amplitude can be obtained by interpolating the data between the skeletonized lines. However, this straightforward solution can be difficult or even impossible in some cases. If the vibration amplitude is big enough high order Bessel fringes become very dark. Some solution to that problem was proposed with the use of two-frame Hilbert transform processing [24]. After mutual subtraction of both frames resulted pattern is spoiled with constant non-zero background, which produces normalization errors precisely in the location of dark Bessel fringes (contrast reversal zones of carrier cosinusoidal fringes). Subtracting the results of unspoiled and spoiled pattern normalization contouring of dark Bessel fringes can be obtained.

Another approach to the TAI-recorded interferogram analysis is to use the Hilbert-Huang transform both for Besselogram and the vibration amplitude distribution estimation steps [25] The two-dimensional empirical mode decomposition algorithm is deployed for operator-defined pre-filtering of both interferograms and the Besselogram, increasing their signal to noise ratio and preparing for efficient calculation of 2D Hilbert transform, employed for amplitude and phase demodulation (of filtered interferogram and filtered Besselogram, respectively). In this paper, we propose Hilbert-Variational Vibration Amplitude Demodulation (HV2AD) algorithm, where unsupervised variational image decomposition (uVID) [26] and Hilbert spiral transform (HST) [27] are tailored and combined to devise novel end-to-end Besselogram demodulation scheme.

The paper is structured as follows: section 2 introduces proposed HV2AD scheme and describes two steps of the TAI-recorded interferogram analysis, section 3 contains numerical evaluation of the proposed novel HV2AD algorithmic solution for full-field vibration profilometry using simulated and experimental data comparing the obtained results with the reference HHT-based ones; section 4 concludes the paper.

2. Principles of Hilbert-Variational Vibration Amplitude Demodulation

The Hilbert spiral transform (HST) [27], when appropriately implemented, introduces numerically the phase shift equal to π/2 (quadrature transform) between its input s(x, y) and output s_H(x, y) signals - interferograms in our case. It means that if we describe the HST input signal as:

s({x, y} )= b({x, y} )\cos ({p({x, y} )} ),

then the HST output signal can be described as:

{s_H}({x, y} )={-} b({x, y} )\sin ({p({x, y} )} ).

Having those two we can establish the analytic signal in the following form:

{s_A}({x, y} )= s({x, y} )+ i{s_H}({x, y} ),

with easy access to the input signal amplitude (b(x, y)) by calculating the modulus of analytic signal and phase (p(x, y)) by calculating the angle of the analytic signal. However there are some requirements for the HST to perform well: (1) the modulating function b(x, y) needs to be slowly varying (Bedrosian theorem), (2) the input data needs to be zero mean valued, and (3) the noise needs to be minimized. Two latter assumptions can be achieved by successful background and noise removal with the use of interferogram preprocessing algorithm. The example of such an algorithm is unsupervised variational image decomposition (uVID) [26], where the input signal is automatically decomposed into three components: structure (which can be considered as interferogram or Besselogram slowly varying background term), texture (which can be considered as interferogram or Besselogram information carrying fringe-part) and noise. The uVID filtration algorithm may be used for a variety of image characteristics (in the meaning of the spatial frequency complexity) which expands its application area. However the main limitation of uVID is the calculation time. Generally, the more complicated analyzed image is the more inner iterations are needed and therefore the calculation time is longer. Another important aspect in the case of the HST calculation is the local fringe direction map estimation. In our approach the local fringe orientation map is firstly calculated using combined gradient and plane-fitting method [28], and then unwrapped [29] to the local fringe direction map. Based on the HST and uVID algorithmic solutions we develop HV2AD as a tool for full-field vibration profilometry. The whole HV2AD calculation path is presented in Fig. 1 and discussed in sub sections 2.1 and 2.2.

Fig. 1. The calculation scheme of the proposed HV2AD approach.

Download Full Size | PDF

2.1. Interferogram processing to calculate its amplitude modulation distribution (the Besselogram)

The first step of TAI-registered interferogram processing, Fig. 1, is the calculation of its amplitude modulation distribution (Besselogram), which is described by the Bessel function of zero order (first kind). We need to acquire two phase-shifted TAI-registered interferograms with phase shift according to Eq. (3). The input pattern used for further processing is obtained by their subsequent subtraction. The main advantage of this approach is the amplitude modulation augmentation dependent on phase-shift value (see Eq. (5)). If we assume that background is stable (not changing in time) during the phase-shifting registration then the subtraction result can be described as (background can be simply neglected by subtraction of two phase-shifted interferograms):

(5)$${I_{subtr}}({x, y} )= 2 \cdot {C_{stat}}({x, y} ){J_0}\left( {\frac{{4\pi }}{\lambda }{a_0}({x, y} )} \right)\sin \left( {\frac{{\Delta \varphi }}{2}} \right)\sin \left( {{\varphi_{vibr}}({x, y} )+ \frac{{\Delta \varphi }}{2}} \right),$$

where Δφ describes phase-shift value between two registered interferograms. It can be clearly seen that the most suitable Δφ value is π, but significantly larger regime (π/3:4π/3) is profitable [24,30–32]. The influence of different phase step values was already analyzed in [25] and in this work, for the numerical analysis, we simply chose the best possible phase-shift value without reproducing the computations presented in [25]. The subtraction step can be considered as physical (experimental) interferogram signal-to-noise ratio improvement. In the HV2AD procedure, in order to improve SNR even more, we added the uVID filtration step. This way any residual background and noise can be successfully removed, which has crucial impact on the robustness of the proposed approach. Having the zero mean value signal calculated the Hilbert spiral transform can be computed and the analytic signal (analytic-interferogram) can be established successfully. Accessing its amplitude term by calculating its modulus the Besselogram is obtained. The Besselogram in turn is phase modulated by the spatial distribution of the vibration amplitude of interest.

2.2. Besselogram processing to calculate its argument for obtaining the vibration amplitude distribution

The next step of the HV2AD processing path is connected with the Besselogram analysis, Fig. 1. Firstly, we need to remove its background component in order to successfully employ the HST. To serve this purpose the uVID filtration is used. The additional advantage of this approach is error minimization connected with classifying some errors introduced in the previous HV2AD steps as noise by another uVID procedure. Similarly to the first step of interfeogram processing the HST is computed and now the analytic-Besselogram is established. Vibration amplitude distribution is encoded in the Bessel function argument thus the angle (phase) of analytic-Besselogram needs to be obtained. After unwrapping [29] the HV2AD scheme is finalized and the vibration amplitude distribution is performed.

It is worth to notice that we are using the algorithmic solution (HST) especially tailored for dealing with cosinusoidal signal while Bessel function can be considered as highly non-cosinusioidal one. In consequence the fringe like error appearance in the resultant distribution of vibration amplitude is observed [19,33,34]. The problem will be discussed in the next section and the numerical error correction method will be proposed.

3. Experiments and results

In this section the numerical analyses of the proposed HV2AD algorithmic solution are conducted. The first and the second step of the calculation path are discussed separately employing both simulated and experimental data.

3.1. Numerical evaluation

Initially we are going to find how the additional uVID filtration of input pattern is significant for the further data processing or whether it can be neglected due to previous subtraction of two phase-shifted time-averaged interferograms. In Fig. 2 the results of Besselogram determination (TAI-registered interferogram amplitude demodulation) without uVID filtration are presented. Two cases with low and high vibration amplitude values were analyzed. The higher the vibration amplitude value the denser the Bessel fringes are. In general single-frame fringe pattern analysis techniques cope better with dense fringes [35–41]. There is some limit to Bessel fringe frequency, however, as with the increase of the measured vibration amplitude the contrast of Bessel fringes decreases (signal to noise ratio is lower). Additionally we need to remember that according to the Bedrosian’s theorem for Hilbert transform the modulation function needs to be slowly varying, which means it has no frequency content above the TAI-registered interferogram carrier frequency. In order to measure the accuracy of estimated Bessel patterns the RMS error values were calculated. The RMS value for the result presented in Fig. 2(d) is equal to 0.056 (no uVID pre-filtration) and for the result presented in Fig. 2(i) is equal to 0.05 (no uVID pre-filtration). In Fig. 3 the Besselograms estimated with uVID pre-filtration are presented. One can notice that by simply adding the uVID filtration step, the estimated results may be significantly improved (compare error maps in Figs. 2(e) and 2(j) with Figs. 3(e) and 3(j)). The RMS values can be lowered employing the uVID pre-filtration to 0.01 for the result presented in Fig. 3(d) and 0.013 for the result presented in Fig. 3(i). Pre-filtration can both minimize noise and uneven background and therefore it is an important calculation step for accurate interferogram amplitude demodulation (Besselogram calculation).

Fig. 2. Interferogram amplitude estimation without uVID prefiltration: (a),(b),(f),(g) two π-phase-shifted interferograms, (c),(h) input patterns obtained by subtraction of two interferograms, (d),(i) estimated interferogram amplitude distributions and (e),(j) error maps; upper row: lower vibration amplitude, lower row: higher vibration amplitude, colorbar: arbitrary units.

Download Full Size | PDF

Fig. 3. Interferogram amplitude estimation with uVID prefiltration: (a),(f) residual background terms, (b),(g) filtered fringes, (c),(h) noise components, (d),(i) estimated interferogram amplitude distributions and (e),(j) error maps; upper row: lower vibration amplitude, lower row: higher vibration amplitude, colorbar: arbitrary units.

Download Full Size | PDF

Although it was proven that the uVID pre-filtration positively affects the retrieved Besselogram accuracy the ultimate goal of TAI-registered interferogram analysis is to calculate the vibration amplitude distribution from the argument of the Bessel function. In Figs. 4 and 5 the vibration amplitude demodulation results are presented for lower and higher vibration amplitude values, respectively. The proposed HV2AD calculation path is compared with the Hilbert-Huang transform (HHT) [25], which has been already successfully used for both interferogram and Besselogram analysis. Analyzing the Besselogram filtration results one can notice some errors in the corners of the calculated background term. However, these errors do not transfer to the vibration amplitude map, since using the HST we are calculating the angle (phase) of the Bessel function not its amplitude modulation. Unification of the HST input Besselogram contrast modulation may be considered as an advantage. In the case of the lower vibration amplitude the RMS error value is equal to 0.3 for HHT (Fig. 4(d)) and 0.2 for HV2AD (Fig. 4(i)). The difference between two methods is more visible in the case of the higher vibration amplitude, where the RMS error value is equal to 1.27 for HHT (Fig. 5(d)) and 0.27 for HV2AD (Fig. 5(i)). Exceptionally high RMS error value in the case of the HHT calculation path is connected with discontinuities introduced in the unwrapping step and clearly visible in the error map (Fig. 5(e)). This effect is caused by the noise present in the analyzed Besselogram. Additionally, due to the high vibration amplitude, the contrast of Bessel fringes in the central part of Besselogram (around maximum of vibration amplitude) is low, which in combination with high noise significantly lowers the signal to noise ratio. This is another argument proving why successful noise removal is important in the case of TAI-registered interferogram analysis.

Fig. 4. Besselogram phase estimation – lower vibration amplitude case: (a), (f) estimated Besselogram background terms, (b), (g) filtered Bessel patterns, (c), (h) noise components, (d), (i) estimated vibration amplitude distributions and (e), (j) error maps; upper row: the HHT calculation path, lower row: the HV2AD calculation path, colorbar: (a-c), (f-h) arbitrary units, (d),(e),(i),(j) radians.

Download Full Size | PDF

Fig. 5. Besselogram phase estimation – higher vibration amplitude case: (a), (f) estimated Besselogram background terms, (b), (g) filtered Bessel patterns, (c), (h) noise components, (d), (i) estimated vibration amplitude distributions and (e), (j) error maps; upper row: the HHT calculation path, lower row: the HV2AD calculation path, colorbar: (a-c), (f-h) arbitrary units, (d),(e),(i),(j) radians.

Download Full Size | PDF

The Bessel fringes contrast decrease is connected with our previous statement that if the vibration amplitude is big enough high order Bessel fringes become very dark, which makes the task of amplitude demodulation very challenging (see Figs. 2(i), 3(i) and 5). Nevertheless in both cases (low and high vibration amplitude) similarly for the HHT and the HV2AD the fringe-like error in obtained phase map is clearly visible. As we mentioned before it is connected with highly non-cosinusoidal nature of the Bessel function. In the next subsection the possible solutions to this problem will be discussed and the numerical error correction method will be proposed.

3.2. Phase error minimization

Capable solutions for elimination of phase fringe-like errors caused by noncosinusolidal characteristic of the analyzed function were already proposed in literature. The representatives of mentioned methods are: cubic spline-smooting method [42] or iterative solution for phase-shifting method [43]. However, Bessel fringes cannot be treated only as a sinusoidal pattern with higher harmonics causing signal nonlinearity. Comparing Bessel function J₀ and cosinusoidal function, one may clearly see the discrepancy of the extreme positions of both functions, especially for the small argument values. Moreover, the zeros of Bessel function are not equally spaced. Therefore the error related to aforementioned Bessel function property has to be minimized in the phase-shifting based Besselogram analysis [11,44,45]. During Bessel fringe pattern processing one needs to apply an additional phase correction routine. This routine uses mainly a look-up-table, to correct the obtained phase values and estimate the real ones. This process provides satisfactory results, however it is regarded as time-consuming and cumbersome especially under noisy conditions.

According to [11] argument of the Bessel function (${\Omega }({\textrm{x},\textrm{y}} )= \frac{{4\pi }}{\lambda }{a_0}({x, y} )$) using temporal phase shifting (TPS) scheme may be calculated using the equation:

(6)$${\Omega^{\prime}}({x, y} )= \frac{1}{2}{\tan ^{ - 1}}\left\{ {\left[ {\frac{{1 - \cos ({2B} )}}{{\sin ({2B} )}}} \right]\frac{{J_0^2[{|{{\Omega }({x, y} )+ B} |} ]- J_0^2[{|{{\Omega }({x, y} )- B} |} ]}}{{2J_0^2[{|{{\Omega }({x, y} )} |} ]- J_0^2[{|{{\Omega }({x, y} )- B} |} ]- J_0^2[{|{{\Omega }({x, y} )+ B} |} ]}}} \right\},$$

where B is a Besselogram phase shift value. In our simulations, Fig. 6, the value B=π/3 is used resulting in phase-shifted Besselogram series (Figs. 6(a)–6(c)). As was mentioned before, due to the difference between Bessel and cosine function the phase map calculated with the use of Eq. (6) is flawed with a fringe-like error. The resultant phase and error maps for a TPS algorithm are presented in Fig. 6(d) and Fig. 6(e), respectively. The nature of the fringe-like error is slightly different in the case of single-frame Besselogram analysis. In this paper in order to retrieve both the amplitude modulation of the TAI-registered interferogram and the argument of the Bessel function the Hilbert spiral transform is used. It is worth to notice that after the interferogram amplitude demodulation with the use of HST we are not receiving the Bessel function but its modulus, Fig. 6(f). One of the requirements for HST to perform well is the zero-mean value signal as an input. It means that the DC component needs to be removed from the estimated interferogram amplitude map. The Besselogram HST input may be described as:

(7)$$f({x, y} )= |{{J_0}({{\Omega }({x, y} )} )} |- DC.$$

Fig. 6. Fringe-like error in multi- and single-frame Besselogram analysis; Besselogram phase-shifted series: (a) B = 0, (b) B =π/3, (c) B = -π/3, (d) estimated vibration amplitude modulation by TPS and (e) error map, (f) input pattern (modulus of the Bessel function) for a single-frame Besselogram analysis algorithm, (g) modulus of the Bessel function after DC component removal by uVID, (h) HST result of Fig. 6(g), (i) estimated vibration amplitude modulation by HV2AD and (j) error map, (k) modulus of the Bessel function after DC component removal by empirical mode decomposition, (l) HST result of Fig. 6(k), (m) estimated vibration amplitude modulation by HHT and (n) error map.

Download Full Size | PDF

The argument of the Bessel function using HST may be calculated using the equation:

(8)$${\Omega^{\prime}}({x, y} )= {\tan ^{ - 1}}\left[ {\frac{{HST({f({x, y} )} )}}{{f({x, y} )}}} \right].$$

The results of the Besselogram filtration using uVID and empirical mode decomposition are presented in Fig. 6(g) and Fig. 6(k), respectively. Theoretically the HST should introduce the phase shift equal to π/2 to its input pattern. The results retrieved for the Besselogram analysis are presented in Fig. 6(h) and Fig. 6(l) for HV2AD and HHT, respectively. The resultant phase maps and error maps are presented in Figs. 6(i)–6(j) and Figs. 6(m)–6(n). One can clearly notice that the fringe-like error characteristic is different in the case of multi- and single-frame Besselogram analysis. The RMS error values for the vibration amplitude map demodulation are: 0.19 rad for the TPS, 0.17 rad for the HV2AD and 0.21 rad for the HHT.

Our idea for numerical error correction method is based on an inherent feature of the uVID method (or variational image decomposition in general). It is an excellent algorithmic solution for detecting the oscillatory parts of the image. Since the uVID can successfully distinguish the oscillatory part of the image from the cartoon part of the image it can be used for removing the fringe-like error from the ground-truth phase function (related to the vibration distribution map) we want to determine. The variational-based decomposition is calculated from erroneous phase map we want to correct. The uVID enables the image decomposition into three components: structure – previously considered as interferogram/Besselogram background term, texture – previously considered as interferogram/Besselogram fringes and noise. In our correction method the structure may be considered as ground-truth phase function while texture and noise as oscillatory, fringe-like error. Straightforward removal of those two decomposition components improves calculated result. Similar approach was already applied for the phase enhancement in single-frame fringe pattern analysis [46] and, as we show below, can be successfully adapted to the Bessel-function-based full-field vibration profilometry.

The performance of the proposed solution is presented in Fig. 7 using the simulated data. The vibration amplitude was simulated as a spherical function with the maximum located in the center of the simulation grid. By the change of the vibration amplitude map maximum value the density of retrieved Bessel fringes was changed, where the higher the simulated amplitude the denser the fringes. One can notice that in general the HV2AD results are improving with the increase of the vibration amplitude maximum. It is connected with the fact that the single-frame fringe pattern analysis techniques typically cope better with the dense fringes [35–41]. Basing on the analysis of the results presented in Fig. 7 it was proven that the uVID-based phase correction method can minimize the fringe-like error in the resultant phase map. In order to highlight this feature of our algorithm the error maps for the uncorrected and corrected vibration amplitude are shown in inserts of Fig. 7. The three cases were considered in detail: they correspond to low, medium and high vibration amplitude values. In each of the highlighted results the significant improvement may be observed.

Fig. 7. Relation between simulated vibration amplitude and the RMS error value for the uncorrected and corrected Besselogram demodulation cases, colorbar: radians.

Download Full Size | PDF

3.3. Experimental data processing

After the tests on the simulated data the performance of the proposed HV2AD algorithm is verified in a step-by-step manner exploiting the experimental data. The detailed description of the experimental setup is given in [2]. As an object under test flat square silicon membrane with lateral size of 3.48 × 3.48 mm and thickness of 15 µm, excited externally by PZT transducer, was investigated. For the first vibration mode of the membrane a frequency of 19.8 kHz was found. The investigations of resonant frequencies as well as experimental data acquisition, for further processing with proposed methods, were performed using standard TAI method (no Bessel fringes shifting present). In Fig. 8 the results of Besselogram estimation using HHT and HV2AD are presented. The input pattern (Fig. 8(c)) is calculated by the subtraction of two phase-shifted interferograms (Figs. 8(a) and 8(b)). Firstly, the input pattern amplitude modulation (Fig. 8(d)) is determined with straightforward use of the Hilbert spiral transform, following the HHT calculation scheme. One can notice that the result calculated this way is very noisy, which suggests that the input pattern pre-filtration is needed. The uVID decomposition result is presented in Figs. 8(e)–8(g). After the uVID pre-filtration the amplitude modulation function (Fig. 8(h)) is calculated with the use of the Hilbert spiral transform (Fig. 8(f)), following the HV2AD calculation scheme. This way the quality of the result can be visibly improved.

Fig. 8. Real interferogram amplitude estimation: (a),(b) two π-phase-shifted interferograms, (c) input pattern obtained by the subtraction of two interferograms, (d) estimated interferogram amplitude without uVID prefiltration, (e) residual background, (f) fringes and (g) noise estimated by uVID, (h) estimated interferogram amplitude with uVID prefiltration, colorbar: arbitrary units.

Download Full Size | PDF

The second step of the TAI-registered interferogram is connected with the vibration amplitude map demodulation. The results of the Besselogram demodulation using the HHT and the HV2AD are shown in Fig. 9. The Besselogram from Fig. 8(d) was filtered using empirical mode decomposition (Huang’s part of the HHT) and the result of the filtration may be seen in Figs. 9(a)–9(c), where Fig. 9(a) depicts background, Fig. 9(b) presents the filtered Bessel pattern and Fig. 9(c) shows the noise. It is worth to emphasize that the result of the decomposition using empirical mode decomposition highly depends on the proper classification of retrieved bidimensional intrinsic mode functions (modes) into three categories: background, fringes and noise, while the decomposition by uVID is done automatically. Here, we simply defined the first mode as noise, background was determined with the use of scheme proposed in [47] and the rest of modes was classified as fringes. The vibration amplitude map demodulation calculated by the HHT is presented in Fig. 9(d). In a similar way the Besselogram from Fig. 8(h) was filtered using the uVID algorithm. The uVID-calculated background is shown in Fig. 9(e), the uVID-calculated noise in Fig. 9(g) and the uVID-calculated Bessel pattern in Fig. 9(f). The resultant HV2AD vibration amplitude map is presented in Fig. 9(h). It can be seen that the result estimated by the HHT is more noisy than the one estimated by the HV2AD. Additionally there is some error present in the central part of the image in the case of the HHT-based vibration amplitude map caused by the unwrapping algorithm noise-induced failure.

Fig. 9. Real Besselogram phase estimation: (a), (e) estimated Besselogram background terms, (b), (f) filtered Bessel patterns, (c), (g) noise components and (d), (h) estimated vibration amplitude distributions; upper row: the HHT calculation path, lower row: the HV2AD calculation path, colorbar: (a-c), (e-g) arbitrary units, (d),(h) radians.

Download Full Size | PDF

In Figs. 8(e)–8(h) and Figs. 9(e)–9(g) the values in the colorbar are significantly smaller comparing with Figs. 8(a)–8(d) and Figs. 9(a)–9(c), because initial image pixel values needed to be normalized to the range <0,1 > before calculating image decomposition using uVID calculation scheme. However, this operation does not affect the vibration amplitude range and therefore quality of the measurement.

In Fig. 10 the results of the proposed fringe-like error correction method are presented. The improvement in the vibration amplitude map quality may be observed comparing the HV2AD result before (Fig. 10(a)) and after correction (Fig. 10(b)) or the HHT result before (Fig. 10(d)) and after correction (Fig. 10(e)). In both cases in the part removed during the correction scheme (Figs. 10(c) and 10(f)) the fringe-like error and some discontinuities introduced upon unwrapping step are present. Additionally, in the case of the HHT-based result, the noise was also minimized thanks to the uVID phase filtration. It is also worth to notice that the gain in the form of the error minimization was not accomplished at a cost of any information loss. The part removed during the correction oscillates around zero and therefore no part of vibration amplitude was removed. Lack of the vibration amplitude loss can be also observed in the cross sections of the results presented in Fig. 10(g).

Fig. 10. Vibration amplitude map (Besselogram phase map) correction results: vibration amplitude map calculated by the HV2AD (a) before uVID correction, (b) after uVID correction, (c) difference between Fig. 9(a) and Fig. 9(b); vibration amplitude map calculated by the HHT (d) before uVID correction, (e) after uVID correction, (f) difference between Fig. 9(d) and Fig. 9(e); and (g) cross sections of the results by the central point, colorbar: radians.

Download Full Size | PDF

4. Conclusions

The main novelty of this paper resides in a new algorithmic solution for coherent time-averaged vibration amplitude map determination called Hilbert-Variational Vibration Amplitude Demodulation (HV2AD). Using HV2AD one can successfully estimate the full-field vibration amplitude map having only two phase-shifted TAI-registered interferograms. Described solution directly applicable to MEMS elements with specular reflection compares favorably with the already well-developed and efficient Hilbert-Huang transform [25]. The improvement of the retrieved vibration amplitude map accuracy was accomplished thanks to additional interferogram pre-filtering by unsupervised variational image decomposition. Successful noise removal prevents unwrapping algorithm failure, eliminating possibly jeopardizing discontinuities in the resultant vibration amplitude map. High level of noise may also make it impossible to retrieve any useful information using Hilbert spiral transform. It is also very important to remove residual nonzero background, which may appear in HV2AD input pattern, e.g., due to the change of illumination during two interferograms recording [24]. Uneven background introduces the fringe like error to the Hilbert transform-based Besselogram demodulation.

The same set of novel tools is applied to pre-filter the Besselogram and calculate its phase distribution representing the vibration amplitude map of interest. Another novel aspect is related to the uVID algorithm employed for the correction of the fringe-like errors in the Bessel-function-based full-field vibration profilometry. These errors are induced by strong discrepancy between the Bessel function and the harmonic (cosine) one. The usefulness and efficiency of the proposed end-to-end HV2AD method was corroborated utilizing the simulated and experimental data.

We consider the following aspects as main novelties of the paper:

(1) a new automatic end-to-end data processing path for time-averaged interferometry with limited requirements on recorded interferograms is proposed,
(2) the use of a posteriori automatic filtration for Besselogram phase map is proposed,
(3) reported approach provides novel results which quantitatively and qualitatively outperform the ones obtainable by the Hilbert-Huang transform approach [25] and are numerically comparable with results of the phase shifting technique (with significantly increased number of interferograms/Besselograms; their obtaining complicates the setup).

It is important to note that in the HHT method [25] the pre-filtering, thus the selection of modes (bidimensional intrinsic mode functions), was performed based on the visual judgement of the operator. This approach can be automated for a certain class of predefined patterns, nonetheless, in general it lacks full adaptivity. The HV2AD works completely automatically taking advantage of the adaptive Chambolle projection algorithm which fuels the variational analysis. In our previous work [41] we have presented automated approach tailored mainly for consinusoidal fringes pre-filtering to enhance their subsequent phase demodulation; in this contribution, for the first time to the best of our knowledge, fully automated variational approach is presented for consinusoidal fringes pre-filtering to enhance their amplitude demodulation and Bessel fringes pre-filtering to augment their phase demodulation. Additional a posteriori Bessel phase map filtering is to be emphasized as well.

In the future work we will focus on the acceleration of the proposed algorithmic solution. We plan to conduct the separate comparative studies considering discussion of mixing the HV2AD and HHT algorithms for interferogram and Besselogram analysis in order to make our calculation path versatile and robust. It is possible that one algorithm would perform better in the case of interferogram analysis while the other in the case of Besselogram analysis. Codes for HV2AD will be made available upon reasonable request.

Funding

Narodowe Centrum Nauki (OPUS 13 2017/25/B/ST7/02049); Politechnika Warszawska ((Dean funds), Faculty of Mechatronics statutory funds.

Acknowledgment

We acknowledge the support of Polish National Agency for Academic Exchange funds.

References

1. A. Bosseboeuf and S. Petitgrand, “Application of microscopic interferometry techniques in the MEMS field,” Proc. SPIE 5145, 1–16 (2003). [CrossRef]

2. L. Sałbut, K. Patorski, M. Józwik, J. Kacperski, C. Gorecki, A. Jacobelli, and T. Dean, “Active micro-elements testing by interferometry using time-average and quasi-stroboscopic techniques,” Proc. SPIE 5145, 23–32 (2003). [CrossRef]

3. P. Picart, J. Leval, D. Mounier, and S. Gougeon, “Time-averaged digital holography,” Opt. Lett. 28(20), 1900–1902 (2003). [CrossRef]

4. C. Gorecki, M. Jozwik, and L. Salbut, “Multifunctional interferometric platform for on-chip testing the micromechanical properties of MEMS/MOEMS,” J. Micro/Nanolithogr., MEMS, MOEMS 4(4), 041402 (2005). [CrossRef]

5. E. Hong, S. Trolier-McKinstry, R. Smith, S. V. Krishnaswamy, and C. B. Freidhoff, “Vibration of micromachined circular piezoelectric diaphragms,” IEEE Trans. Sonics Ultrason. 53(4), 697–706 (2006). [CrossRef]

6. M. Olfatnia, V. R. Singh, T. Xu, J. M. Miao, and L. S. Ong, “Analysis of the vibration modes of piezoelectric circular microdiaphragms,” J. Micromech. Microeng. 20(8), 085013 (2010). [CrossRef]

7. S. Ellingsrud and G. O. Rosvold, “Analysis of data-based TV-holography system used to measure small vibration amplitudes,” J. Opt. Soc. Am. A 9(2), 237–251 (1992). [CrossRef]

8. T. Statsenko, V. Chatziioannou, T. Moore, and W. Kausel, “Deformation reconstruction by means of surface optimization. Part I: Time-averaged electronic speckle pattern interferometry,” Appl. Opt. 56(3), 654–661 (2017). [CrossRef]

9. K. Patorski and A. Styk, “Interferogram intensity modulation calculations using temporal phase shifting: error analysis,” Opt. Eng. 45(8), 085602 (2006). [CrossRef]

10. M. Kirkove, S. Guérit, L. Jacques, C. Loffet, F. Languy, J.-F. Vandenrijt, and M. Georges, “Determination of vibration amplitudes from binary phase patterns obtained by phase-shifting time-averaged speckle shearing interferometry,” Appl. Opt. 57(27), 8065–8077 (2018). [CrossRef]

11. R. J. Pryputniewicz and K. Stetson, “Measurement Of Vibration Patterns Using Electro-Optic Holography,” Proc. SPIE 1162, 456–467 (1990). [CrossRef]

12. A. Styk and M. Józwik, “Recent Advances in Phase Shifted Time Averaging and Stroboscopic Interferometry,” Proc. SPIE 9960, 99600A (2016). [CrossRef]

13. A. Styk and H. Dziubecka, “Phase shift strategies in phase shifting time averaging interferometry for harmonic motion measurements,” Proc. SPIE 10834, 108342C (2018). [CrossRef]

14. A. Styk and K. Patorski, “Analysis of systematic errors in spatial carrier phase shifting applied to interferogram intensity modulation determination,” Appl. Opt. 46(21), 4613–4624 (2007). [CrossRef]

15. S. Petitgrand, R. Yahiaoui, A. Bosseboeuf, and K. Danaie, “Quantitative time-averaged microscopic interferometry for micromechanical device vibration mode characterization,” Proc. SPIE 4400, 51–60 (2001). [CrossRef]

16. K. Pokorski and K. Patorski, “Visualization of additive-type moiréand time-average fringe patterns using the continuous wavelet transform,” Appl. Opt. 49(19), 3640–3651 (2010). [CrossRef]

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

18. D. N. Borza, “Full-field vibration amplitude recovery from high-resolution time-averaged speckle interferograms and digital holograms by regional inverting of the Bessel function,” Opt. Laser Eng. 44(8), 747–770 (2006). [CrossRef]

19. U. P. Kumar, N. K. Mohan, and M. P. Kothiyal, “Time average vibration fringe analysis using Hilbert transformation,” Appl. Opt. 49(30), 5777–5786 (2010). [CrossRef]

20. A. Styk and M. Brzezinski, “Vibration amplitude recovery from time averaged interferograms using the directional spatial carrier phase shifting metod,” Proc. SPIE 8082, 80821X (2011). [CrossRef]

21. T. Statsenko, V. Chatziioannou, T. Moore, and W. Kausel, “Methods of phase reconstruction for time-averaging electronic speckle pattern interferometry,” Appl. Opt. 55(8), 1913–1919 (2016). [CrossRef]

22. B. Deepan, C. Quan, and C. J. Tay, “Quantitative vibration analysis using a single fringe pattern in time-average speckle interferometry,” Appl. Opt. 55(22), 5876–5883 (2016). [CrossRef]

23. P. A. A. M. Somers and N. Bhattacharya, “A new Method for Processing Time Averaged Vibration Patterns: Linear Regression,” Strain 52(4), 264–275 (2016). [CrossRef]

24. K. Patorski and M. Trusiak, “Highly contrasted Bessel fringe minima visualization for time-averaged vibration profilometry using Hilbert transform two-frame processing,” Opt. Express 21(14), 16863–16881 (2013). [CrossRef]

25. M. Trusiak, A. Styk, and K. Patorski, “Hilbert–Huang transform based advanced Bessel fringe generation and demodulation for full-field vibration studies of specular reflection micro-objects,” Opt. Laser Eng. 110, 100–112 (2018). [CrossRef]

26. M. Cywińska, M. Trusiak, and K. Patorski, “Automatized fringe pattern preprocessing using unsupervised variational image decomposition,” Opt. Express 27(16), 22542–22562 (2019). [CrossRef]

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

28. X. Yang, Q. Yu, and S. Fu, “A combined method for obtaining fringe orientations of ESPI,” Opt. Commun. 273(1), 60–66 (2007). [CrossRef]

29. M. A. Herráez, D. R. Burton, M. J. Lalor, and M. A. Gdeisat, “Fast two-dimensional phase-unwrap- ping algorithm based on sorting by reliability following a noncontinuous path,” Appl. Opt. 41(35), 7437–7444 (2002). [CrossRef]

30. J. Na, W. J. Choi, E. S. Choi, S. Y. Ryu, and B. H. Lee, “Image restoration method based on Hilbert transform for full-field optical coherence tomography,” Appl. Opt. 47(3), 459–466 (2008). [CrossRef]

31. M. S. Hrebesh, “Full-field and single shot full-field optical coherence tomography: a novel technique for biomedical imaging applications,” Adv. Opt. Technol. 2012, 435408 (2012). [CrossRef]

32. Z. Sunderland, K. Patorski, and M. Trusiak, “Subtractive two-frame three-beam phase-stepping interferometry for testing surface shape of quasi-parallel plates,” Opt. Express 24(26), 30505–30513 (2016). [CrossRef]

33. Z. Cai, X. Liu, H. Jiang, D. He, X. Peng, S. Huang, and Z. Zhang, “Flexible phase error compensation based on Hilbert transform in phase shifting profilometry,” Opt. Express 23(19), 25171–25181 (2015). [CrossRef]

34. Y. Wang, Z. Liu, C. Jiang, and S. Zhang, “Motion induced phase error reduction using a Hilbert transform,” Opt. Express 26(26), 34224–34235 (2018). [CrossRef]

35. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]

36. Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” Appl. Opt. 43(13), 2695–2702 (2004). [CrossRef]

37. L. R. Watkins, S. M. Tan, and T. H. Barnes, “Determination of interferometer phase distributions by use of wavelets,” Opt. Lett. 24(13), 905–907 (1999). [CrossRef]

38. M. Servin, J. L. Marroquin, and F. J. Cuevas, “Demodulation of a single interferogram by use of a twodimensional regularized phase-tracking technique,” Appl. Opt. 36(19), 4540–4548 (1997). [CrossRef]

39. M. Kujawinska and J. Wojciak, “Spatial-carrier phase-shifting technique of fringe pattern analysis,” Proc. SPIE 1508, 61–67 (1991). [CrossRef]

40. X. Zhu, Z. Chen, and C. Tang, “Variational image decomposition for automatic background and noise removal of fringe patterns,” Opt. Lett. 38(3), 275–277 (2013). [CrossRef]

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

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

43. B. Pan, Q. Kemao, L. Huang, and A. Asundi, “Phase error analysis and compensation for nonsinusoidal waveforms in phase-shifting digital fringe projection profilometry,” Opt. Lett. 34(4), 416–418 (2009). [CrossRef]

44. S. H. Baik, S. K. Park, C. J. Kimi, and S. Y. Kim, “Analysis of phase measurement errors in electro-optic holography,” Opt. Rev. 8(1), 26–31 (2001). [CrossRef]

45. U. P. Kumar, Y. Kalyani, N. K. Mohan, and M. P. Kothiyal, “Time-average TV holography for vibration fringe analysis,” Appl. Opt. 48(16), 3094–3101 (2009). [CrossRef]

46. M. Cywińska, M. Trusiak, C. Zuo, and K. Patorski, “Enhancing single-shot fringe pattern phase demodulation using advanced variational image demposition,” J. Opt. 21(4), 045702 (2019). [CrossRef]

47. D. Saide, M. Trusiak, and K. Patorski, “Evaluation of adaptively enhanced two-shot fringe pattern phase and amplitude demodulation methods,” Appl. Opt. 56(19), 5489–5500 (2017). [CrossRef]

Full-field vibration profilometry using time-averaged interference microscopy aided by variational analysis

Abstract

1. Introduction

2. Principles of Hilbert-Variational Vibration Amplitude Demodulation

2.1. Interferogram processing to calculate its amplitude modulation distribution (the Besselogram)

2.2. Besselogram processing to calculate its argument for obtaining the vibration amplitude distribution

3. Experiments and results

3.1. Numerical evaluation

3.2. Phase error minimization

3.3. Experimental data processing

4. Conclusions

Funding

Acknowledgment

References

Cited By

Figures (10)

Equations (11)

Optics Express