1. INTRODUCTION

Quantification of non-harmonic deformations by optical experimental methods has been an attractive research topic for decades. One process that can accomplish this is electronic speckle pattern interferometry (ESPI), which is often used for full-field dynamic deformation measurements [1]. ESPI captures whole deforming surfaces in a single recording and provides vibrational analysis with high spatial resolution, allowing real-time measurements [2–4]. Furthermore, ESPI can be applied to visualize transient events in the deforming surfaces [5].

Early time-resolved ESPI (TRESPI) arrangements involved stroboscopic and pulsed laser illumination, which required the repetition of the observed event. Using a high-speed camera for the ESPI measurements is more suitable for non-repeatable events [5,6], and initially the process was based on phase-shifting interferometry, which introduces temporal phase modulation into the experimental arrangement [7]. This method requires a series of phase-shifted images and may not be robust enough to analyze dynamic phenomena [8]. Another approach, spatial phase shifting [9,10], is a promising tool, which introduces spatial shift between two interfering light beams, but yet does not tackle high-speed events.

Recent advances in TRESPI involve advanced mathematical tools such as Hilbert transform [11,12], Fourier transform [13], and empirical mode decomposition (EMD) [14]. Hilbert transform and EMD can be combined for phase extraction in dynamic measurements [13,15], but this method does not perform well when the deformation is non-monotonic [8,16]. EMD is a powerful but not fully developed technique which has several drawbacks, such as non-uniqueness of the decomposition and boundary issues. An a priori knowledge is required to resolve sign ambiguity and to process the images spatially and temporally at the same time [17].

Most of the existing phase recovery algorithms used with ESPI interferograms are aimed at direct processing of the pixel intensity to solve the inverse problem, i.e., using the experimental results to reconstruct the parameters that describe the system under investigation. Pixel processing may be temporal, spatial or spatiotemporal [18]. Another method is to apply a surface optimization algorithm. Surface optimization has been successfully applied by the authors for the case of time-averaged ESPI (TAESPI) in the first part of the paper [19]. It was shown that the intensities can be calculated from the predicted surface deformation shape. The goal of the optimization is to iteratively change the calculated deformation profile to obtain the minimum of an objective function, which characterizes the mismatch between the calculated and the recorded intensity. This is the second part of the two-part investigation of deformation reconstruction by means of surface optimization.

In the past 10 years, ESPI methods have been widely used in the area of music acoustics [20,21]. One of the applications is the quality control for manufacturers. Another application is comparison between musical instruments with equal geometrical characteristics but made from different materials. The scope of the study goes beyond music acoustics and targets research areas in which the reproducibility of the experiments that involve transient response can be complex, such as automotive design.

In what follows we describe the TRESPI methodology for out-of-plane deformation recovery. The aim of this work is to recover the temporal out-of-plane deformation of the surface of a Brazilian cuíca drum under various types of excitation from recorded high-speed interferometric images. The method does not apply phase unwrapping for a single differential frame nor does it deal with the whole frame sequence. It is an iterative method with a spatial processing algorithm intended for phase recovery with minimal cumulative error. The computation is performed using MATLAB (MATLAB and Image Processing Toolbox R2014a, The MathWorks, Inc., Natick, Massachusetts, United States).

Section 2 describes the features of the experimental arrangement for TRESPI. Section 3 presents necessary mathematical derivations. Section 4 discusses the optimization algorithm for the temporal deformation recovery in the case of a harmonically excited object. Section 5 shows how the total deformation can be recovered. Section 6 presents the results of the deformation reconstruction for several excitation conditions, including excitation by hand. Section 7 discusses possible limitations of the system. Section 8 discusses the performance of the methodology, possible extensions and applications of the developed method.

2. EXPERIMENTAL ARRANGEMENT

For the experiments, a diode-pumped solid-state laser Coherent Verdi V5 with the main frequency $\lambda =532\mathrm{nm}$ and maximum output power of 5 W was utilized. To achieve high signal with a short exposure time, the laser power was set to 1.5 W. Recording was performed with a CMOS high-speed camera Phantom v12.1 in a 12-bit mode. The pixel size of the camera sensor is 20 μm.

Figure 1 shows the experimental arrangement for the TRESPI measurement. A polarizing beam splitter with a split ratio 50/50 divides a continuous-wave laser beam into two beams. One of the beams is expanded by a microscopic objective and illuminates the object. The second beam is expanded onto a piece of ground glass and interferes with the light coming from the object after reflection from a beam splitter placed in front of a camera. A high-speed camera records the interferometric images at a frame rate of 21 K frames per second at a resolution of 512 by 512 pixels with exposure time 47 μs. The small angle between the illumination and the observation direction was achieved by placing the microscopic objective as close to the camera viewing angle as possible. In general, the maximum frame rate depends on the resolution of the camera. In TAESPI, the decorrelation of the reference beam is necessary to reduce the random speckle and to decrease the influence of the low-frequency environmental noise [19]. This can be achieved by modulating the reference beam using a piezoelectrically driven mirror. During the high-speed measurements, the interference pattern is resolved in time, so any decorrelation of consecutive frames will add to the relative deformation of the surface. Thus, the only differences between the TAESPI and TRESPI experimental arrangements are that the reference beam is not decorrelated during the recording and the image is captured on a time scale that is short compared to the characteristic time scale of the motion of the object.

 

Fig. 1. Experimental arrangement for a TRESPI measurement. BS, beam splitter; PBS, polarizing beam splitter; M, mirror; L, lens (microscopic objective); $\lambda /2$, half-wave plate; G, ground glass.

Download Full Size | PDF

The object chosen for the experiments described here is a Brazilian cuíca drum from Meinl Percussion, shown in Fig. 2. The membrane is 15 cm in diameter and is made of goat skin. A wooden stick is attached to the center of the clamped membrane. The excitation is accomplished using an electrodynamic shaker (Brüel and Klaer, 4810) or by rubbing the wooden stick by hand. The distance between the object and the camera was approximately 60 cm.

 

Fig. 2. Brazilian cuíca.

Download Full Size | PDF

3. THEORY

This section describes the theoretical aspects of the phase deformation recovery. Some of the details have been covered in the first part of the paper [19].

A. Intensity Normalization

The interference of the light reflected from the object and the reference beam from the same source results in the intensity at a time $t$ being given by

 

Fig. 3. Recorded intensity of a single pixel during the deformation of the surface with sample rate of 21,000 frames per second. Intensity is presented in analog-to-digital units (ADU).

Download Full Size | PDF

In TAESPI the exposure time is much longer than the period of vibration of the object, and the random phase is eliminated by temporal, spatial, and frequency filtering [4]. The recording of TRESPI images is a time-resolved process, hence the intensity is affected by noise. The noise consists of both external environmental noise and internal camera noise. Unfortunately, any filtering to remove the noise affects the information used to calculate the displacement and therefore is not applied.

The first step of the image processing is the normalization of the intensity to remove coefficients $A$ and $B$ for every pixel. In the case of harmonic excitation (under the assumption of a linear system), the coefficients $A$ and $B$ can be considered constant in time during the image acquisition and are given by

However, if the total deformation during the entire oscillation is less than half of the wavelength of illuminating light $\lambda /2$, which corresponds to the phase shift $2\pi$, the normalization and detection of the coefficients $A$ and $B$ cannot be performed with Eqs. (2) and (3). In this case the periodicity in the intensity cannot be tracked, and the upper and lower values of the temporal intensity do not represent the pixel modulation. Therefore, measurements from regions of such a low deformation are not considered to be valid and the modulation value is assigned to the noise level. During the image processing, the deformation values at those invalid pixels are interpolated. One of the approaches to this issue is to introduce a temporal carrier in the first milliseconds of the recording time by inducing a phase shift in the reference beam that is greater than $2\pi$ [18]. Alternatively, an EMD could be applied [13–15], but this would change the phase recovery procedure and therefore it is not applied in the current method.

If the frame rate is large enough so that the deformation can be considered linear in a certain time interval, the intensity on a single pixel is a cosine oscillation in time with varying modulation, which comes from the speckle drift. Figure 3(a) shows that for harmonic oscillation the background intensity and modulation can be regarded constant for a small period of time (here around 5 ms, around two periods for a given frequency). When the phase of the deformation exceeds $2\pi$ during the recording time, the coefficients $A$ and $B$ can be calculated using Eqs. (2) and (3). However, due to the presence of camera noise it is more convenient to use average peak extrema instead of global extrema. This is achieved by calculating a histogram of the intensity values for a given time interval and assigning the lower and upper bin histogram values to the maximum and minimum intensity [23]. With given assumptions, the analysis can be applied to transients as well. Applying surface optimization would reduce inference of camera noise and when the proper boundary conditions are set, temporal deformation can be recovered for the whole surface.

Figure 4 shows recorded pixel intensity in time for three neighboring pixels under harmonic excitation at 627 Hz. Average intensity and modulation vary from pixel to pixel, but the temporal profile shows common repetitive behavior. Intensity peaks are shifted with respect to each other due to random initial phase. However, the effect of camera noise is smaller in the middle graph due to higher modulation, and the information from this particular pixel would be more valuable than from the others.

 

Fig. 4. Pixel intensity from three neighboring pixels chosen in the middle part of the drum harmonically excited at 627 Hz. Intensity is measured in ADU.

Download Full Size | PDF

B. Pixel Quality

The quality of the signal from every pixel is defined by the matrix $Q$ given by

The camera noise is assumed to be the standard deviation of the amplitude of the background in the image sample over time. If $Q$ is less than 1 for a given pixel, the pixel is considered to be not valid and is discarded. The camera noise was calculated at a 20-by-20 pixel region in the background near the drum area. An example of the matrix $Q$ is shown in Fig. 5. To increase the visibility during the measurements the object has been painted matte white.

 

Fig. 5. Pixel quality matrix $Q$.

Download Full Size | PDF

C. Calculation of the Phase Difference

The normalized intensity $I_j$ of a single image recorded at a time $t_j$ for every pixel is

The noisy speckled structure comes from different values of the random speckle phase for every pixel. The intensity is additionally affected by camera noise but pixel-by-pixel temporal filtering (averaging, Golay filter, etc.) reduces peak heights and affects the normalizing procedure [19], so no temporal filtering is applied.

The first frame of the recorded sequence is the starting point of the calculation of the deformation. The speckle phase is defined as ${\varphi }_0$ and the corresponding intensity is $I_0$. The phase recovery can start at any frame and move forward or backward in time. The modulus of the speckle phase can be evaluated as the inverse cosine of the first frame,

The modulus of the random phase is bounded by 0 and $\pi$, which limits the recovered phase to only positive values. At the first glance, if the deformation direction during the time between the first two frames is constant and the phase difference is less than $2\pi$, the random speckle phase can be estimated as

However, if the absolute deformation at a particular pixel exceeds twice the initial speckle phase, Eq. (8) is not valid. The equation holds only for cases when $\mathrm{\Delta }\varphi$ and ${\varphi }_0$ have the same sign and their absolute values are lower than $\pi /2$. Possible outcomes are presented in Table 1. More detailed analysis can be found in [23].

Table 1. Conditions on Deformation Phase and Relation Between the Sign of the Initial Speckle Phase and the Intensity Difference

View Table

The qualitative assessment of the deformation shape at a certain moment in time is accomplished by subtracting the image of interest from the first image. The intensity difference between the first and the $j$th image is

In the subtracted image black fringes are fringes of equal displacement and represent a deformation of $\lambda /2$, corresponding to a $2\pi$ phase difference [Eq. (6)]. Figure 6 shows the result of such a subtraction of images of the cuíca drum excited harmonically at 627 Hz by a shaker. The black lines can be regarded as isolines, i.e., lines of equal displacement. This image is used to verify the recovered phase—the black lines in the experimentally obtained image should coincide with those in the reconstructed image.

 

Fig. 6. Normalized ESPI image of the cuíca during harmonic excitation at 627 Hz. Subtraction intensity corresponds to a time difference of 1.2 ms (25 frames).

Download Full Size | PDF

It is not always possible to predict the phase difference from only two frames, because the sign of ${\varphi }_0$ is unknown. The $2\pi$ phase difference can be mistaken for a nodal line and vice versa, and the black lines in the images are noisy. In this case direct recovery of the phase will always result in a non-zero amplitude. Assessment of a series of frames helps to resolve the $2\pi$ ambiguity. The algorithm for deformation recovery between two subsequent frames is described in the following section, followed by a discussion of the total phase recovery.

4. SINGLE-PHASE DIFFERENCE RECOVERY

This section describes the algorithm used for the recovery of a single-phase difference. The problem is ill-posed since the direction of the deformation cannot be obtained from the interferometric image. However, it may be determined by considering the mechanical properties of the object. The speckle noise and internal camera noise create additional problems.

A. Surface Representation

The surface under study is represented using an adapted grid based on the concept of quadtree division [24]. The implementation of such an approach is discussed in detail in the first part of this work [19]. The surface deformation is represented by a grid which is adapted during the optimization procedure (a feature that was not required in the work reported in [19]) so it becomes more dense in the areas where more precision is required. The grid point values are the optimization parameters and the surface map is produced by bilinear interpolation of the deformation values at the grid points.

B. Adaptive Time Step

During the motion some parts of the surface may deform much more slowly than the rest of the object, so the intensity difference in these regions between adjacent images may not exceed the camera noise and the recovery procedure will fail. If one part of the object deforms noticeably after one recorded frame while the other part moves insignificantly, the intensity variation of part of the image is treated as noise and it is not possible to determine the deformation. This means that more images are required to analyze certain regions of the object. In this case, the surface can be divided into regions of similar “deformation speed” and the frame after which the intensity difference from the first frame is greater than the camera noise (or exceeds a certain signal-to-noise ratio) is determined for every pixel. The matrix of the arguments is denoised, spatially averaged, and clustered. An example of such an output matrix is shown in Fig. 7. The objective function $F_{\text{obj}}$ is found for every region separately and is summed afterwards; therefore, different frames participate in the formation of the mismatch error and the objective function.

 

Fig. 7. (a) Matrix with frame numbers at which the intensity difference exceeds the camera noise, and (b) matrix of the segmented regions with indicated frame numbers. Dark blue is the background which is cut off by masking.

Download Full Size | PDF

In the example shown in Fig. 7(b) the lowest value is 2, which means that the deformation between frames 1 and 3 is recovered over the entire object. The deformation is considered linear between frames 1 and 3 and the next optimization procedure starts at the frame 3.

C. Optimization

To understand the optimization procedure it is useful to consider three consecutive frames $I_0,I_1,I_2$. Assuming a linear out-of-plane displacement of the surface during the time required to capture the three images, the intensities are described by

The Rosenbrock optimization compares the simulated intensity patterns with the measured intensities $I_1$ and $I_2$. The objective function contains two error functions $E_1,E_2$, the pixel quality matrix $Q$, and boundary conditions $C$ and is given by

The matrix $C$ is calculated for every pixel and contains boundary conditions that assume the phase difference is bounded from the upper and lower sides by given values $U$, $L$ such that

The total objective function is a sum over all the pixel values normalized to the total number of pixels:

The initial guess for the optimization is the flat phase map with four nodes on the edges. Thus, the starting direction of the Rosenbrock search is described by the 4-by-4 unit matrix (4 nodes) and all the nodes have the same weight in the optimization algorithm. The values at the nodes are changed according to the Rosenbrock algorithm until every point has at least one successful and one unsuccessful optimization step. There is no basis rotation as in the original Rosenbrock method to facilitate the computation. Testing has shown that this simplification has no effect on the resulting calculation of the surface deformation. The grid is rearranged and is divided as a quadtree when a better spatial resolution is necessary, i.e., when the local objective function is higher than $F_{\mathrm{obj}}^{\mathrm{\Sigma }}/e$ ($e\sim 2.71$).

After the quadtree division, the number of parameters of the optimization changes. The algorithm can be generalized and applied to any number of frames with the assumption of a linear deformation in time. In the case where multiple frames are considered, the standard deviation is used as the error function instead of the difference between $E_1$ and $E_2$. Taking a series of frames helps to recover nodal lines and speeds up the process. The result of the single optimization ($\mathrm{\Delta }{\varphi }_{\mathrm{1,2}}$) is shown in Fig. 8.

 

Fig. 8. Optimized phase difference between two frames of the recorded sequence. The optimization procedure recovers only positive values of the phase.

Download Full Size | PDF

D. Phase Sign Detection

The optimization algorithm produces only positive deformation values; therefore, a sign detection procedure is necessary. The proposed method for the sign detection is based on the assumption that the phase of regions separated by a nodal line will have different signs. The sign detection procedure is fully presented in the first part of this work [19]. The nodal lines are detected by thresholding and skeletonization, the nearby regions are separated, and all possible combinations of the positive and negative phase values are considered. The criteria for the mechanically stable solution is that the sum of the pixel values near the nodal lines should be minimal since the nodal line should separate regions that move out-of-phase with respect to each other. Naturally, two possible complementary solutions will appear, so the total resulting phase may be inverted. The ambiguity can be solved by placing a laser Doppler vibrometer behind the object to determine the correct absolute phase at a single point.

After the sign of the phase is corrected for the first recovered phase difference, ${\varphi }_0$ is found using Eq. (8). The result of the deformation recovery between two frames is shown in Fig. 9 together with the resulting experimental and calculated intensities. In Fig. 9(c) the lines of zero displacement are clearly defined and do not contain speckle noise. This method allows the deformation recovery between consecutive frames with the correct location of the nodal lines from a series of interferometric images with an assumption of linear displacement between frames.

 

Fig. 9. (a) Result of the phase sign correction, (b) experimental subtraction intensity, and (c) calculated intensity.

Download Full Size | PDF

5. TOTAL PHASE RECOVERY

During the optimization the grid is adapted spatially during the calculation of the phase difference between two moments in time. Temporal step is calculated before the spatial optimization procedure and can be updated after each optimization step.

When all of the phase differences are found, it is possible to determine the total phase between the first and the $j$th frame by simple summation of the phases:

However, direct summation has several drawbacks. First, it may result in a cumulative error, especially at the boundaries where small values can add up to violate the fixed boundary condition. Second, the sign detection algorithm may fail in the regions where the deformation changes the sign. Finally, the addition of any incorrect phase difference will result in an incorrect value of the final phase.

An alternative method to construct the total phase is to use the initial phase value that is found during the first phase difference recovery with a proper sign. Instead of using Eq. (8) to find ${\varphi }_1$ from the intensities, ${\varphi }_1$ is composed of the previously defined values:

Thus, the next phase difference is added to the previous phase and becomes a starting phase for the next series of images. This procedure provides the certainty in the sign of the phase as long as the first frame is correctly determined. If the phase of the first frame is correct, the following frames do not need an additional sign detection. After the phase differences are obtained for every frame, they are summed to determine the absolute value of the phase. This method requires an assumption that coefficients $A$ and $B$ are considered constant for a particular time interval. The coefficients can be recalculated for higher precision at any given time step. Then the normalized intensity and the phase sign must be updated, and a new ${\varphi }_1$ must be established.

6. RESULTS

In this section results of the experiments performed under two types of excitation are presented. The drum has been either excited by a shaker or played by hand. During harmonic excitation, the electrodynamic shaker was attached to the wooden stick behind the drum membrane. The temporal resolution of the high-speed camera was 50 μs. Figures 10 and 11 show the results of the temporal displacement recovery at two frequencies. The time between the selected frames is approximately 0.2 ms. For analysis of the surface deformation under playing conditions, the drum was excited by slightly rubbing the wooden stick that is attached to the center of the membrane by hand. Total relative deformation is shown in Fig. 12 with a 5 ms interval between images.

 

Fig. 10. Results of temporal reconstruction during the harmonic excitation of the cuíca drum at 627 Hz with interval of 0.2 ms.

Download Full Size | PDF

 

Fig. 11. Results of temporal reconstruction during the harmonic excitation of the cuíca drum at 787 Hz with interval of 0.2 ms.

Download Full Size | PDF

 

Fig. 12. Results of temporal reconstruction during the excitation of the cuíca drum by hand with interval of 5 ms.

Download Full Size | PDF

In Figs. 10–12 the fringe patterns $(c)$ are directly calculated from the deformations $(a)$. The fringe pattern $I_{\mathrm{calc}}$ is calculated at each pixel as

When the optimization method is applied to noisy interferograms there is a large error in the reconstructed image. The error appears to be larger in the areas of low surface velocity, where the camera noise exceeds the intensity difference caused by motion. However, the areas of equal displacement are correctly defined and the calculated intensity patterns show good agreement with experimentally obtained interferograms.

In the current experiments the phase could be retrieved without recalculating the initial speckle phase for 20 frames which corresponds to 1 ms, and up to 20 ms in the case of a very small deformation during a slight rubbing of the wooden stick of the drum.

Figure 13 shows a block diagram with the main steps of the image processing routine for the total phase l. The step “intensity normalization” is described in Section 3, the steps preceding the iterative phase recovery process refer to Section 4, and “total phase” calculation follows from Section 5.

 

Fig. 13. Block diagram describing the main steps of the processing for the temporal displacement retrieval.

Download Full Size | PDF

Optimization stops when it hits the plateau of the objective function or the surface grid division cannot go lower than the particular size (in this case it is 32 or 64 pixels). The number of iterations (e.g., calculation of the optimization function and following parameter update) is limited to 5000.

Single-phase difference processing takes 10–15 min on 512 by 512 pixels images on a standard laptop computer (Intel Core i7-7500U with 8 Gb RAM). The processed datasets consist of 20–60 images. When the deformation is very low, frames for processing can be taken with a larger temporal interval between them (in the case of hand excitation, every fifth frame was taken).

7. DISCUSSION

As in every experimental method, certain limitations can be identified. The main limiting factor for the proposed deformation retrieval method is the surface velocity with respect to the camera frame rate. At large displacements speckles can drift in space so that the initial phase value is no longer valid and the objective function increases. Another limitation could be the uncertainty in definition of the phase sign, which is obtained using assumptions on the mechanical properties of the object under deformation. The deformation between frames is considered linear in a sense that the direction of the deformation does not change between two exposures. The phase differences are recovered at discrete times, and temporal deformation of a single pixel is represented by a function that is not inherently smooth. This can be corrected by a polynomial fitting of the deformation in time, a process that may also increase the performance of the optimization algorithm. Furthermore, it may be required to paint an object to remove specular reflections and provide an even reflection of the laser (this can be performed with dissolving “snow frost ”).

The limitation of the linearity between more than two consequent frames may be resolved by acquiring a camera with a higher frame rate. In the high-speed camera used in the experiments, the frame rate was related to its spatial resolution. The limit of the 21 K frames per second comes from the fact that the resolution of the image was kept 512 by 512 pixels. Nevertheless, the proposed TRESPI method appears to be an extension of time-averaged ESPI with opportunities for further enhancement.

8. CONCLUSIONS

The results reported here show that the optimization method can be applied in TRESPI for out-of-plane deformation analysis with both harmonic and non-harmonic excitation. The phase recovery procedure uses a flexible time step to increase the precision of the optimization and quadtree spatial grid division to reduce the number of optimization parameters. Future improvements involve full automation of the deformation recovery procedure, application to high-amplitude cases (more than $2\pi$ phase difference between the frames), reduction of the cumulative error, and extensional error analysis.