1. INTRODUCTION

Nondestructive optical testing techniques are commonly used for high-precision measurements of surface deformations [1]. These methods are primarily based on interferometric measurement techniques, which have evolved significantly since the invention of the first laser [2]. One such method is electronic speckle pattern interferometry (ESPI), which can be used to visualize both static and dynamic displacements of an object in time [3].

ESPI encompasses a class of full-field experimental techniques that use intensity information from an object to detect surface deformations using properties of the optical speckle. The speckles appear due to the interference of laser light that is diffusely reflected from the rough surfaces. For every point in space, reflected wavefronts, each with a different amplitude and phase, sum resulting in a random intensity distribution on the surface of the detector. The information about the surface displacement between two moments in time is encoded in the intensity pattern in the form of correlation fringes, which appear due to the change in the intensity pattern of the speckle that occurs when the object is deformed. These temporal intensity images are processed to determine the absolute displacement from the phase difference between two interfering beams, one reflected from the object under study and the other serving as a reference beam. The displacement and the phase difference are related linearly. Thus, knowledge of the phase difference automatically provides information about the displacement [4]. Therefore, the main objective of the ESPI is phase recovery. Depending on the illumination arrangement, ESPI can detect in-plane [5–7] or out-of-plane [8,9] deformations. An arrangement designed to detect out-of-plane displacement is insensitive to an in-plane motion and vice versa. These two methods can be combined to build a three-dimensional measurement system [10,11].

ESPI is usually applied for vibrational analysis since the experimental arrangement is not complicated and uses a conventional digital camera [12]. The method is referred to as time-averaged ESPI (TAESPI) because the exposure time of the camera is much longer than the period of vibration of the object. Moreover, the correlation fringes can be observed in real time. Therefore, the out-of-plane configuration is widely used for the investigation of the operating deflection shapes of vibrating objects. For more complicated deflection shapes, a high-speed camera is necessary [13]. In some cases, the correlation fringe patterns give enough qualitative information, and the only task is to filter the images to eliminate the speckle [14]. However, a measure of the deformation is of interest to many researchers [15], and this is not a trivial task due to the phase ambiguity, poor temporal and spatial resolution, and external noise.

Most ESPI processing methods are based on the spatial transformation of the recorded data, so each pixel is processed independently to obtain a phase map. If the phase difference between two recordings exceeds $2\pi$, the phase map becomes spatially periodic (i.e., wrapped). Phase unwrapping is then required to recover the true phase values from a wrapped phase distribution. This is generally accomplished by an iterative assessment of local discontinuities of the phase maps, obtained from the recorded intensity frames [16].

Another method using global multiobjective optimization was proposed recently by the authors [17]. The surface is discretized with a given spatial step and the corresponding intensity values are simulated. Operating on a grid was chosen in order to minimize the computing time. The idea behind the solution of the inverse problem is to predict the amplitude values that minimize the mismatch between experimentally obtained intensity and simulated patterns. The surface parameters are changed smoothly to decrease the difference between the experimental and the simulated intensity at each point. The optimization is based on the Rosenbrock algorithm. It is a derivative-free direct search method that approximates a gradient search [18]. However, using the conflicting objective functions does not lead to a fully meaningful solution. The weighting of the multiobjective parameters must be adjusted for every case, so the multiobjective optimization may fail for complicated deflection shapes. Additionally, not all the pixels are taken into account due to the long computing time, so when the grid is coarse, some surface features may be missed. Nevertheless, the method may be improved to overcome these obstacles.

The goal of this work is to develop a simple experimental method and a robust image processing routine that is applicable to both steady-state and transient deformation analysis without using phase-shifting algorithms and any devices that would require exceptional stability of the laboratory arrangement. Taking into account physical properties of the investigated object may also benefit the image processing—knowledge of the boundary conditions, maximum attainable amplitude, and the excitation point are useful for meaningful displacement amplitude retrieval.

The aim of this work is to recover the out-of-plane surface deformation profile of the harmonic vibration from the recorded interferometric images. This first part of a two-part investigation is focused on TAESPI and the image processing methods to retrieve the out-of-plane surface displacement for the case of harmonic excitation. The theoretical section is adopted from previous work [17], but the processing has been significantly improved. The current method can be applied to more cases and the relative phase is also retrieved. This method does not require an unwrapping procedure in its classical sense because the local amplitude range is assigned before the phase retrieval and is not ambiguous. Image processing is accomplished using MATLAB software (MATLAB and Image Processing Toolbox R2014a, The MathWorks, Inc., Natick, Massachusetts, USA).

Section 2 summarizes the previous work on the phase reconstruction methodology for time-averaged speckle pattern interferometry, including the experimental arrangement and relevant formulas. Section 3 discusses the modified optimization algorithm for the recovery of the amplitude of the deformation phase for the case of a harmonically excited object. Section 4 presents the results of the deformation reconstruction for different excitation frequencies. Section 5 concludes the paper and notes the possible extension of the methodology to 3D imaging and transient analysis.

2. PRELIMINARIES

A. Experimental Arrangement

The experimental arrangement for the out-of-plane displacement measurements described in detail in the preceding papers [19,20] is shown in Fig. 1. A continuous wave laser beam with a wavelength $\lambda =532\mathrm{nm}$ is split in two beams, one illuminating the object and the other serving as a reference beam. A CCD camera records the intensity of the two interfering beams at a rate of 10 frames per second. A piezoelectrically driven mirror (PZT-mirror), producing a linear temporal phase shift in the reference arm, is synchronized with the camera. The slow decorrelation of the reference beam enhances the sensitivity compared to the arrangement without a moving mirror [21]. The reference beam is projected onto a ground glass, which is imaged through a beam splitter. The object for the investigation is a Brazilian cuíca drum. It is fixed on the optical table and is harmonically excited by an electromechanical shaker at 757 Hz.

 

Fig. 1. Experimental arrangement for a TAESPI measurement. A polarizing beam splitter (PBS) divides the laser light into an objective and a reference beam. The first half-wave plate ($\lambda /2$) controls the relative intensity of the beams, and another $\lambda /2$ plate controls the intensity of the reference beam. Lenses ($L$) expand the beams. A PZT-actuated mirror adds a linear phase shift to the reference beam. A ground glass ($G$) randomizes the phase of the reference beam. The two beams are combined by a 50/50 beam splitter (BS) and the intensity is recorded by the camera.

Download Full Size | PDF

B. Theory

If an object linearly responds to a harmonic excitation, the resulting phase difference between the object and the reference beam is described as

The previous work [17] introduced two methods to determine vibration amplitude from a series of ESPI images: peak direct inverse method and the optimization method. The peak direct inverse method determines local extrema of the fringes in the ESPI images and associates the arguments of the Bessel function at maxima and minima. The optimization method is applied on a grid of values and takes advantage of the continuity of the object surface. It calculates an objective optimization function that must be minimized upon a change of the pixel values. The optimization method was shown to be more efficient than the peak direct inverse method in terms of spatial resolution and performance time. In the current manuscript, the image processing routine has been modified to assess several different deflection shapes, incorporate boundary conditions, and determine relative displacement of the in-phase and out-of-phase regions of the surface.

C. Filtering

The filtering procedure involves temporal, spatial, and frequency filtering to remove speckle effects from the interferometric images (see [17] for details). First, the TAESPI images of the vibrating object and the reference frames of the nondeformed object are temporally averaged over all the recorded frames. Therefore, the image processing only includes two frames: the deformation frame and the reference frame. This temporal averaging decreases the effect of the random phase of the speckle. Usually, 30–50 frames are recorded for each case so that the recording time does not exceed 10 s, including the time necessary to store the images. Both resulting images are spatially filtered using averaging and median filtering on a $3\times 3$ pixel neighborhood. Then the intensity of the image of the vibrating state is divided by the reference image to remove the coefficient $\beta M$ for every pixel, resulting in an image with normalized intensity. Localized Fourier filtering (LFF) is then applied to reduce speckle noise, which involves Gaussian smoothing of the entire image and windowed Fourier filtering [22]. The size of the block is $12\times 12\mathrm{pixels}$ and the filter parameter is 0.1. As in [17], the LFF is applied to a complement image of the normalized intensity and the image is inverted after the filtering procedure and then normalized.

D. Bessel Fringe Orders

Direct inversion of the modulus of the quasi-periodic Bessel function results in an ambiguous phase detection. Therefore, every pixel should be limited to a possible range of values to resolve the uncertainty. The filtered image is thresholded to detect and separate regions near the roots of the modulus of the Bessel function (dark fringes) and regions near the peaks (bright fringes). These fringes are independently processed and the phase intervals are assigned for every fringe. For the bright fringes, the algorithm searches for the peak values close to the values of the peaks in the $|J_0|$ function. The dark fringes are then assessed depending on the neighboring bright fringes. Each pixel is assigned a value that corresponds to a phase interval depending on whether it belongs to a bright or a dark fringe (see [17] for details). This procedure removes the uncertainty in phase reconstruction, resulting in a map of the Bessel fringe orders (not to be confused with the orders of the Bessel function), which shows the discrete bands corresponding to the accepted values of the phase. Each bright fringe is normalized on the peak value of the Bessel region between the two closest roots to facilitate phase recovery. An experimentally obtained single TAESPI image and the resulting Bessel fringe map are shown in Fig. 2. The fringe map gives a qualitative assessment of the vibration amplitude since every level reflects the corresponding phase interval. The Bessel fringe orders correspond to positive values of the phase. Thus, it is not possible to differentiate between in-phase and out-of-phase regions.

 

Fig. 2. (a) Experimental TAESPI image after filtering, and (b) assigned Bessel fringe orders. The cuíca drum is harmonically excited by the shaker at 757 Hz.

Download Full Size | PDF

3. DEFORMATION RECOVERY

The optimization method is based on the solution of the inverse problem by the minimization of an objective function during gradual deformation of the simulated displacement surface using the Rosenbrock algorithm. At every step of the optimization, the phase is changed according to the Rosenbrock approach, and the intensity is simulated by applying the modulus of the Bessel function. The objective function is formed from the total pixel-by-pixel difference between the simulated intensity pattern and the experimentally obtained normalized intensity, with specified boundary conditions.

A. Surface Representation

Using the phase values of all the pixels as the variable parameters is not advisable due to the large amount of memory and time required for computation. It is desirable to have a grid that adapts to the scale of the deformation amplitude of the surface so that it is dense in the areas where the gradient of the deformation amplitude is high. On the other hand, almost flat surfaces can be described using a small number of reference points, which can be processed quickly. Thus, the surface is spatially discretized using an approach based on the quadtree division [23], which is a tree data structure used to separate the regions by recursive subdivision in the spatial domain. It is generally used in image analysis to iteratively divide an image into smaller equal-sized blocks so that all blocks meet a specific homogeneity criterion.

Using the quadtree division scheme, the whole surface area is represented as square blocks with nodes on the corners. If an objective function taken from a single block does not meet a specific criterion, the block is separated into four. Thus, the roughness of the object surface and the precision of the optimization procedure are governed by the distance between the nodes. An example of a single block and the block after two divisions is shown in Fig. 3.

 

Fig. 3. Operating points on the grid.

Download Full Size | PDF

The phase amplitude values at the nodes are the parameters of the optimization that change during the processing. The amplitude of the deformation phase $\mathrm{\Phi }$ is constructed by bilinear interpolation of the values on the grid points. For every square block, the pixel value inside the block is estimated as the distance-weighted average with the respect to the four corner pixels [24]. However, the quadtree grid is not homogeneous and if the method is applied directly, it results in discontinuities on the edges of the blocks. Therefore, the grid size is at first scaled down to the smallest distance between the nodes obtained by the quadtree division, and the values on the new nodes are found by bilinear interpolation. Subsequently, the phase surface is obtained, again using bilinear interpolation. An example of the grid values and the resulting phase based on the grid from Fig. 3(b) is shown in Fig. 4.

 

Fig. 4. (a) Quadtree division points from Fig. 3(b), (b) the modified grid, and (c) the interpolated phase.

Download Full Size | PDF

The quadtree division reduces the total number of optimization parameters and significantly speeds the recovery algorithm. Physically, this division corresponds to the gradual change of the stiffness and number of degrees of freedom of the discretized object.

B. Optimization

The optimization routine changes the values on the nodes of the grid, produces values of $\mathrm{\Phi }$ over the whole surface in the form of a matrix, and compares the calculated intensity pattern $|J_0(\mathrm{\Phi })|$ with the normalized experimental intensity $I_{\text{norm}}=|J_0(\xi )|$ to form the error function $E$:

The total objective function is a sum of $F_{\mathrm{obj}}$ over all pixel values of a $k\times l$ image normalized on the total number of pixels:

The initial guess for the optimization is the phase map with the values from the middle interval points of the solutions of the Bessel inverse function, based on the matrix of the assigned Bessel fringe orders. Thus, the correct Bessel fringe order assignment is the most important part of the algorithm. If it fails, the method will optimize the phase values to the wrong local solutions. The initial phase guess is shown in Fig. 5.

 

Fig. 5. Initial phase amplitude guess for the optimization.

Download Full Size | PDF

The quadtree-based grid has to be adjusted to the initial phase so that the boundary condition $C$ is met. Prior to the Rosenbrock optimization, the blocks are iteratively divided in order to meet the condition

 

Fig. 6. Quadtree-based grid adapted to the scale of the initial guess of the deformation amplitude. Black points indicate the nodes.

Download Full Size | PDF

For $n$ grid points, the starting direction of the Rosenbrock search is described by an $n$-dimensional unit matrix. The matrix consists of the orthonormal basis vectors, and all of the nodes have the same weight. The values at the nodes are changed according to the Rosenbrock algorithm with respect to the basis vectors until every point has at least one successful and one unsuccessful optimization step. Thereafter, the matrix of the basis vectors is rotated, changing the influence of every node in a way that the most successful direction becomes the first vector for the new basis. This operation is defined as the “basis rotation” [18]. The rest of the basis is built orthogonally to the first vector, and the optimization procedure continues with the new basis. The criterion for stopping the optimization is either exceeding the maximum number of basis rotations or the maximum number of iteration steps assigned by the practitioner.

In the case of a simple deformation profile, such as tilt, the phase amplitude can be represented as a polynomial model. This procedure drastically reduces the number of parameters, but it can be applied in only a few situations.

C. Phase Sign Detection

The nodal lines are detected by MATLAB built-in morphological operations, such as thresholding and skeletonization (Image Processing Toolbox). The nearby regions are separated and all possible solutions of the positive and negative phase values are produced. The criterion for the mechanically stable solution is that the sum of the pixel values near the nodal lines should be the minimum of all possible combinations for in-phase and out-of-phase motion of neighboring areas. This approach does not necessarily hold for transient events, but it works well for steady-state vibrational analysis. The procedure is performed only once after the optimization, so it does not significantly affect performance.

The absolute phase resulting from the optimization procedure is shown in Fig. 7(a). The skeleton shown in Fig. 7(b) separates two areas within the object. Changing the signs of these two areas independently leads to four possible solutions, two of which are symmetric with respect to the undeformed object. The parts with opposite sign are chosen and shown in Fig. 7(c). Even though the nodal lines are not clearly defined, small branches that appear from the image processing do not affect the phase sign recovery procedure. Figure 7(d) shows the interpolated result, which is smoothed using a least-square smoothing algorithm [25] to reflect the surface continuity. A block diagram of the process is shown in Fig. 8.

 

Fig. 7. (a) Optimized absolute phase amplitude distribution, (b) the skeleton showing the nodal lines, (c) corrected separated phase regions, and (d) the resulting phase distribution.

Download Full Size | PDF

 

Fig. 8. Block-diagram describing the main steps of the experimental procedure and image processing for the displacement retrieval.

Download Full Size | PDF

4. RESULTS

The TAESPI image, the reconstructed correlation image, and the out-of-plane vibration amplitude $A$ with recovered relative phase for harmonic excitation ${\omega }_0=757\mathrm{Hz}$ are shown in Fig. 9. The smallest grid size was 16 pixels, the number of iterations was 20,000, and the number of basis rotations was 3. The vibration amplitude is obtained using Eq. (2). The experimental and intensity images are in good agreement. The brightest areas correspond to the deformation in the range from 0 to 102 nm, including the nodal lines, which corresponds to the first zero of the Bessel function. The distance between two black fringes corresponds to a difference of approximately 131 nm. The vibration can be temporally reconstructed as $A\sin ({\omega }_{0}t)$.

 

Fig. 9. (a) Experimental and (b) reconstructed TAESPI images for the excitation frequency 757 Hz and (c) the out-of-plane deformation amplitude.

Download Full Size | PDF

The comparison between the intensities of a line in the middle of the image is shown in Fig. 10. It is possible to calculate the average pixel error between the normalized experimental intensity and the calculated intensity. However, the black fringes ideally should pass through zero, which does not happen in the experimental image due to the presence of noise. Therefore, the pixel error does not provide a relevant evaluation of the reconstruction method. The procedure of matching the black fringes from the calculated intensity to one from the experiment is more beneficial. The profiles are very similar and the peak positions coincide. The recovered profile is smooth and noiseless, which is useful for the visualization of the interferometric data.

 

Fig. 10. Profiles of the experimentally obtained normalized interferometric image and reconstructed interferometric pattern, taken as the absolute value of the zero-order Bessel function of the phase map.

Download Full Size | PDF

A. Phase Recovery for Various Frequencies

This section presents results on the recovered vibration amplitude for several applied excitation frequencies (Fig. 11). The first column shows the experimentally obtained TAESPI image, the second column is the calculated intensity, the third column is a 3D representation of the out-of-plane displacement. The smallest grid size varies between 16 and 64 pixels, depending on the intensity profile, and the optimization runs for a maximum 20,000 iteration steps as in the previous case. The operating deflection shapes become more complicated with increasing excitation frequency. However, this does not adversely affect the analysis. The nodal lines are clearer in the simulated patterns and the deformation amplitude reflects the out-of-plane vibrational behavior of the surface. These results show robustness for various intensity distributions, independent of the fringe orientation. The results shown in Fig. 11 were obtained using the same parameters (threshold values for the Bessel fringe order assignment, filtering parameters, threshold values for the nodal line detection, etc.) for all frequencies, and the relative phase of the displacement is computed automatically after the amplitude recovery.

 

Fig. 11. Results for four different frequencies of excitation of the cuíca drum.

Download Full Size | PDF

The validity of the optimization method has been previously assessed using the laser Doppler vibrometry technique [20]. Single point measurements were performed at several positions on the drum membrane, and the results were compared to the recovered values of the deformation amplitude. The error between the measured displacement and the reconstruction is around $\lambda /15$.

The random phase is a part of the $M$ coefficient, which is removed during the normalization procedure described in Section 2. The coefficient $\beta M$ can be estimated for every frame through the division of the subtracted recorded frames $I_{\mathrm{sub}}(\xi )$ by the calculated intensity $|J_0(\xi )|$. For every pixel, the standard deviation $\sigma$ of the value $\beta M$ over time shows how much the intensity pattern is affected by noise. The mean value of $\sigma$ over all pixels estimates the average intensity noise error. For the given excitation frequencies, the mean error does not exceed 4% in each case. However, the dark fringe intensity areas were excluded from the calculation. A threshold of 0.05 was applied because the division by zero intensity (the roots of the Bessel function) would produce inaccurate data.

B. Changing the Optimization Parameters

A small grid size and large number of iterations make the optimization performance time-consuming. The process can be accelerated by increasing the minimum distance between the nodes and decreasing the number of the iteration steps. These two parameters are crucial in determining the efficiency of the method.

Figure 12 shows the comparison between the calculated intensity patterns achieved using two different sets of parameters. The first image in Fig. 12(a) depicts the previously obtained result shown in Fig. 11(d). In Fig. 12(b), the grid size is 32 pixels and the number of iterations is decreased to 1500. The amplitude appears to be a bit lower, with the error due to the least-square smoothing and roughness of the quadtree-based grid. In the first case, the simulation matches the experimentally obtained intensity pattern more closely. However, the processing time increases linearly with the number of iterations and becomes the major inconvenience of the method.

 

Fig. 12. Comparison of the calculated correlation images obtained with different optimization parameters.

Download Full Size | PDF

5. CONCLUSIONS

The optimization method described above has been applied in TAESPI for out-of-plane vibration analysis. It is an effective and robust method that requires neither phase-shifting techniques nor sophisticated experimental arrangement. The method does not depend on the fringe orientation and is almost fully automated. The reduction of the deformation amplitude error can be achieved by an increase in the number of optimization parameters, with a commensurate increase in the computational time. Furthermore, the proposed optimization algorithm can be implemented in three-dimensional vibrational analysis in the existing experimental arrangement.

The second part of this work will report on the application of this optimization methodology to transient deformations.